Simplify the removal of random effects

This is achieved by introducing the argument 'no_random_effect' to the
saem function.

5 files changed, 94 insertions, 29 deletions

diff --git a/R/saem.R b/R/saem.R
index 8929cf6f..05cce682 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -33,6 +33,12 @@ utils::globalVariables(c("predicted", "std"))
 #'   for parameter that are tested if requested by 'test_log_parms'.
 #' @param solution_type Possibility to specify the solution type in case the
 #'   automatic choice is not desired
+#' @param no_random_effect Character vector of degradation parameters for
+#'   which there should be no variability over the groups. Only used
+#'   if the covariance model is not explicitly specified.
+#' @param covariance.model Will be passed to [saemix::SaemixModel()]. Per
+#'   default, uncorrelated random effects are specified for all degradation
+#'   parameters.
 #' @param fail_with_errors Should a failure to compute standard errors
 #'   from the inverse of the Fisher Information Matrix be a failure?
 #' @param quiet Should we suppress the messages saemix prints at the beginning
@@ -121,6 +127,8 @@ saem.mmkin <- function(object,
   test_log_parms = TRUE,
   conf.level = 0.6,
   solution_type = "auto",
+  covariance.model = "auto",
+  no_random_effect = NULL,
   nbiter.saemix = c(300, 100),
   control = list(displayProgress = FALSE, print = FALSE,
     nbiter.saemix = nbiter.saemix,
@@ -134,7 +142,10 @@ saem.mmkin <- function(object,
     degparms_start = degparms_start,
     test_log_parms = test_log_parms, conf.level = conf.level,
     solution_type = solution_type,
-    transformations = transformations, ...)
+    transformations = transformations,
+    covariance.model = covariance.model,
+    no_random_effect = no_random_effect,
+    ...)
   d_saemix <- saemix_data(object, verbose = verbose)
 
   fit_failed <- FALSE
@@ -253,14 +264,22 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
 #' @rdname saem
 #' @return An [saemix::SaemixModel] object.
 #' @export
-saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"),
-  degparms_start = numeric(), test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...)
+saemix_model <- function(object, solution_type = "auto",
+  transformations = c("mkin", "saemix"), degparms_start = numeric(),
+  covariance.model = "auto", no_random_effect = NULL,
+  test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...)
 {
   if (nrow(object) > 1) stop("Only row objects allowed")
 
   mkin_model <- object[[1]]$mkinmod
 
   degparms_optim <-  mean_degparms(object, test_log_parms = test_log_parms)
+  na_degparms <- names(which(is.na(degparms_optim)))
+  if (length(na_degparms) > 0) {
+    message("Did not find valid starting values for ", paste(na_degparms, collapse = ", "), "\n",
+      "Now trying with test_log_parms = FALSE")
+    degparms_optim <-  mean_degparms(object, test_log_parms = FALSE)
+  }
   if (transformations == "saemix") {
     degparms_optim <- backtransform_odeparms(degparms_optim,
       object[[1]]$mkinmod,
@@ -378,7 +397,7 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
             }
           } else {
             model_function <- function(psi, id, xidep) {
-              tb <- exp(psi[id, 4])
+              tb <- psi[id, 4]
               t <- xidep[, "time"]
               psi[id, 1] * ifelse(t <= tb,
                 exp(- psi[id, 2] * t),
@@ -560,12 +579,22 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
   psi0_matrix <- matrix(degparms_psi0, nrow = 1)
   colnames(psi0_matrix) <- names(degparms_psi0)
 
+  if (covariance.model[1] == "auto") {
+    covariance_diagonal <- rep(1, length(degparms_optim))
+    if (!is.null(no_random_effect)) {
+      degparms_no_random <- which(names(degparms_psi0) %in% no_random_effect)
+      covariance_diagonal[degparms_no_random] <- 0
+    }
+    covariance.model = diag(covariance_diagonal)
+  }
+
   res <- saemix::saemixModel(model_function,
     psi0 = psi0_matrix,
     "Mixed model generated from mmkin object",
     transform.par = transform.par,
     error.model = error.model,
     verbose = verbose,
+    covariance.model = covariance.model,
     ...
   )
   attr(res, "mean_dp_start") <- degparms_optim
diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml
index ee93ec22..22ed5aa7 100644
--- a/docs/dev/pkgdown.yml
+++ b/docs/dev/pkgdown.yml
@@ -12,7 +12,7 @@ articles:
   compiled_models: web_only/compiled_models.html
   dimethenamid_2018: web_only/dimethenamid_2018.html
   multistart: web_only/multistart.html
-last_built: 2022-09-28T14:45Z
+last_built: 2022-09-30T10:35Z
 urls:
   reference: https://pkgdown.jrwb.de/mkin/reference
   article: https://pkgdown.jrwb.de/mkin/articles
diff --git a/log/test.log b/log/test.log
index df663253..942ed50d 100644
--- a/log/test.log
+++ b/log/test.log
@@ -5,7 +5,7 @@
 ✔ |         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 [31.6s]
 ────────────────────────────────────────────────────────────────────────────────
 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
@@ -16,7 +16,7 @@ Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
 ✔ |        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]
-✔ |         7 | Batch fitting and diagnosing hierarchical kinetic models [14.3s]
+✔ |         7 | Batch fitting and diagnosing hierarchical kinetic models [14.5s]
 ✔ |     1  12 | Nonlinear mixed-effects models [0.3s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_mixed.R:68:3): saemix results are reproducible for biphasic fits
@@ -26,26 +26,26 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |        10 | Special cases of mkinfit calls [0.4s]
 ✔ |         3 | mkinfit features [0.7s]
 ✔ |         8 | mkinmod model generation and printing [0.2s]
-✔ |         3 | Model predictions with mkinpredict [0.4s]
-✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.8s]
+✔ |         3 | Model predictions with mkinpredict [0.3s]
+✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.7s]
 ✔ |         9 | Nonlinear mixed-effects models with nlme [8.6s]
-✔ |        16 | Plotting [9.9s]
+✔ |        16 | Plotting [9.7s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |        28 | saemix parent models [181.0s]
+✔ |        35 | saemix parent models [191.8s]
 ✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
 ✔ |        11 | Processing of residue series
 ✔ |         7 | Fitting the SFORB model [3.8s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         5 | Summary [0.2s]
-✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1s]
-✔ |         9 | Hypothesis tests [8.0s]
+✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
+✔ |         9 | Hypothesis tests [7.8s]
 ✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 279.0 s
+Duration: 288.9 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)
 
-[ FAIL 0 | WARN 0 | SKIP 2 | PASS 246 ]
+[ FAIL 0 | WARN 0 | SKIP 2 | PASS 253 ]
diff --git a/man/saem.Rd b/man/saem.Rd
index dd787475..a9c4a1bb 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -18,6 +18,8 @@ saem(object, ...)
   test_log_parms = TRUE,
   conf.level = 0.6,
   solution_type = "auto",
+  covariance.model = "auto",
+  no_random_effect = NULL,
   nbiter.saemix = c(300, 100),
   control = list(displayProgress = FALSE, print = FALSE, nbiter.saemix = nbiter.saemix,
     save = FALSE, save.graphs = FALSE),
@@ -34,6 +36,8 @@ saemix_model(
   solution_type = "auto",
   transformations = c("mkin", "saemix"),
   degparms_start = numeric(),
+  covariance.model = "auto",
+  no_random_effect = NULL,
   test_log_parms = FALSE,
   conf.level = 0.6,
   verbose = FALSE,
@@ -70,6 +74,14 @@ for parameter that are tested if requested by 'test_log_parms'.}
 \item{solution_type}{Possibility to specify the solution type in case the
 automatic choice is not desired}
 
+\item{covariance.model}{Will be passed to \code{\link[saemix:SaemixModel-class]{saemix::SaemixModel()}}. Per
+default, uncorrelated random effects are specified for all degradation
+parameters.}
+
+\item{no_random_effect}{Character vector of degradation parameters for
+which there should be no variability over the groups. Only used
+if the covariance model is not explicitly specified.}
+
 \item{nbiter.saemix}{Convenience option to increase the number of
 iterations}
 
diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R
index 5d9a01de..ce776bf7 100644
--- a/tests/testthat/test_saemix_parent.R
+++ b/tests/testthat/test_saemix_parent.R
@@ -4,30 +4,54 @@ test_that("Parent fits using saemix are correctly implemented", {
 
   skip_on_cran()
   expect_error(saem(fits), "Only row objects")
-  # Some fits were done in the setup script
+
+  # mmkin_sfo_1 was generated in the setup script
+  # We did not introduce variance of parent_0 in the data generation
+  # This is correctly detected
+  expect_equal(illparms(sfo_saem_1), "sd(parent_0)")
+  # So we remove this variance
+  sfo_saem_1_reduced <- update(sfo_saem_1, no_random_effect = "parent_0")
+  expect_equal(illparms(sfo_saem_1_reduced), character(0))
+
+  # We cannot currently do the fit with completely fixed initial values
   mmkin_sfo_2 <- update(mmkin_sfo_1, fixed_initials = c(parent = 100))
+  sfo_saem_3 <- expect_error(saem(mmkin_sfo_2, quiet = TRUE), "at least two parameters")
+
+  # We get an error if we do not supply a suitable model specification
   expect_error(update(mmkin_sfo_1, models = c("SFOOO")), "Please supply models.*")
 
-  sfo_saem_2 <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "mkin")
-  sfo_saem_3 <- expect_error(saem(mmkin_sfo_2, quiet = TRUE), "at least two parameters")
-  expect_equal(endpoints(sfo_saem_1), endpoints(sfo_saem_2), tol = 0.01)
-  s_sfo_s1 <- summary(sfo_saem_1)
-  s_sfo_s2 <- summary(sfo_saem_2)
+  sfo_saem_1_mkin <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "mkin")
+  expect_equal(illparms(sfo_saem_1_mkin), "sd(parent_0)")
+  sfo_saem_1_reduced_mkin <- update(sfo_saem_1_mkin, no_random_effect = "parent_0")
+
+  # The endpoints obtained do not depend on the transformation
+  expect_equal(endpoints(sfo_saem_1), endpoints(sfo_saem_1_mkin), tol = 0.01)
+  expect_equal(endpoints(sfo_saem_1_reduced), endpoints(sfo_saem_1_reduced_mkin), tol = 0.01)
+
+  s_sfo_saem_1 <- summary(sfo_saem_1)
+  s_sfo_saem_1_reduced <- summary(sfo_saem_1_reduced)
+  s_sfo_saem_1_mkin <- summary(sfo_saem_1_mkin)
+  s_sfo_saem_1_reduced_mkin <- summary(sfo_saem_1_reduced_mkin)
 
   sfo_nlme_1 <- expect_warning(nlme(mmkin_sfo_1), "not converge")
-  s_sfo_n <- summary(sfo_nlme_1)
+  s_sfo_nlme_1 <- summary(sfo_nlme_1)
 
   # Compare with input
-  expect_equal(round(s_sfo_s2$confint_ranef["SD.log_k_parent", "est."], 1), 0.3)
+  expect_equal(round(s_sfo_saem_1$confint_ranef["SD.k_parent", "est."], 1), 0.3)
+  expect_equal(round(s_sfo_saem_1_mkin$confint_ranef["SD.log_k_parent", "est."], 1), 0.3)
   # k_parent is a bit different from input 0.03 here
-  expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), 0.035)
-  expect_equal(round(s_sfo_s2$confint_back["k_parent", "est."], 3), 0.035)
+  expect_equal(round(s_sfo_saem_1$confint_back["k_parent", "est."], 3), 0.035)
+  expect_equal(round(s_sfo_saem_1_mkin$confint_back["k_parent", "est."], 3), 0.035)
 
   # But the result is pretty unanimous between methods
-  expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3),
-    round(s_sfo_s2$confint_back["k_parent", "est."], 3))
-  expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3),
-    round(s_sfo_n$confint_back["k_parent", "est."], 3))
+  expect_equal(round(s_sfo_saem_1_reduced$confint_back["k_parent", "est."], 3),
+    round(s_sfo_saem_1$confint_back["k_parent", "est."], 3))
+  expect_equal(round(s_sfo_saem_1_mkin$confint_back["k_parent", "est."], 3),
+    round(s_sfo_saem_1$confint_back["k_parent", "est."], 3))
+  expect_equal(round(s_sfo_saem_1_reduced_mkin$confint_back["k_parent", "est."], 3),
+    round(s_sfo_saem_1$confint_back["k_parent", "est."], 3))
+  expect_equal(round(s_sfo_nlme_1$confint_back["k_parent", "est."], 3),
+    round(s_sfo_saem_1$confint_back["k_parent", "est."], 3))
 
   mmkin_fomc_1 <- mmkin("FOMC", ds_fomc, quiet = TRUE, error_model = "tc", cores = n_cores)
   fomc_saem_1 <- saem(mmkin_fomc_1, quiet = TRUE, transformations = "saemix")