aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2020-11-30 14:50:33 +0100
committerJohannes Ranke <jranke@uni-bremen.de>2020-11-30 14:50:33 +0100
commit78884beed74c18c99521b9ceeaa643e13cf94c06 (patch)
treeaf1ebfd975cac837a6bf15c86638299a67a0e017
parentb3bebc06061cc77b6d549f7b2760afe0b9489bf7 (diff)
Log-Cholesky parameterisation as default in nlme.mmkin
-rw-r--r--R/nlme.R1
-rw-r--r--R/nlme.mmkin.R47
-rw-r--r--R/summary.nlme.mmkin.R1
-rw-r--r--man/mkinmod.Rd6
-rw-r--r--man/nlme.Rd1
-rw-r--r--man/nlme.mmkin.Rd45
-rw-r--r--man/summary.nlme.mmkin.Rd1
-rw-r--r--test.log22
-rw-r--r--tests/testthat/test_nlme.R20
9 files changed, 86 insertions, 58 deletions
diff --git a/R/nlme.R b/R/nlme.R
index 8810fab3..9215aab0 100644
--- a/R/nlme.R
+++ b/R/nlme.R
@@ -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),
diff --git a/test.log b/test.log
index 67fd66c4..84a8748e 100644
--- a/test.log
+++ b/test.log
@@ -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))

Contact - Imprint