Attempt at automatic setting of random effects

Based on parameters in the separate fits that fail the t-test.

7 files changed, 85 insertions, 30 deletions

diff --git a/R/mhmkin.R b/R/mhmkin.R
index de41c733..1f29dc40 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -12,6 +12,13 @@
 #' degradation models to the same data
 #' @param backend The backend to be used for fitting. Currently, only saemix is
 #' supported
+#' @param no_random_effect Default is NULL and will be passed to [saem]. If
+#' you specify "auto", random effects are only included if the number
+#' of datasets in which the parameter passed the t-test is at least 'auto_ranef_threshold'.
+#' Beware that while this may make for convenient model reduction or even
+#' numerical stability of the algorithm, it will likely lead to
+#' underparameterised models.
+#' @param auto_ranef_threshold See 'no_random_effect.
 #' @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.
@@ -31,7 +38,7 @@
 #' @author Johannes Ranke
 #' @seealso \code{\link{[.mhmkin}} for subsetting [mhmkin] objects
 #' @export
-mhmkin <- function(objects, backend = "saemix", algorithm = "saem", ...) {
+mhmkin <- function(objects, ...) {
   UseMethod("mhmkin")
 }
 
@@ -43,7 +50,8 @@ mhmkin.mmkin <- function(objects, ...) {
 
 #' @export
 #' @rdname mhmkin
-mhmkin.list <- function(objects, backend = "saemix",
+mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
+  no_random_effect = NULL, auto_ranef_threshold = 3,
   ...,
   cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), cluster = NULL)
 {
@@ -94,7 +102,19 @@ mhmkin.list <- function(objects, backend = "saemix",
     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)))
+    if (identical(no_random_effect, "auto")) {
+      ip <- illparms(mmkin_row)
+      ipt <- table(unlist(lapply(as.list(ip), strsplit, ", ")))
+      no_ranef <- names(ipt)[which((length(ds) - ipt) <= auto_ranef_threshold)]
+      transparms <- rownames(mmkin_row[[1]]$start_transformed)
+      bparms <- rownames(mmkin_row[[1]]$start)
+      names(transparms) <- bparms
+      no_ranef_trans <- transparms[no_ranef]
+    } else {
+      no_ranef_trans <- no_random_effect
+    }
+    res <- try(do.call(backend_function,
+        args = c(list(mmkin_row), dot_args, list(no_random_effect = no_ranef_trans))))
     return(res)
   }
 
@@ -207,7 +227,7 @@ anova.mhmkin <- function(object, ...,
   if (identical(model.names, "auto")) {
     model.names <- outer(rownames(object), colnames(object), paste)
   }
-  rlang::inject(anova(!!!(object), method = method, test = test, 
+  rlang::inject(anova(!!!(object), method = method, test = test,
       model.names = model.names))
 }
 
diff --git a/log/check.log b/log/check.log
index 3fea2ec6..d8225158 100644
--- a/log/check.log
+++ b/log/check.log
@@ -57,12 +57,12 @@ 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 ... [15s/15s] OK
+* checking examples ... [36s/46s] OK
 * checking for unstated dependencies in ‘tests’ ... OK
 * checking tests ... SKIPPED
 * checking for unstated dependencies in vignettes ... OK
 * checking package vignettes in ‘inst/doc’ ... OK
-* checking re-building of vignette outputs ... OK
+* checking re-building of vignette outputs ... [26s/41s] OK
 * checking PDF version of manual ... OK
 * checking HTML version of manual ... OK
 * checking for non-standard things in the check directory ... OK
diff --git a/log/test.log b/log/test.log
index 611b66b6..95d50017 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,57 +1,57 @@
 ℹ Testing mkin
 ✔ | F W S  OK | Context
 ✔ |         5 | AIC calculation
-✔ |         5 | Analytical solutions for coupled models [3.3s]
+✔ |         5 | Analytical solutions for coupled models [3.4s]
 ✔ |         5 | Calculation of Akaike weights
 ✔ |         3 | Export dataset for reading into CAKE
 ✔ |        12 | Confidence intervals and p-values [1.0s]
-✔ |     1  12 | Dimethenamid data from 2018 [31.9s]
+✔ |     1  12 | Dimethenamid data from 2018 [33.0s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_dmta.R:98:3): 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 [5.6s]
 ✔ |         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]
+✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.9s]
+✔ |         4 | Test fitting the decline of metabolites from their maximum [0.4s]
 ✔ |         1 | Fitting the logistic model [0.2s]
-✔ |         8 | Batch fitting and diagnosing hierarchical kinetic models [14.5s]
+✔ |        10 | Batch fitting and diagnosing hierarchical kinetic models [26.9s]
 ✔ |     1  12 | Nonlinear mixed-effects models [0.3s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_mixed.R:74:3): 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.6s]
+✔ |         3 | mkinfit features [0.8s]
+✔ |         8 | mkinmod model generation and printing [0.3s]
 ✔ |         3 | Model predictions with mkinpredict [0.4s]
-✔ |         7 | Multistart method for saem.mmkin models [36.2s]
-✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.4s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [8.8s]
-✔ |        16 | Plotting [10.0s]
+✔ |         7 | Multistart method for saem.mmkin models [44.7s]
+✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.8s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [9.9s]
+✔ |        16 | Plotting [10.5s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |     1  37 | saemix parent models [71.4s]
+✔ |     1  37 | saemix parent models [76.1s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_saemix_parent.R:153:3): 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.5s]
 ✔ |        11 | Processing of residue series
-✔ |         7 | Fitting the SFORB model [3.7s]
+✔ |         7 | Fitting the SFORB model [3.9s]
 ✔ |         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.3s]
+✔ |         9 | Hypothesis tests [8.4s]
 ✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 206.7 s
+Duration: 237.3 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 263 ]
+[ FAIL 0 | WARN 0 | SKIP 3 | PASS 265 ]
diff --git a/man/mhmkin.Rd b/man/mhmkin.Rd
index 1a6e3869..597fa8ac 100644
--- a/man/mhmkin.Rd
+++ b/man/mhmkin.Rd
@@ -9,13 +9,16 @@
 \title{Fit nonlinear mixed-effects models built from one or more kinetic
 degradation models and one or more error models}
 \usage{
-mhmkin(objects, backend = "saemix", algorithm = "saem", ...)
+mhmkin(objects, ...)
 
 \method{mhmkin}{mmkin}(objects, ...)
 
 \method{mhmkin}{list}(
   objects,
   backend = "saemix",
+  algorithm = "saem",
+  no_random_effect = NULL,
+  auto_ranef_threshold = 3,
   ...,
   cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(),
   cluster = NULL
@@ -31,13 +34,22 @@ degradation models to the same data, but using different error models.
 Alternatively, a single \link{mmkin} object containing fits of several
 degradation models to the same data}
 
+\item{\dots}{Further arguments that will be passed to the nonlinear mixed-effects
+model fitting function.}
+
 \item{backend}{The backend to be used for fitting. Currently, only saemix is
 supported}
 
 \item{algorithm}{The algorithm to be used for fitting (currently not used)}
 
-\item{\dots}{Further arguments that will be passed to the nonlinear mixed-effects
-model fitting function.}
+\item{no_random_effect}{Default is NULL and will be passed to \link{saem}. If
+you specify "auto", random effects are only included if the number
+of datasets in which the parameter passed the t-test is at least 'auto_ranef_threshold'.
+Beware that while this may make for convenient model reduction or even
+numerical stability of the algorithm, it will likely lead to
+underparameterised models.}
+
+\item{auto_ranef_threshold}{See 'no_random_effect.}
 
 \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
diff --git a/tests/testthat/illparms_hfits_synth_no_ranef_auto.txt b/tests/testthat/illparms_hfits_synth_no_ranef_auto.txt
new file mode 100644
index 00000000..a64ed222
--- /dev/null
+++ b/tests/testthat/illparms_hfits_synth_no_ranef_auto.txt
@@ -0,0 +1,4 @@
+           error
+degradation const tc          
+       SFO        sd(parent_0)
+       FOMC       b.1         
diff --git a/tests/testthat/print_hfits_synth_no_ranef_auto.txt b/tests/testthat/print_hfits_synth_no_ranef_auto.txt
new file mode 100644
index 00000000..9af1cbcd
--- /dev/null
+++ b/tests/testthat/print_hfits_synth_no_ranef_auto.txt
@@ -0,0 +1,9 @@
+<mhmkin> object
+Status of individual fits:
+
+           error
+degradation const tc
+       SFO  OK    OK
+       FOMC OK    OK
+
+OK: Fit terminated successfully
diff --git a/tests/testthat/test_mhmkin.R b/tests/testthat/test_mhmkin.R
index 9fd29c50..93333ac1 100644
--- a/tests/testthat/test_mhmkin.R
+++ b/tests/testthat/test_mhmkin.R
@@ -38,12 +38,22 @@ 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
-  # SFO datasets fail to converge with FOMC
+  # 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_no_ranef_auto <- update(hfits, no_random_effect = "auto", auto_ranef_threshold = 2)
+
+  expect_known_output(
+    print(hfits_no_ranef_auto),
+    "print_hfits_synth_no_ranef_auto.txt")
+
+  expect_known_output(
+    print(illparms(hfits_no_ranef_auto)),
+    "illparms_hfits_synth_no_ranef_auto.txt")
+
 })