From ae8ba4b0e52aae9b317b0244e7162037bee9d27b Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 2 Dec 2020 07:53:37 +0100 Subject: Possibility to specify random effects structures The default is pdDiag again, as we often have a small number of datasets in degradation kinetics. --- R/nlme.mmkin.R | 30 +++++++++++++----------------- man/nlme.mmkin.Rd | 9 +++++---- test.log | 22 +++++++++------------- tests/testthat/test_nlme.R | 12 ++++++------ 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)) -- cgit v1.2.1