Possibility to manually specify no_random_effects in mhmkin

9 files changed, 187 insertions, 94 deletions

diff --git a/NEWS.md b/NEWS.md
index 35f499d5..65df4191 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,6 @@
 # mkin 1.2.2
 
-- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well.
+- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object or a list with suitable dimensions as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well.
 
 - 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters
 
diff --git a/R/mhmkin.R b/R/mhmkin.R
index 6a61e8bb..a9798ff4 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -14,11 +14,12 @@
 #' supported
 #' @param no_random_effect Default is NULL and will be passed to [saem]. If a
 #' character vector is supplied, it will be passed to all calls to [saem],
-#' regardless if the corresponding parameter is in the model. Alternatively,
-#' an object of class [illparms.mhmkin] can be specified. This has to have
-#' the same dimensions as the return object of the current call. In this way,
-#' ill-defined parameters found in corresponding simpler versions of the
-#' degradation model can be specified.
+#' which will exclude random effects for all matching parameters. Alternatively,
+#' a list of character vectors or an object of class [illparms.mhmkin] can be
+#' specified. They have to have the same dimensions that the return object of
+#' the current call will have, i.e. the number of rows must match the number
+#' of degradation models in the mmkin object(s), and the number of columns must
+#' match the number of error models used in the mmkin object(s).
 #' @param algorithm The algorithm to be used for fitting (currently not used)
 #' @param \dots Further arguments that will be passed to the nonlinear mixed-effects
 #' model fitting function.
@@ -50,6 +51,42 @@ mhmkin.mmkin <- function(objects, ...) {
 
 #' @export
 #' @rdname mhmkin
+#' @examples
+#' \dontrun{
+#' # We start with separate evaluations of all the first six datasets with two
+#' # degradation models and two error models
+#' f_sep_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = 2, quiet = TRUE)
+#' f_sep_tc <- update(f_sep_const, error_model = "tc")
+#' # The mhmkin function sets up hierarchical degradation models aka
+#' # nonlinear mixed-effects models for all four combinations, specifying
+#' # uncorrelated random effects for all degradation parameters
+#' f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cores = 2)
+#' status(f_saem_1)
+#' # The 'illparms' function shows that in all hierarchical fits, at least
+#' # one random effect is ill-defined (the confidence interval for the
+#' # random effect expressed as standard deviation includes zero)
+#' illparms(f_saem_1)
+#' # Therefore we repeat the fits, excluding the ill-defined random effects
+#' f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1))
+#' status(f_saem_2)
+#' illparms(f_saem_2)
+#' # Model comparisons show that FOMC with two-component error is preferable,
+#' # and confirms our reduction of the default parameter model
+#' anova(f_saem_1)
+#' anova(f_saem_2)
+#' # The convergence plot for the selected model looks fine
+#' saemix::plot(f_saem_2[["FOMC", "tc"]]$so, plot.type = "convergence")
+#' # The plot of predictions versus data shows that we have a pretty data-rich
+#' # situation with homogeneous distribution of residuals, because we used the
+#' # same degradation model, error model and parameter distribution model that
+#' # was used in the data generation.
+#' plot(f_saem_2[["FOMC", "tc"]])
+#' # We can specify the same parameter model reductions manually
+#' no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta"))
+#' dim(no_ranef) <- c(2, 2)
+#' f_saem_2m <- update(f_saem_1, no_random_effect = no_ranef)
+#' anova(f_saem_2m)
+#' }
 mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
   no_random_effect = NULL,
   ...,
@@ -97,25 +134,38 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
   dimnames(fit_indices) <- list(degradation = names(deg_models),
                                 error = error_models)
 
-  fit_function <- function(fit_index) {
-    w <- which(fit_indices == fit_index, arr.ind = TRUE)
-    deg_index <- w[1]
-    error_index <- w[2]
-    mmkin_row <- objects[[error_index]][deg_index, ]
+  if (is.null(no_random_effect) || length(dim(no_random_effect)) == 1) {
+    no_ranef <- rep(list(no_random_effect), n.fits)
+    dim(no_ranef) <- dim(fit_indices)
+  } else {
+    if (!identical(dim(no_random_effect), dim(fit_indices))) {
+      stop("Dimensions of argument 'no_random_effect' are not suitable")
+    }
     if (is(no_random_effect, "illparms.mhmkin")) {
-      if (identical(dim(no_random_effect), dim(fit_indices))) {
-        no_ranef_split <- strsplit(no_random_effect[[fit_index]], ", ")
-        no_ranef <- sapply(no_ranef_split, function(x) {
-          gsub("sd\\((.*)\\)", "\\1", x)
+      no_ranef_dim <- dim(no_random_effect)
+      no_ranef <- lapply(no_random_effect, function(x) {
+        no_ranef_split <- strsplit(x, ", ")
+        ret <- sapply(no_ranef_split, function(y) {
+          gsub("sd\\((.*)\\)", "\\1", y)
         })
-      } else {
-        stop("Dimensions of illparms.mhmkin object given as 'no_random_effect' are not suitable")
-      }
+        return(ret)
+      })
+      dim(no_ranef) <- no_ranef_dim
     } else {
       no_ranef <- no_random_effect
     }
+  }
+
+  fit_function <- function(fit_index) {
+    w <- which(fit_indices == fit_index, arr.ind = TRUE)
+    deg_index <- w[1]
+    error_index <- w[2]
+    mmkin_row <- objects[[error_index]][deg_index, ]
     res <- try(do.call(backend_function,
-        args = c(list(mmkin_row), dot_args, list(no_random_effect = no_ranef))))
+        args = c(
+          list(mmkin_row),
+          dot_args,
+          list(no_random_effect = no_ranef[[deg_index, error_index]]))))
     return(res)
   }
 
@@ -145,15 +195,16 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
 #' @param j Column index selecting the fits to specific datasets
 #' @param drop If FALSE, the method always returns an mhmkin object, otherwise
 #'   either a list of fit objects or a single fit object.
-#' @return An object of class \code{\link{mhmkin}}.
+#' @return An object inheriting from \code{\link{mhmkin}}.
 #' @rdname mhmkin
 #' @export
 `[.mhmkin` <- function(x, i, j, ..., drop = FALSE) {
+  original_class <- class(x)
   class(x) <- NULL
   x_sub <- x[i, j, drop = drop]
 
   if (!drop) {
-    class(x_sub) <- "mhmkin"
+    class(x_sub) <- original_class
   }
   return(x_sub)
 }
diff --git a/log/check.log b/log/check.log
index 31fc31eb..42365918 100644
--- a/log/check.log
+++ b/log/check.log
@@ -5,7 +5,7 @@
 * using options ‘--no-tests --as-cran’
 * checking for file ‘mkin/DESCRIPTION’ ... OK
 * checking extension type ... Package
-* this is package ‘mkin’ version ‘1.2.1’
+* this is package ‘mkin’ version ‘1.2.2’
 * package encoding: UTF-8
 * checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
 Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’
@@ -18,7 +18,7 @@ Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’
 * checking for portable file names ... OK
 * checking for sufficient/correct file permissions ... OK
 * checking serialization versions ... OK
-* checking whether package ‘mkin’ can be installed ... OK
+* checking whether package ‘mkin’ can be installed ... [11s/11s] OK
 * checking installed package size ... OK
 * checking package directory ... OK
 * checking for future file timestamps ... OK
@@ -41,7 +41,7 @@ Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’
 * checking S3 generic/method consistency ... OK
 * checking replacement functions ... OK
 * checking foreign function calls ... OK
-* checking R code for possible problems ... [17s/17s] OK
+* checking R code for possible problems ... [19s/19s] OK
 * checking Rd files ... OK
 * checking Rd metadata ... OK
 * checking Rd line widths ... OK
@@ -57,7 +57,7 @@ Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’
 * checking data for ASCII and uncompressed saves ... OK
 * checking installed files from ‘inst/doc’ ... OK
 * checking files in ‘vignettes’ ... OK
-* checking examples ... [20s/20s] OK
+* checking examples ... [24s/24s] OK
 * checking for unstated dependencies in ‘tests’ ... OK
 * checking tests ... SKIPPED
 * checking for unstated dependencies in vignettes ... OK
diff --git a/log/test.log b/log/test.log
index e17ecc1f..84fa49b9 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,58 +1,58 @@
 ℹ Testing mkin
 ✔ | F W S  OK | Context
 ✔ |         5 | AIC calculation
-✔ |         5 | Analytical solutions for coupled models [3.3s]
+✔ |         5 | Analytical solutions for coupled models [4.2s]
 ✔ |         5 | Calculation of Akaike weights
 ✔ |         3 | Export dataset for reading into CAKE
-✔ |        12 | Confidence intervals and p-values [1.1s]
-✔ |     1  12 | Dimethenamid data from 2018 [32.2s]
+✔ |        12 | Confidence intervals and p-values [1.2s]
+✔ |     1  12 | Dimethenamid data from 2018 [42.0s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_dmta.R:98'): Different backends get consistent results for SFO-SFO3+, dimethenamid data
 Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
 ────────────────────────────────────────────────────────────────────────────────
-✔ |        14 | Error model fitting [4.9s]
+✔ |        14 | Error model fitting [6.5s]
 ✔ |         5 | Time step normalisation
-✔ |         4 | Calculation of FOCUS chi2 error levels [0.6s]
-✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8s]
-✔ |         4 | Test fitting the decline of metabolites from their maximum [0.3s]
-✔ |         1 | Fitting the logistic model [0.2s]
-✔ |         8 | Batch fitting and diagnosing hierarchical kinetic models [14.5s]
-✔ |     1  11 | Nonlinear mixed-effects models [13.1s]
+✔ |         4 | Calculation of FOCUS chi2 error levels [0.7s]
+✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.1s]
+✔ |         4 | Test fitting the decline of metabolites from their maximum [0.5s]
+✔ |         1 | Fitting the logistic model [0.3s]
+✔ |        10 | Batch fitting and diagnosing hierarchical kinetic models [54.1s]
+✔ |     1  11 | Nonlinear mixed-effects models [14.3s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_mixed.R:78'): saemix results are reproducible for biphasic fits
 Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ────────────────────────────────────────────────────────────────────────────────
 ✔ |         3 | Test dataset classes mkinds and mkindsg
-✔ |        10 | Special cases of mkinfit calls [0.4s]
-✔ |         3 | mkinfit features [0.7s]
-✔ |         8 | mkinmod model generation and printing [0.2s]
+✔ |        10 | Special cases of mkinfit calls [0.8s]
+✔ |         3 | mkinfit features [0.9s]
+✔ |         8 | mkinmod model generation and printing [0.3s]
 ✔ |         3 | Model predictions with mkinpredict [0.3s]
-✔ |        12 | Multistart method for saem.mmkin models [50.1s]
-✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.2s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [8.7s]
-✔ |        15 | Plotting [10.2s]
+✔ |        12 | Multistart method for saem.mmkin models [80.1s]
+✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.8s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [11.4s]
+✔ |        15 | Plotting [12.1s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |     1  36 | saemix parent models [103.8s]
+✔ |     1  36 | saemix parent models [85.9s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_saemix_parent.R:143'): We can also use mkin solution methods for saem
 Reason: This still takes almost 2.5 minutes although we do not solve ODEs
 ────────────────────────────────────────────────────────────────────────────────
-✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
+✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.9s]
 ✔ |        11 | Processing of residue series
-✔ |        10 | Fitting the SFORB model [3.8s]
+✔ |        10 | Fitting the SFORB model [4.6s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         5 | Summary [0.2s]
-✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
-✔ |         9 | Hypothesis tests [8.1s]
+✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.9s]
+✔ |         9 | Hypothesis tests [11.0s]
 ✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 266.0 s
+Duration: 342.6 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
 • Fitting this ODE model with saemix takes about 15 minutes on my system (1)
 • Fitting with saemix takes around 10 minutes when using deSolve (1)
 • This still takes almost 2.5 minutes although we do not solve ODEs (1)
 
-[ FAIL 0 | WARN 0 | SKIP 3 | PASS 268 ]
+[ FAIL 0 | WARN 0 | SKIP 3 | PASS 270 ]
 Error while shutting down parallel: unable to terminate some child processes
diff --git a/man/mhmkin.Rd b/man/mhmkin.Rd
index 4230e44f..c77f4289 100644
--- a/man/mhmkin.Rd
+++ b/man/mhmkin.Rd
@@ -43,11 +43,12 @@ supported}
 
 \item{no_random_effect}{Default is NULL and will be passed to \link{saem}. If a
 character vector is supplied, it will be passed to all calls to \link{saem},
-regardless if the corresponding parameter is in the model. Alternatively,
-an object of class \link{illparms.mhmkin} can be specified. This has to have
-the same dimensions as the return object of the current call. In this way,
-ill-defined parameters found in corresponding simpler versions of the
-degradation model can be specified.}
+which will exclude random effects for all matching parameters. Alternatively,
+a list of character vectors or an object of class \link{illparms.mhmkin} can be
+specified. They have to have the same dimensions that the return object of
+the current call will have, i.e. the number of rows must match the number
+of degradation models in the mmkin object(s), and the number of columns must
+match the number of error models used in the mmkin object(s).}
 
 \item{cores}{The number of cores to be used for multicore processing. This
 is only used when the \code{cluster} argument is \code{NULL}. On Windows
@@ -74,7 +75,7 @@ be indexed using the degradation model names for the first index (row index)
 and the error model names for the second index (column index), with class
 attribute 'mhmkin'.
 
-An object of class \code{\link{mhmkin}}.
+An object inheriting from \code{\link{mhmkin}}.
 }
 \description{
 The name of the methods expresses that (\strong{m}ultiple) \strong{h}ierarchichal
@@ -82,6 +83,43 @@ The name of the methods expresses that (\strong{m}ultiple) \strong{h}ierarchicha
 fitted. Our kinetic models are nonlinear, so we can use various nonlinear
 mixed-effects model fitting functions.
 }
+\examples{
+\dontrun{
+# We start with separate evaluations of all the first six datasets with two
+# degradation models and two error models
+f_sep_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = 2, quiet = TRUE)
+f_sep_tc <- update(f_sep_const, error_model = "tc")
+# The mhmkin function sets up hierarchical degradation models aka
+# nonlinear mixed-effects models for all four combinations, specifying
+# uncorrelated random effects for all degradation parameters
+f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cores = 2)
+status(f_saem_1)
+# The 'illparms' function shows that in all hierarchical fits, at least
+# one random effect is ill-defined (the confidence interval for the
+# random effect expressed as standard deviation includes zero)
+illparms(f_saem_1)
+# Therefore we repeat the fits, excluding the ill-defined random effects
+f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1))
+status(f_saem_2)
+illparms(f_saem_2)
+# Model comparisons show that FOMC with two-component error is preferable,
+# and confirms our reduction of the default parameter model
+anova(f_saem_1)
+anova(f_saem_2)
+# The convergence plot for the selected model looks fine
+saemix::plot(f_saem_2[["FOMC", "tc"]]$so, plot.type = "convergence")
+# The plot of predictions versus data shows that we have a pretty data-rich
+# situation with homogeneous distribution of residuals, because we used the
+# same degradation model, error model and parameter distribution model that
+# was used in the data generation.
+plot(f_saem_2[["FOMC", "tc"]])
+# We can specify the same parameter model reductions manually
+no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta"))
+dim(no_ranef) <- c(2, 2)
+f_saem_2m <- update(f_saem_1, no_random_effect = no_ranef)
+anova(f_saem_2m)
+}
+}
 \seealso{
 \code{\link{[.mhmkin}} for subsetting \link{mhmkin} objects
 }
diff --git a/tests/testthat/illparms_hfits_synth.txt b/tests/testthat/illparms_hfits_synth.txt
index affd1318..7a69645b 100644
--- a/tests/testthat/illparms_hfits_synth.txt
+++ b/tests/testthat/illparms_hfits_synth.txt
@@ -1,8 +1,4 @@
            error
-degradation const                      
-       SFO                             
-       FOMC sd(log_alpha), sd(log_beta)
-           error
-degradation tc                                       
-       SFO  sd(parent_0)                             
-       FOMC sd(parent_0), sd(log_alpha), sd(log_beta)
+degradation const        tc                        
+       SFO  sd(parent_0) sd(parent_0)              
+       FOMC sd(log_beta) sd(parent_0), sd(log_beta)
diff --git a/tests/testthat/print_fits_synth_const.txt b/tests/testthat/print_fits_synth_const.txt
index b4bbe6ca..5d076d3d 100644
--- a/tests/testthat/print_fits_synth_const.txt
+++ b/tests/testthat/print_fits_synth_const.txt
@@ -4,8 +4,6 @@ Status of individual fits:
       dataset
 model  1  2  3  4  5  6 
   SFO  OK OK OK OK OK OK
-  FOMC C  OK OK OK OK C 
+  FOMC OK OK OK OK OK OK
 
-C: Optimisation did not converge:
-false convergence (8)
 OK: No warnings
diff --git a/tests/testthat/summary_hfit_sfo_tc.txt b/tests/testthat/summary_hfit_sfo_tc.txt
index 0a61f75f..0618c715 100644
--- a/tests/testthat/summary_hfit_sfo_tc.txt
+++ b/tests/testthat/summary_hfit_sfo_tc.txt
@@ -8,7 +8,7 @@ Equations:
 d_parent/dt = - k_parent * parent
 
 Data:
-104 observations of 1 variable(s) grouped in 6 datasets
+95 observations of 1 variable(s) grouped in 6 datasets
 
 Model predictions using solution type analytical 
 
@@ -19,7 +19,7 @@ Variance model: Two-component variance function
 
 Starting values for degradation parameters:
     parent_0 log_k_parent 
-         101           -3 
+          94           -2 
 
 Fixed degradation parameter values:
 None
@@ -27,7 +27,7 @@ None
 Starting values for random effects (square root of initial entries in omega):
              parent_0 log_k_parent
 parent_0            4          0.0
-log_k_parent        0          0.4
+log_k_parent        0          0.7
 
 Starting values for error model parameters:
 a.1 b.1 
@@ -37,15 +37,15 @@ Results:
 
 Likelihood computed by importance sampling
   AIC BIC logLik
-  524 523   -257
+  542 541   -266
 
 Optimised parameters:
-                  est. lower  upper
-parent_0        100.68 99.27 102.08
-log_k_parent     -3.38 -3.55  -3.21
-a.1               0.87  0.59   1.14
-b.1               0.05  0.04   0.06
-SD.log_k_parent   0.21  0.09   0.33
+                 est. lower upper
+parent_0        92.52 89.11  95.9
+log_k_parent    -1.66 -2.07  -1.3
+a.1              2.03  1.60   2.5
+b.1              0.09  0.07   0.1
+SD.log_k_parent  0.51  0.22   0.8
 
 Correlation: 
              pr_0
@@ -53,18 +53,18 @@ log_k_parent 0.1
 
 Random effects:
                 est. lower upper
-SD.log_k_parent  0.2  0.09   0.3
+SD.log_k_parent  0.5   0.2   0.8
 
 Variance model:
     est. lower upper
-a.1 0.87  0.59  1.14
-b.1 0.05  0.04  0.06
+a.1 2.03  1.60   2.5
+b.1 0.09  0.07   0.1
 
 Backtransformed parameters:
-          est. lower upper
-parent_0 1e+02 99.27 1e+02
-k_parent 3e-02  0.03 4e-02
+         est. lower upper
+parent_0 92.5  89.1  95.9
+k_parent  0.2   0.1   0.3
 
 Estimated disappearance times:
        DT50 DT90
-parent   20   68
+parent    4   12
diff --git a/tests/testthat/test_mhmkin.R b/tests/testthat/test_mhmkin.R
index e2339f28..da063326 100644
--- a/tests/testthat/test_mhmkin.R
+++ b/tests/testthat/test_mhmkin.R
@@ -3,8 +3,11 @@ context("Batch fitting and diagnosing hierarchical kinetic models")
 test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", {
 
   skip_on_cran()
-  fits_synth_const <- suppressWarnings(
-    mmkin(c("SFO", "FOMC"), ds_sfo[1:6], cores = n_cores, quiet = TRUE))
+  fits_synth_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = n_cores, quiet = TRUE)
+
+  expect_known_output(
+    print(fits_synth_const),
+    "print_fits_synth_const.txt")
 
   fits_synth_tc <- suppressWarnings(
     update(fits_synth_const, error_model = "tc"))
@@ -19,8 +22,8 @@ test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", {
     print(illparms(hfits)),
     "illparms_hfits_synth.txt")
 
-  expect_equal(which.min(AIC(hfits)), 3)
-  expect_equal(which.min(BIC(hfits)), 3)
+  expect_equal(which.min(AIC(hfits)), 4)
+  expect_equal(which.min(BIC(hfits)), 4)
 
   hfit_sfo_tc <- update(hfits[["SFO", "tc"]],
     covariance.model = diag(c(0, 1)))
@@ -38,12 +41,19 @@ test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", {
   expect_known_output(print(test_summary, digits = 1),
     "summary_hfit_sfo_tc.txt")
 
-  # It depends on the platform exactly which of the datasets fail to converge
-  # with FOMC, because they were generated to be SFO
-  skip_on_travis()
-
-  expect_known_output(
-    print(fits_synth_const),
-    "print_fits_synth_const.txt")
-
+  hfits_sfo_reduced <- update(hfits,
+    no_random_effect = illparms(hfits))
+  expect_equal(
+    as.character(illparms(hfits_sfo_reduced)),
+    rep("", 4))
+
+  # We can also manually set up an object specifying random effects to be
+  # excluded. Entries in the inital list have to be by column
+  no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta"))
+  dim(no_ranef) <- c(2, 2)
+
+  hfits_sfo_reduced_2 <- update(hfits,
+    no_random_effect = no_ranef)
+  expect_equivalent(round(anova(hfits_sfo_reduced), 0),
+    round(anova(hfits_sfo_reduced_2), 0))
 })