Possibility to specify random effects structures

The default is pdDiag again, as we often have a small number of datasets
in degradation kinetics.

4 files changed, 33 insertions, 40 deletions

diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index 0d6e6c37..c6f15c8e 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -130,9 +130,13 @@ get_deg_func <- function() {
 #'   anova(f_nlme_dfop_sfo, f_nlme_dfop_sfo_obs, f_nlme_dfop_sfo_tc)
 #'
 #' }
-nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
-  fixed, random = fixed,
-  groups, start, correlation = NULL, weights = NULL,
+nlme.mmkin <- function(model, data = "auto",
+  fixed = lapply(as.list(names(mean_degparms(model))),
+    function(el) eval(parse(text = paste(el, 1, sep = "~")))),
+  random = pdDiag(fixed),
+  groups,
+  start = mean_degparms(model, random = TRUE),
+  correlation = NULL, weights = NULL,
   subset, method = c("ML", "REML"),
   na.action = na.fail, naPattern,
   control = list(), verbose= FALSE)
@@ -143,8 +147,8 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
 
   # Warn in case arguments were used that are overriden
   if (any(!is.na(match(names(thisCall),
-               c("fixed", "data"))))) {
-    warning("'nlme.mmkin' will redefine 'fixed' and 'data'")
+               c("data"))))) {
+    warning("'nlme.mmkin' will redefine 'data'")
   }
 
   deg_func <- nlme_function(model)
@@ -158,21 +162,13 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
 
   thisCall[["model"]] <- this_model
 
-  mean_dp_start <- mean_degparms(model)
-  dp_names <- names(mean_dp_start)
-
   thisCall[["data"]] <- nlme_data(model)
 
-  if (missing(start)) {
-    thisCall[["start"]] <- mean_degparms(model, random = TRUE)
-  }
+  thisCall[["start"]] <- start
 
-  thisCall[["fixed"]] <- lapply(as.list(dp_names), function(el)
-                                 eval(parse(text = paste(el, 1, sep = "~"))))
+  thisCall[["fixed"]] <- fixed
 
-  if (missing(random)) {
-    thisCall[["random"]] <- pdLogChol(thisCall[["fixed"]])
-  }
+  thisCall[["random"]] <- random
 
   error_model <- model[[1]]$err_mod
 
@@ -198,7 +194,7 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
   val$mkinmod <- model[[1]]$mkinmod
   val$data <- thisCall[["data"]]
   val$mmkin <- model
-  val$mean_dp_start <- mean_dp_start
+  val$mean_dp_start <- start$fixed
   val$transform_rates <- model[[1]]$transform_rates
   val$transform_fractions <- model[[1]]$transform_fractions
   val$solution_type <- model[[1]]$solution_type
diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd
index 7a72a19a..f78256ac 100644
--- a/man/nlme.mmkin.Rd
+++ b/man/nlme.mmkin.Rd
@@ -8,11 +8,12 @@
 \usage{
 \method{nlme}{mmkin}(
   model,
-  data = sys.frame(sys.parent()),
-  fixed,
-  random = fixed,
+  data = "auto",
+  fixed = lapply(as.list(names(mean_degparms(model))), function(el) eval(parse(text =
+    paste(el, 1, sep = "~")))),
+  random = pdDiag(fixed),
   groups,
-  start,
+  start = mean_degparms(model, random = TRUE),
   correlation = NULL,
   weights = NULL,
   subset,
diff --git a/test.log b/test.log
index 84a8748e..c977d6a0 100644
--- a/test.log
+++ b/test.log
@@ -6,17 +6,17 @@ Testing mkin
 ✔ |   2       | Export dataset for reading into CAKE
 ✔ |  14       | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.9 s]
 ✔ |   4       | Calculation of FOCUS chi2 error levels [0.4 s]
-✔ |   7       | Fitting the SFORB model [3.4 s]
+✔ |   7       | Fitting the SFORB model [3.5 s]
 ✔ |   5       | Analytical solutions for coupled models [3.1 s]
 ✔ |   5       | Calculation of Akaike weights
 ✔ |  10       | Confidence intervals and p-values [1.0 s]
-✔ |  14       | Error model fitting [4.4 s]
+✔ |  14       | Error model fitting [4.7 s]
 ✔ |   4       | Test fitting the decline of metabolites from their maximum [0.3 s]
 ✔ |   1       | Fitting the logistic model [0.2 s]
 ✔ |   1       | Test dataset class mkinds used in gmkin
 ✔ |   1       | mkinfit features [0.3 s]
 ✔ |  12       | Special cases of mkinfit calls [0.7 s]
-✔ |   8       | mkinmod model generation and printing [0.2 s]
+✔ |   8       | mkinmod model generation and printing [0.3 s]
 ✔ |   3       | Model predictions with mkinpredict [0.4 s]
 ✔ |  14     2 | Evaluations according to 2015 NAFTA guidance [1.2 s]
 ────────────────────────────────────────────────────────────────────────────────
@@ -26,28 +26,24 @@ Reason: getRversion() < "4.1.0"  is TRUE
 Skip (test_nafta.R:55:5): Test data from Appendix D are correctly evaluated
 Reason: getRversion() < "4.1.0"  is TRUE
 ────────────────────────────────────────────────────────────────────────────────
-✖ |   7 1     | Nonlinear mixed-effects models [8.0 s]
-────────────────────────────────────────────────────────────────────────────────
-FAILURE (test_nlme.R:90:3): nlme_function works correctly
-`tmp <- update(m_nlme_mmkin)` did not produce any warnings.
-────────────────────────────────────────────────────────────────────────────────
+✔ |   9       | Nonlinear mixed-effects models [8.0 s]
 ✔ |   0     1 | Plotting [0.7 s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_plot.R:25:3): Plotting mkinfit and mmkin objects is reproducible
 Reason: getRversion() < "4.1.0"  is TRUE
 ────────────────────────────────────────────────────────────────────────────────
 ✔ |   4       | Residuals extracted from mkinfit models
-✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s]
-✔ |   4       | Summary [0.1 s]
+✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [1.7 s]
+✔ |   4       | Summary [0.2 s]
 ✔ |   1       | Summaries of old mkinfit objects
-✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2 s]
+✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3 s]
 ✔ |   9       | Hypothesis tests [7.2 s]
 ✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 38.9 s
+Duration: 39.8 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
 ● getRversion() < "4.1.0"  is TRUE (3)
 
-[ FAIL 1 | WARN 0 | SKIP 3 | PASS 144 ]
+[ FAIL 0 | WARN 0 | SKIP 3 | PASS 146 ]
diff --git a/tests/testthat/test_nlme.R b/tests/testthat/test_nlme.R
index da994f49..989914da 100644
--- a/tests/testthat/test_nlme.R
+++ b/tests/testthat/test_nlme.R
@@ -56,18 +56,18 @@ test_that("nlme_function works correctly", {
   m_nlme_raw_2 <- nlme(value ~ SSasymp(time, 0, parent_0, log_k_parent),
     data = grouped_data,
     fixed = parent_0 + log_k_parent ~ 1,
-    random = pdLogChol(parent_0 + log_k_parent ~ 1),
+    random = pdDiag(parent_0 + log_k_parent ~ 1),
     start = mean_degparms(f, random = TRUE),
     control = list("msWarnNoConv" = FALSE))
 
   expect_equal(m_nlme_raw_2$coefficients, m_nlme_mmkin$coefficients)
 
-  anova_nlme <- anova(m_nlme_raw, m_nlme_mmkin)
+  anova_nlme <- anova(m_nlme_raw, m_nlme_mkin, m_nlme_raw_2, m_nlme_mmkin)
 
-  # We get a slightly lower AIC with the improved starting values used within
-  # nlme.mmkin, specifying also random effects
-  expect_lt(anova_nlme["m_nlme_mmkin", "AIC"],
+  expect_equal(anova_nlme["m_nlme_mkin", "AIC"],
     anova_nlme["m_nlme_raw", "AIC"])
+  expect_equal(anova_nlme["m_nlme_mmkin", "AIC"],
+    anova_nlme["m_nlme_raw_2", "AIC"])
 
   m_nlme_raw_up_1 <- update(m_nlme_raw, random = log_k_parent_sink ~ 1)
   # The following three calls give an error although they should
@@ -87,7 +87,7 @@ test_that("nlme_function works correctly", {
   m_nlme_mkin_up_2 <- update(m_nlme_mkin, random = parent_0 ~ 1)
   expect_equal(m_nlme_raw_up_2$coefficients, m_nlme_mkin_up_2$coefficients)
 
-  expect_warning(tmp <- update(m_nlme_mmkin), "Iteration 1, LME step")
+  expect_silent(tmp <- update(m_nlme_mmkin))
 
   geomean_dt50_mmkin <- exp(mean(log((sapply(f, function(x) endpoints(x)$distimes["parent", "DT50"])))))
   expect_equal(round(endpoints(m_nlme_mmkin)$distimes["parent", "DT50"]), round(geomean_dt50_mmkin))