diff options
| -rw-r--r-- | R/nlme.R | 1 | ||||
| -rw-r--r-- | R/nlme.mmkin.R | 47 | ||||
| -rw-r--r-- | R/summary.nlme.mmkin.R | 1 | ||||
| -rw-r--r-- | man/mkinmod.Rd | 6 | ||||
| -rw-r--r-- | man/nlme.Rd | 1 | ||||
| -rw-r--r-- | man/nlme.mmkin.Rd | 45 | ||||
| -rw-r--r-- | man/summary.nlme.mmkin.Rd | 1 | ||||
| -rw-r--r-- | test.log | 22 | ||||
| -rw-r--r-- | tests/testthat/test_nlme.R | 20 |
9 files changed, 86 insertions, 58 deletions
@@ -54,7 +54,6 @@ #' # The procedure is greatly simplified by the nlme.mmkin function #' f_nlme <- nlme(f) #' plot(f_nlme) -#' #' @return A function that can be used with nlme #' @export nlme_function <- function(object) { diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index a9e1694f..b434bbf4 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -56,15 +56,15 @@ get_deg_func <- function() { #' f <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, cores = 1) #' library(nlme) #' f_nlme_sfo <- nlme(f["SFO", ]) -#' f_nlme_dfop <- nlme(f["DFOP", ]) -#' AIC(f_nlme_sfo, f_nlme_dfop) -#' print(f_nlme_dfop) -#' plot(f_nlme_dfop) -#' endpoints(f_nlme_dfop) #' #' \dontrun{ -#' f_nlme_2 <- nlme(f["SFO", ], start = c(parent_0 = 100, log_k_parent = 0.1)) -#' update(f_nlme_2, random = parent_0 ~ 1) +#' +#' f_nlme_dfop <- nlme(f["DFOP", ]) +#' anova(f_nlme_sfo, f_nlme_dfop) +#' print(f_nlme_dfop) +#' plot(f_nlme_dfop) +#' endpoints(f_nlme_dfop) +#' #' ds_2 <- lapply(experimental_data_for_UBA_2019[6:10], #' function(x) x$data[c("name", "time", "value")]) #' m_sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), @@ -82,14 +82,15 @@ get_deg_func <- function() { #' f_nlme_sfo_sfo <- nlme(f_2["SFO-SFO", ]) #' plot(f_nlme_sfo_sfo) #' -#' # With formation fractions -#' f_nlme_sfo_sfo_ff <- nlme(f_2["SFO-SFO-ff", ]) -#' plot(f_nlme_sfo_sfo_ff) +#' # With formation fractions this does not coverge with defaults +#' # f_nlme_sfo_sfo_ff <- nlme(f_2["SFO-SFO-ff", ]) +#' #plot(f_nlme_sfo_sfo_ff) #' -#' # For the following fit we need to increase pnlsMaxIter and the tolerance -#' # to get convergence -#' f_nlme_dfop_sfo <- nlme(f_2["DFOP-SFO", ], -#' control = list(pnlsMaxIter = 120, tolerance = 5e-4), verbose = TRUE) +#' # With the log-Cholesky parameterization, this converges in 11 +#' # iterations and around 100 seconds, but without tweaking control +#' # parameters (with pdDiag, increasing the tolerance and pnlsMaxIter was +#' # necessary) +#' f_nlme_dfop_sfo <- nlme(f_2["DFOP-SFO", ]) #' #' plot(f_nlme_dfop_sfo) #' @@ -114,10 +115,18 @@ get_deg_func <- function() { #' ds_2, quiet = TRUE, error_model = "obs") #' f_nlme_sfo_sfo_obs <- nlme(f_2_obs["SFO-SFO", ]) #' print(f_nlme_sfo_sfo_obs) -#' # The same with DFOP-SFO does not converge, apparently the variances of -#' # parent and A1 are too similar in this case, so that the model is -#' # overparameterised -#' #f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ], control = list(maxIter = 100)) +#' f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ]) +#' +#' f_2_tc <- mmkin(list("SFO-SFO" = m_sfo_sfo, +#' "DFOP-SFO" = m_dfop_sfo), +#' ds_2, quiet = TRUE, error_model = "tc") +#' # f_nlme_sfo_sfo_tc <- nlme(f_2_tc["SFO-SFO", ]) # stops with error message +#' f_nlme_dfop_sfo_tc <- nlme(f_2_tc["DFOP-SFO", ]) +#' # We get warnings about false convergence in the LME step in several iterations +#' # but as the last such warning occurs in iteration 25 and we have 28 iterations +#' # we can ignore these +#' 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, @@ -160,7 +169,7 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()), eval(parse(text = paste(el, 1, sep = "~")))) if (missing(random)) { - thisCall[["random"]] <- pdDiag(thisCall[["fixed"]]) + thisCall[["random"]] <- pdLogChol(thisCall[["fixed"]]) } error_model <- model[[1]]$err_mod diff --git a/R/summary.nlme.mmkin.R b/R/summary.nlme.mmkin.R index 42326b39..29f1207b 100644 --- a/R/summary.nlme.mmkin.R +++ b/R/summary.nlme.mmkin.R @@ -55,6 +55,7 @@ #' ds_sfo_mean <- lapply(k_in, pred_sfo) #' names(ds_sfo_mean) <- paste("ds", 1:5) #' +#' set.seed(12345) #' ds_sfo_syn <- lapply(ds_sfo_mean, function(ds) { #' add_err(ds, #' sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2), diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd index 77319aac..95cec09a 100644 --- a/man/mkinmod.Rd +++ b/man/mkinmod.Rd @@ -155,13 +155,17 @@ print(SFO_SFO) fit_sfo_sfo <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE, solution_type = "deSolve") # Now supplying compound names used for plotting, and write to user defined location + # We need to choose a path outside the session tempdir because this gets removed + DLL_dir <- "~/.local/share/mkin" + if (!dir.exists(DLL_dir)) dir.create(DLL_dir) SFO_SFO.2 <- mkinmod( parent = mkinsub("SFO", "m1", full_name = "Test compound"), m1 = mkinsub("SFO", full_name = "Metabolite M1"), - name = "SFO_SFO", dll_dir = "~/dll", unload = TRUE, overwrite = TRUE) + name = "SFO_SFO", dll_dir = DLL_dir, unload = TRUE, overwrite = TRUE) # Now we can save the model and restore it in a new session saveRDS(SFO_SFO.2, file = "~/SFO_SFO.rds") # Terminate the R session here if you would like to check, and then do +library(mkin) SFO_SFO.3 <- readRDS("~/SFO_SFO.rds") fit_sfo_sfo <- mkinfit(SFO_SFO.3, FOCUS_2006_D, quiet = TRUE, solution_type = "deSolve") diff --git a/man/nlme.Rd b/man/nlme.Rd index df721a0f..307cca82 100644 --- a/man/nlme.Rd +++ b/man/nlme.Rd @@ -78,7 +78,6 @@ plot(augPred(m_nlme, level = 0:1), layout = c(3, 1)) # The procedure is greatly simplified by the nlme.mmkin function f_nlme <- nlme(f) plot(f_nlme) - } \seealso{ \code{\link{nlme.mmkin}} diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd index abcd0e81..0a9f6913 100644 --- a/man/nlme.mmkin.Rd +++ b/man/nlme.mmkin.Rd @@ -87,15 +87,15 @@ ds <- lapply(experimental_data_for_UBA_2019[6:10], f <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, cores = 1) library(nlme) f_nlme_sfo <- nlme(f["SFO", ]) -f_nlme_dfop <- nlme(f["DFOP", ]) -AIC(f_nlme_sfo, f_nlme_dfop) -print(f_nlme_dfop) -plot(f_nlme_dfop) -endpoints(f_nlme_dfop) \dontrun{ - f_nlme_2 <- nlme(f["SFO", ], start = c(parent_0 = 100, log_k_parent = 0.1)) - update(f_nlme_2, random = parent_0 ~ 1) + + f_nlme_dfop <- nlme(f["DFOP", ]) + anova(f_nlme_sfo, f_nlme_dfop) + print(f_nlme_dfop) + plot(f_nlme_dfop) + endpoints(f_nlme_dfop) + ds_2 <- lapply(experimental_data_for_UBA_2019[6:10], function(x) x$data[c("name", "time", "value")]) m_sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), @@ -113,14 +113,15 @@ endpoints(f_nlme_dfop) f_nlme_sfo_sfo <- nlme(f_2["SFO-SFO", ]) plot(f_nlme_sfo_sfo) - # With formation fractions - f_nlme_sfo_sfo_ff <- nlme(f_2["SFO-SFO-ff", ]) - plot(f_nlme_sfo_sfo_ff) + # With formation fractions this does not coverge with defaults + # f_nlme_sfo_sfo_ff <- nlme(f_2["SFO-SFO-ff", ]) + #plot(f_nlme_sfo_sfo_ff) - # For the following fit we need to increase pnlsMaxIter and the tolerance - # to get convergence - f_nlme_dfop_sfo <- nlme(f_2["DFOP-SFO", ], - control = list(pnlsMaxIter = 120, tolerance = 5e-4), verbose = TRUE) + # With the log-Cholesky parameterization, this converges in 11 + # iterations and around 100 seconds, but without tweaking control + # parameters (with pdDiag, increasing the tolerance and pnlsMaxIter was + # necessary) + f_nlme_dfop_sfo <- nlme(f_2["DFOP-SFO", ]) plot(f_nlme_dfop_sfo) @@ -145,10 +146,18 @@ endpoints(f_nlme_dfop) ds_2, quiet = TRUE, error_model = "obs") f_nlme_sfo_sfo_obs <- nlme(f_2_obs["SFO-SFO", ]) print(f_nlme_sfo_sfo_obs) - # The same with DFOP-SFO does not converge, apparently the variances of - # parent and A1 are too similar in this case, so that the model is - # overparameterised - #f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ], control = list(maxIter = 100)) + f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ]) + + f_2_tc <- mmkin(list("SFO-SFO" = m_sfo_sfo, + "DFOP-SFO" = m_dfop_sfo), + ds_2, quiet = TRUE, error_model = "tc") + # f_nlme_sfo_sfo_tc <- nlme(f_2_tc["SFO-SFO", ]) # stops with error message + f_nlme_dfop_sfo_tc <- nlme(f_2_tc["DFOP-SFO", ]) + # We get warnings about false convergence in the LME step in several iterations + # but as the last such warning occurs in iteration 25 and we have 28 iterations + # we can ignore these + anova(f_nlme_dfop_sfo, f_nlme_dfop_sfo_obs, f_nlme_dfop_sfo_tc) + } } \seealso{ diff --git a/man/summary.nlme.mmkin.Rd b/man/summary.nlme.mmkin.Rd index ea625dd7..d7e61074 100644 --- a/man/summary.nlme.mmkin.Rd +++ b/man/summary.nlme.mmkin.Rd @@ -80,6 +80,7 @@ pred_sfo <- function(k) { ds_sfo_mean <- lapply(k_in, pred_sfo) names(ds_sfo_mean) <- paste("ds", 1:5) +set.seed(12345) ds_sfo_syn <- lapply(ds_sfo_mean, function(ds) { add_err(ds, sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2), @@ -6,11 +6,11 @@ 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.2 s] -✔ | 5 | Analytical solutions for coupled models [2.9 s] +✔ | 7 | Fitting the SFORB model [3.4 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.3 s] +✔ | 14 | Error model fitting [4.4 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 @@ -26,24 +26,28 @@ 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 ──────────────────────────────────────────────────────────────────────────────── -✔ | 9 | Nonlinear mixed-effects models [7.8 s] +✖ | 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. +──────────────────────────────────────────────────────────────────────────────── ✔ | 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.4 s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s] ✔ | 4 | Summary [0.1 s] ✔ | 1 | Summaries of old mkinfit objects -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.0 s] -✔ | 9 | Hypothesis tests [6.4 s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2 s] +✔ | 9 | Hypothesis tests [7.2 s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 36.9 s +Duration: 38.9 s ── Skipped tests ────────────────────────────────────────────────────────────── ● getRversion() < "4.1.0" is TRUE (3) -[ FAIL 0 | WARN 0 | SKIP 3 | PASS 146 ] +[ FAIL 1 | WARN 0 | SKIP 3 | PASS 144 ] diff --git a/tests/testthat/test_nlme.R b/tests/testthat/test_nlme.R index 4fd8e53f..da994f49 100644 --- a/tests/testthat/test_nlme.R +++ b/tests/testthat/test_nlme.R @@ -38,24 +38,27 @@ test_that("nlme_function works correctly", { m_nlme_raw <- nlme(value ~ SSasymp(time, 0, parent_0, log_k_parent_sink), data = grouped_data, fixed = parent_0 + log_k_parent_sink ~ 1, - random = pdDiag(parent_0 + log_k_parent_sink ~ 1), - start = mean_dp) + random = pdLogChol(parent_0 + log_k_parent_sink ~ 1), + start = mean_dp, + control = list("msWarnNoConv" = FALSE)) m_nlme_mkin <- nlme(value ~ nlme_f(name, time, parent_0, log_k_parent_sink), data = grouped_data, fixed = parent_0 + log_k_parent_sink ~ 1, - random = pdDiag(parent_0 + log_k_parent_sink ~ 1), - start = mean_dp) + random = pdLogChol(parent_0 + log_k_parent_sink ~ 1), + start = mean_dp, + control = list("msWarnNoConv" = FALSE)) expect_equal(m_nlme_raw$coefficients, m_nlme_mkin$coefficients) - m_nlme_mmkin <- nlme(f) + m_nlme_mmkin <- nlme(f, control = list("msWarnNoConv" = FALSE)) 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 = pdDiag(parent_0 + log_k_parent ~ 1), - start = mean_degparms(f, random = TRUE)) + random = pdLogChol(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) @@ -84,8 +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_silent(tmp <- update(m_nlme_mkin)) - expect_silent(tmp <- update(m_nlme_mmkin)) + expect_warning(tmp <- update(m_nlme_mmkin), "Iteration 1, LME step") 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)) |