From 80ae27dcc4186f03e17164be74158454651474e7 Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Wed, 16 Nov 2022 14:50:12 +0100
Subject: Read in all data per default
---
R/read_spreadsheet.R | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
(limited to 'R')
diff --git a/R/read_spreadsheet.R b/R/read_spreadsheet.R
index a20af6db..7ad09c3e 100644
--- a/R/read_spreadsheet.R
+++ b/R/read_spreadsheet.R
@@ -37,7 +37,7 @@
#' and moisture normalisation factors in the sheet 'Datasets'?
#' @export
read_spreadsheet <- function(path, valid_datasets = "all",
- parent_only = TRUE, normalize = TRUE)
+ parent_only = FALSE, normalize = TRUE)
{
if (!requireNamespace("readxl", quietly = TRUE))
stop("Please install the readxl package to use this function")
--
cgit v1.2.3
From 9db8338d3a9240ed4685fcdd7aab9692031d5a04 Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Fri, 18 Nov 2022 05:56:25 +0100
Subject: Improve logLik.mkinfit to attach nobs attribute
The lack of that attribute lead to a failure to calculate the BIC in
test_AIC.R on R-devel from yesterday.
---
R/logLik.mkinfit.R | 1 +
log/test_dev.log | 61 ++++++++++++++++++++++++++++--------------------------
2 files changed, 33 insertions(+), 29 deletions(-)
(limited to 'R')
diff --git a/R/logLik.mkinfit.R b/R/logLik.mkinfit.R
index 7cc10234..abf8517c 100644
--- a/R/logLik.mkinfit.R
+++ b/R/logLik.mkinfit.R
@@ -37,6 +37,7 @@ logLik.mkinfit <- function(object, ...) {
val <- object$logLik
# Number of estimated parameters
attr(val, "df") <- length(object$bparms.optim) + length(object$errparms)
+ attr(val, "nobs") <- nobs(object)
class(val) <- "logLik"
return(val)
}
diff --git a/log/test_dev.log b/log/test_dev.log
index 24905a1a..527d28ed 100644
--- a/log/test_dev.log
+++ b/log/test_dev.log
@@ -1,54 +1,57 @@
-ℹ Loading mkin
-Loading required package: parallel
ℹ Testing mkin
✔ | F W S OK | Context
✔ | 5 | AIC calculation
-✔ | 5 | Analytical solutions for coupled models [14.6s]
+✔ | 5 | Analytical solutions for coupled models [3.2s]
✔ | 5 | Calculation of Akaike weights
-✔ | 2 | Export dataset for reading into CAKE
+✔ | 3 | Export dataset for reading into CAKE
✔ | 12 | Confidence intervals and p-values [1.0s]
-⠋ | 1 | Dimethenamid data from 2018
-✔ | 1 27 | Dimethenamid data from 2018 [116.1s]
+✔ | 1 12 | Dimethenamid data from 2018 [31.3s]
────────────────────────────────────────────────────────────────────────────────
-Skip (test_dmta.R:164:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data
+Skip ('test_dmta.R:98'): Different backends get consistent results for SFO-SFO3+, dimethenamid data
Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
────────────────────────────────────────────────────────────────────────────────
-✔ | 14 | Error model fitting [6.6s]
+✔ | 14 | Error model fitting [5.2s]
✔ | 5 | Time step normalisation
-✔ | 4 | Calculation of FOCUS chi2 error levels [0.8s]
-✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.5s]
-✔ | 4 | Test fitting the decline of metabolites from their maximum [0.6s]
-✔ | 1 | Fitting the logistic model [0.3s]
-⠋ | 11 | Nonlinear mixed-effects models
-✔ | 1 14 | Nonlinear mixed-effects models [1.3s]
+✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s]
+✔ | 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]
+✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [23.5s]
+✔ | 1 12 | Nonlinear mixed-effects models [0.3s]
────────────────────────────────────────────────────────────────────────────────
-Skip (test_mixed.R:68:3): saemix results are reproducible for biphasic fits
+Skip ('test_mixed.R:74'): saemix results are reproducible for biphasic fits
Reason: Fitting with saemix takes around 10 minutes when using deSolve
────────────────────────────────────────────────────────────────────────────────
✔ | 3 | Test dataset classes mkinds and mkindsg
-✔ | 10 | Special cases of mkinfit calls [0.6s]
-✔ | 3 | mkinfit features [1.1s]
+✔ | 10 | Special cases of mkinfit calls [0.5s]
+✔ | 3 | mkinfit features [0.7s]
✔ | 8 | mkinmod model generation and printing [0.2s]
✔ | 3 | Model predictions with mkinpredict [0.3s]
-✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.1s]
-✔ | 9 | Nonlinear mixed-effects models with nlme [8.7s]
-✔ | 16 | Plotting [1.4s]
+✔ | 9 | Multistart method for saem.mmkin models [35.2s]
+✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.5s]
+✔ | 9 | Nonlinear mixed-effects models with nlme [8.3s]
+✔ | 16 | Plotting [10.0s]
✔ | 4 | Residuals extracted from mkinfit models
-✔ | 23 | saemix parent models [28.4s]
-✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [12.0s]
-✔ | 7 | Fitting the SFORB model [16.9s]
+✔ | 1 36 | saemix parent models [63.3s]
+────────────────────────────────────────────────────────────────────────────────
+Skip ('test_saemix_parent.R:143'): We can also use mkin solution methods for saem
+Reason: This still takes almost 2.5 minutes although we do not solve ODEs
+────────────────────────────────────────────────────────────────────────────────
+✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.3s]
+✔ | 11 | Processing of residue series
+✔ | 10 | Fitting the SFORB model [3.5s]
✔ | 1 | Summaries of old mkinfit objects
-✔ | 4 | Summary [0.1s]
-✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [18.1s]
-✔ | 9 | Hypothesis tests [78.9s]
-✔ | 2 | tffm0
+✔ | 5 | Summary [0.2s]
+✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1s]
+✔ | 9 | Hypothesis tests [7.5s]
✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.0s]
══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 315.9 s
+Duration: 204.4 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)
+• This still takes almost 2.5 minutes although we do not solve ODEs (1)
-[ FAIL 0 | WARN 0 | SKIP 2 | PASS 240 ]
+[ FAIL 0 | WARN 0 | SKIP 3 | PASS 269 ]
--
cgit v1.2.3
From df0cff4b829f1abf62f037591a24a8019174dd0a Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Fri, 18 Nov 2022 08:37:40 +0100
Subject: Pass error.init to saemix_model, show in parplot
Due to an oversight, error.init was not really passed to saemix_model in
saem.mmkin. The new initial values were reverted to c(1, 1), in order to
avoid changing the test results. Initial values for error model
parameters are now shown in parplot.multistart.
---
R/parplot.R | 19 ++++++++---
R/saem.R | 3 +-
log/build.log | 1 +
log/test.log | 38 +++++++++++-----------
man/parplot.Rd | 5 +++
man/saem.Rd | 2 +-
.../multistart/parplot-for-biphasic-saemix-fit.svg | 2 ++
.../_snaps/multistart/parplot-for-sfo-fit.svg | 2 ++
tests/testthat/test_multistart.R | 8 ++---
9 files changed, 51 insertions(+), 29 deletions(-)
(limited to 'R')
diff --git a/R/parplot.R b/R/parplot.R
index 627a4eb8..98579779 100644
--- a/R/parplot.R
+++ b/R/parplot.R
@@ -4,6 +4,10 @@
#' either by the parameters of the run with the highest likelihood,
#' or by their medians as proposed in the paper by Duchesne et al. (2021).
#'
+#' Starting values of degradation model parameters and error model parameters
+#' are shown as green circles. The results obtained in the original run
+#' are shown as red circles.
+#'
#' @param object The [multistart] object
#' @param llmin The minimum likelihood of objects to be shown
#' @param scale By default, scale parameters using the best available fit.
@@ -32,7 +36,7 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
orig <- attr(object, "orig")
orig_parms <- parms(orig)
- start_parms <- orig$mean_dp_start
+ start_degparms <- orig$mean_dp_start
all_parms <- parms(object)
if (inherits(object, "multistart.saem.mmkin")) {
@@ -49,18 +53,18 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
par(las = 1)
if (orig$transformations == "mkin") {
- degparm_names_transformed <- names(start_parms)
+ degparm_names_transformed <- names(start_degparms)
degparm_index <- which(names(orig_parms) %in% degparm_names_transformed)
orig_parms[degparm_names_transformed] <- backtransform_odeparms(
orig_parms[degparm_names_transformed],
orig$mmkin[[1]]$mkinmod,
transform_rates = orig$mmkin[[1]]$transform_rates,
transform_fractions = orig$mmkin[[1]]$transform_fractions)
- start_parms <- backtransform_odeparms(start_parms,
+ start_degparms <- backtransform_odeparms(start_degparms,
orig$mmkin[[1]]$mkinmod,
transform_rates = orig$mmkin[[1]]$transform_rates,
transform_fractions = orig$mmkin[[1]]$transform_fractions)
- degparm_names <- names(start_parms)
+ degparm_names <- names(start_degparms)
names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index]))
@@ -72,6 +76,13 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
colnames(selected_parms)[1:length(degparm_names)] <- degparm_names
}
+ start_errparms <- orig$so@model@error.init
+ names(start_errparms) <- orig$so@model@name.sigma
+
+ start_omegaparms <- orig$so@model@omega.init
+
+ start_parms <- c(start_degparms, start_errparms)
+
scale <- match.arg(scale)
parm_scale <- switch(scale,
best = selected_parms[which.best(object[selected]), ],
diff --git a/R/saem.R b/R/saem.R
index 696ea0ee..c77ff70f 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -149,7 +149,7 @@ saem.mmkin <- function(object,
covariates = NULL,
covariate_models = NULL,
no_random_effect = NULL,
- error.init = c(3, 0.1),
+ error.init = c(1, 1),
nbiter.saemix = c(300, 100),
control = list(displayProgress = FALSE, print = FALSE,
nbiter.saemix = nbiter.saemix,
@@ -708,6 +708,7 @@ saemix_model <- function(object, solution_type = "auto",
covariance.model = covariance.model,
covariate.model = covariate.model,
omega.init = omega.init,
+ error.init = error.init,
...
)
attr(res, "mean_dp_start") <- degparms_optim
diff --git a/log/build.log b/log/build.log
index a56a64df..c4f9b8a2 100644
--- a/log/build.log
+++ b/log/build.log
@@ -6,3 +6,4 @@
* checking for LF line-endings in source and make files and shell scripts
* checking for empty or unneeded directories
* building ‘mkin_1.2.0.tar.gz’
+
diff --git a/log/test.log b/log/test.log
index b305bf58..af8e52fd 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,53 +1,53 @@
ℹ Testing mkin
✔ | F W S OK | Context
✔ | 5 | AIC calculation
-✔ | 5 | Analytical solutions for coupled models [3.5s]
+✔ | 5 | Analytical solutions for coupled models [3.2s]
✔ | 5 | Calculation of Akaike weights
✔ | 3 | Export dataset for reading into CAKE
-✔ | 12 | Confidence intervals and p-values [1.1s]
-✔ | 1 12 | Dimethenamid data from 2018 [34.5s]
+✔ | 12 | Confidence intervals and p-values [1.0s]
+✔ | 1 12 | Dimethenamid data from 2018 [31.4s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_dmta.R:98'): Different backends get consistent results for SFO-SFO3+, dimethenamid data
Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
────────────────────────────────────────────────────────────────────────────────
-✔ | 14 | Error model fitting [5.3s]
+✔ | 14 | Error model fitting [4.9s]
✔ | 5 | Time step normalisation
✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s]
✔ | 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.4s]
✔ | 1 | Fitting the logistic model [0.2s]
-✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [25.6s]
-✔ | 1 12 | Nonlinear mixed-effects models [0.4s]
+✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [23.9s]
+✔ | 1 12 | Nonlinear mixed-effects models [0.3s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_mixed.R:74'): saemix results are reproducible for biphasic fits
Reason: Fitting with saemix takes around 10 minutes when using deSolve
────────────────────────────────────────────────────────────────────────────────
✔ | 3 | Test dataset classes mkinds and mkindsg
-✔ | 10 | Special cases of mkinfit calls [0.6s]
+✔ | 10 | Special cases of mkinfit calls [0.5s]
✔ | 3 | mkinfit features [0.7s]
✔ | 8 | mkinmod model generation and printing [0.2s]
-✔ | 3 | Model predictions with mkinpredict [0.3s]
-✔ | 9 | Multistart method for saem.mmkin models [40.2s]
-✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.6s]
-✔ | 9 | Nonlinear mixed-effects models with nlme [9.5s]
-✔ | 16 | Plotting [10.5s]
+✔ | 3 | Model predictions with mkinpredict [0.4s]
+✔ | 9 | Multistart method for saem.mmkin models [37.0s]
+✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.5s]
+✔ | 9 | Nonlinear mixed-effects models with nlme [8.8s]
+✔ | 16 | Plotting [10.1s]
✔ | 4 | Residuals extracted from mkinfit models
-✔ | 1 36 | saemix parent models [69.8s]
+✔ | 1 36 | saemix parent models [66.0s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_saemix_parent.R:143'): We can also use mkin solution methods for saem
Reason: This still takes almost 2.5 minutes although we do not solve ODEs
────────────────────────────────────────────────────────────────────────────────
-✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.6s]
+✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
✔ | 11 | Processing of residue series
-✔ | 10 | Fitting the SFORB model [3.9s]
+✔ | 10 | Fitting the SFORB model [3.7s]
✔ | 1 | Summaries of old mkinfit objects
✔ | 5 | Summary [0.2s]
-✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3s]
-✔ | 9 | Hypothesis tests [8.5s]
-✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.3s]
+✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
+✔ | 9 | Hypothesis tests [8.0s]
+✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 225.9 s
+Duration: 211.3 s
── Skipped tests ──────────────────────────────────────────────────────────────
• Fitting this ODE model with saemix takes about 15 minutes on my system (1)
diff --git a/man/parplot.Rd b/man/parplot.Rd
index 37c5841d..ac9e02cf 100644
--- a/man/parplot.Rd
+++ b/man/parplot.Rd
@@ -35,6 +35,11 @@ Produces a boxplot with all parameters from the multiple runs, scaled
either by the parameters of the run with the highest likelihood,
or by their medians as proposed in the paper by Duchesne et al. (2021).
}
+\details{
+Starting values of degradation model parameters and error model parameters
+are shown as green circles. The results obtained in the original run
+are shown as red circles.
+}
\references{
Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical
identifiability in the frame of nonlinear mixed effects models: the example
diff --git a/man/saem.Rd b/man/saem.Rd
index 11463351..984d341b 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -24,7 +24,7 @@ saem(object, ...)
covariates = NULL,
covariate_models = NULL,
no_random_effect = NULL,
- error.init = c(3, 0.1),
+ error.init = c(1, 1),
nbiter.saemix = c(300, 100),
control = list(displayProgress = FALSE, print = FALSE, nbiter.saemix = nbiter.saemix,
save = FALSE, save.graphs = FALSE),
diff --git a/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg
index c0332fd5..7017908e 100644
--- a/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg
+++ b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg
@@ -173,6 +173,8 @@
+
+
diff --git a/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
index f3373901..18eb7fcc 100644
--- a/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
+++ b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
@@ -89,6 +89,8 @@
+
+
diff --git a/tests/testthat/test_multistart.R b/tests/testthat/test_multistart.R
index 98d3fb6d..502cee98 100644
--- a/tests/testthat/test_multistart.R
+++ b/tests/testthat/test_multistart.R
@@ -3,16 +3,16 @@ context("Multistart method for saem.mmkin models")
test_that("multistart works for saem.mmkin models", {
skip_on_cran() # Save CRAN time
set.seed(123456)
- saem_sfo_s_multi <- multistart(sfo_saem_1_reduced, n = 8, cores = n_cores,
- no_random_effect = "parent_0")
+ saem_sfo_s_multi <- multistart(sfo_saem_1_reduced, n = 8, cores = n_cores)
anova_sfo <- anova(sfo_saem_1,
sfo_saem_1_reduced,
best(saem_sfo_s_multi),
test = TRUE
)
# On winbuilder, sfo_saem_1 gives an AIC of 1310.8, while we get 1311.7
- # locally on Linux and Windows. The other, well-determined fits
- # both give 1309.7
+ # locally (using saemix 3.2, which likely makes the difference due to the
+ # error parameter patch) on Linux and Windows. The other, well-determined
+ # fits both give 1309.7.
expect_equal(round(anova_sfo, 1)["sfo_saem_1_reduced", "AIC"], 1309.7)
expect_equal(round(anova_sfo, 1)["best(saem_sfo_s_multi)", "AIC"], 1309.7)
expect_true(anova_sfo[3, "Pr(>Chisq)"] > 0.2) # Local: 1, CRAN: 0.4
--
cgit v1.2.3
From 5364f037a72863ef5ba81e14ba4417f68fd389f9 Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Fri, 18 Nov 2022 19:14:47 +0100
Subject: Make mixed model test data permanent to ensure reproducibility
To ensure that tests on different platforms work on the same data, the
mixed modelling test data previosly generated in
tests/testthat/setup_script.R were generated once using the script in
inst/dataset/generation/ds_mixed.R, and are now distributed with the
package.
---
R/ds_mixed.R | 17 ++
_pkgdown.yml | 3 +-
data/ds_mixed.rda | Bin 0 -> 6935 bytes
inst/dataset_generation/ds_mixed.R | 105 +++++++++++
log/check.log | 4 +-
log/test.log | 28 +--
man/ds_mixed.Rd | 24 +++
.../multistart/llhist-for-biphasic-saemix-fit.svg | 62 -------
.../_snaps/multistart/llhist-for-dfop-sfo-fit.svg | 62 +++++++
.../_snaps/multistart/llhist-for-sfo-fit.svg | 30 +--
.../multistart/parplot-for-biphasic-saemix-fit.svg | 201 ---------------------
.../_snaps/multistart/parplot-for-dfop-sfo-fit.svg | 201 +++++++++++++++++++++
.../_snaps/multistart/parplot-for-sfo-fit.svg | 90 ++++-----
tests/testthat/anova_sfo_saem.txt | 10 +-
tests/testthat/print_dfop_saem_1.txt | 23 +++
tests/testthat/print_dfop_saemix_1.txt | 23 ---
tests/testthat/print_fits_synth_const.txt | 2 +-
tests/testthat/print_mmkin_sfo_1_mixed.txt | 4 +-
tests/testthat/print_multistart_biphasic.txt | 4 -
tests/testthat/print_multistart_dfop_sfo.txt | 4 +
tests/testthat/print_sfo_saem_1_reduced.txt | 12 +-
tests/testthat/setup_script.R | 99 +---------
tests/testthat/summary_hfit_sfo_tc.txt | 26 +--
tests/testthat/summary_saem_biphasic_s.txt | 87 ---------
tests/testthat/summary_saem_dfop_sfo_s.txt | 87 +++++++++
tests/testthat/test_mixed.R | 32 ++--
tests/testthat/test_multistart.R | 30 ++-
tests/testthat/test_plot.R | 22 +--
tests/testthat/test_saemix_parent.R | 40 ++--
29 files changed, 697 insertions(+), 635 deletions(-)
create mode 100644 R/ds_mixed.R
create mode 100644 data/ds_mixed.rda
create mode 100644 inst/dataset_generation/ds_mixed.R
create mode 100644 man/ds_mixed.Rd
delete mode 100644 tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg
create mode 100644 tests/testthat/_snaps/multistart/llhist-for-dfop-sfo-fit.svg
delete mode 100644 tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg
create mode 100644 tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg
create mode 100644 tests/testthat/print_dfop_saem_1.txt
delete mode 100644 tests/testthat/print_dfop_saemix_1.txt
delete mode 100644 tests/testthat/print_multistart_biphasic.txt
create mode 100644 tests/testthat/print_multistart_dfop_sfo.txt
delete mode 100644 tests/testthat/summary_saem_biphasic_s.txt
create mode 100644 tests/testthat/summary_saem_dfop_sfo_s.txt
(limited to 'R')
diff --git a/R/ds_mixed.R b/R/ds_mixed.R
new file mode 100644
index 00000000..c5055712
--- /dev/null
+++ b/R/ds_mixed.R
@@ -0,0 +1,17 @@
+#' Synthetic data for hierarchical kinetic degradation models
+#'
+#' The R code used to create this data object is installed with this package in
+#' the 'dataset_generation' directory.
+#'
+#' @name ds_mixed
+#' @aliases ds_sfo ds_fomc ds_dfop ds_hs ds_dfop_sfo
+#' @examples
+#' \dontrun{
+#' sfo_mmkin <- mmkin("SFO", ds_sfo, quiet = TRUE, error_model = "tc", cores = 15)
+#' sfo_saem <- saem(sfo_mmkin, no_random_effect = "parent_0")
+#' plot(sfo_saem)
+#' }
+#'
+#' # This is the code used to generate the datasets
+#' cat(readLines(system.file("dataset_generation/ds_mixed.R", package = "mkin")), sep = "\n")
+NULL
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 8aa44ed4..98ed84f2 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -67,7 +67,7 @@ reference:
- parplot
- title: Datasets and known results
contents:
- - focus_soil_moisture
+ - ds_mixed
- D24_2014
- dimethenamid_2018
- FOCUS_2006_A
@@ -82,6 +82,7 @@ reference:
- synthetic_data_for_UBA_2014
- experimental_data_for_UBA_2019
- test_data_from_UBA_2014
+ - focus_soil_moisture
- mkinds
- mkindsg
- title: NAFTA guidance
diff --git a/data/ds_mixed.rda b/data/ds_mixed.rda
new file mode 100644
index 00000000..cf9f6463
Binary files /dev/null and b/data/ds_mixed.rda differ
diff --git a/inst/dataset_generation/ds_mixed.R b/inst/dataset_generation/ds_mixed.R
new file mode 100644
index 00000000..f2ae6e7e
--- /dev/null
+++ b/inst/dataset_generation/ds_mixed.R
@@ -0,0 +1,105 @@
+# Synthetic data for hierarchical kinetic models
+# Refactored version of the code previously in tests/testthat/setup_script.R
+# The number of datasets was 3 for FOMC, and 10 for HS in that script, now it
+# is always 15 for consistency
+
+library(mkin) # We use mkinmod and mkinpredict
+sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+n <- 15
+log_sd <- 0.3
+err_1 = list(const = 1, prop = 0.05)
+tc <- function(value) sigma_twocomp(value, err_1$const, err_1$prop)
+const <- function(value) 2
+
+set.seed(123456)
+SFO <- mkinmod(parent = mkinsub("SFO"))
+sfo_pop <- list(parent_0 = 100, k_parent = 0.03)
+sfo_parms <- as.matrix(data.frame(
+ k_parent = rlnorm(n, log(sfo_pop$k_parent), log_sd)))
+set.seed(123456)
+ds_sfo <- lapply(1:n, function(i) {
+ ds_mean <- mkinpredict(SFO, sfo_parms[i, ],
+ c(parent = sfo_pop$parent_0), sampling_times)
+ add_err(ds_mean, tc, n = 1)[[1]]
+})
+attr(ds_sfo, "pop") <- sfo_pop
+attr(ds_sfo, "parms") <- sfo_parms
+
+set.seed(123456)
+FOMC <- mkinmod(parent = mkinsub("FOMC"))
+fomc_pop <- list(parent_0 = 100, alpha = 2, beta = 8)
+fomc_parms <- as.matrix(data.frame(
+ alpha = rlnorm(n, log(fomc_pop$alpha), 0.4),
+ beta = rlnorm(n, log(fomc_pop$beta), 0.2)))
+set.seed(123456)
+ds_fomc <- lapply(1:n, function(i) {
+ ds_mean <- mkinpredict(FOMC, fomc_parms[i, ],
+ c(parent = fomc_pop$parent_0), sampling_times)
+ add_err(ds_mean, tc, n = 1)[[1]]
+})
+attr(ds_fomc, "pop") <- fomc_pop
+attr(ds_fomc, "parms") <- fomc_parms
+
+set.seed(123456)
+DFOP <- mkinmod(parent = mkinsub("DFOP"))
+dfop_pop <- list(parent_0 = 100, k1 = 0.06, k2 = 0.015, g = 0.4)
+dfop_parms <- as.matrix(data.frame(
+ k1 = rlnorm(n, log(dfop_pop$k1), log_sd),
+ k2 = rlnorm(n, log(dfop_pop$k2), log_sd),
+ g = plogis(rnorm(n, qlogis(dfop_pop$g), log_sd))))
+set.seed(123456)
+ds_dfop <- lapply(1:n, function(i) {
+ ds_mean <- mkinpredict(DFOP, dfop_parms[i, ],
+ c(parent = dfop_pop$parent_0), sampling_times)
+ add_err(ds_mean, tc, n = 1)[[1]]
+})
+attr(ds_dfop, "pop") <- dfop_pop
+attr(ds_dfop, "parms") <- dfop_parms
+
+set.seed(123456)
+HS <- mkinmod(parent = mkinsub("HS"))
+hs_pop <- list(parent_0 = 100, k1 = 0.08, k2 = 0.01, tb = 15)
+hs_parms <- as.matrix(data.frame(
+ k1 = rlnorm(n, log(hs_pop$k1), log_sd),
+ k2 = rlnorm(n, log(hs_pop$k2), log_sd),
+ tb = rlnorm(n, log(hs_pop$tb), 0.1)))
+set.seed(123456)
+ds_hs <- lapply(1:n, function(i) {
+ ds_mean <- mkinpredict(HS, hs_parms[i, ],
+ c(parent = hs_pop$parent_0), sampling_times)
+ add_err(ds_mean, const, n = 1)[[1]]
+})
+attr(ds_hs, "pop") <- hs_pop
+attr(ds_hs, "parms") <- hs_parms
+
+set.seed(123456)
+DFOP_SFO <- mkinmod(
+ parent = mkinsub("DFOP", "m1"),
+ m1 = mkinsub("SFO"),
+ quiet = TRUE)
+dfop_sfo_pop <- list(parent_0 = 100,
+ k_m1 = 0.007, f_parent_to_m1 = 0.5,
+ k1 = 0.1, k2 = 0.02, g = 0.5)
+dfop_sfo_parms <- as.matrix(data.frame(
+ k1 = rlnorm(n, log(dfop_sfo_pop$k1), log_sd),
+ k2 = rlnorm(n, log(dfop_sfo_pop$k2), log_sd),
+ g = plogis(rnorm(n, qlogis(dfop_sfo_pop$g), log_sd)),
+ f_parent_to_m1 = plogis(rnorm(n,
+ qlogis(dfop_sfo_pop$f_parent_to_m1), log_sd)),
+ k_m1 = rlnorm(n, log(dfop_sfo_pop$k_m1), log_sd)))
+ds_dfop_sfo_mean <- lapply(1:n,
+ function(i) {
+ mkinpredict(DFOP_SFO, dfop_sfo_parms[i, ],
+ c(parent = dfop_sfo_pop$parent_0, m1 = 0), sampling_times)
+ }
+)
+set.seed(123456)
+ds_dfop_sfo <- lapply(ds_dfop_sfo_mean, function(ds) {
+ add_err(ds,
+ sdfunc = function(value) sqrt(err_1$const^2 + value^2 * err_1$prop^2),
+ n = 1, secondary = "m1")[[1]]
+})
+attr(ds_dfop_sfo, "pop") <- dfop_sfo_pop
+attr(ds_dfop_sfo, "parms") <- dfop_sfo_parms
+
+#save(ds_sfo, ds_fomc, ds_dfop, ds_hs, ds_dfop_sfo, file = "data/ds_mixed.rda", version = 2)
diff --git a/log/check.log b/log/check.log
index 7aa4610c..31fc31eb 100644
--- a/log/check.log
+++ b/log/check.log
@@ -5,7 +5,7 @@
* using options ‘--no-tests --as-cran’
* checking for file ‘mkin/DESCRIPTION’ ... OK
* checking extension type ... Package
-* this is package ‘mkin’ version ‘1.2.0’
+* this is package ‘mkin’ version ‘1.2.1’
* package encoding: UTF-8
* checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
Maintainer: ‘Johannes Ranke ’
@@ -41,7 +41,7 @@ Maintainer: ‘Johannes Ranke ’
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
-* checking R code for possible problems ... [18s/18s] OK
+* checking R code for possible problems ... [17s/17s] OK
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd line widths ... OK
diff --git a/log/test.log b/log/test.log
index 87d24690..d1de270e 100644
--- a/log/test.log
+++ b/log/test.log
@@ -5,19 +5,19 @@
✔ | 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.5s]
+✔ | 1 12 | Dimethenamid data from 2018 [32.6s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_dmta.R:98'): Different backends get consistent results for SFO-SFO3+, dimethenamid data
Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
────────────────────────────────────────────────────────────────────────────────
-✔ | 14 | Error model fitting [5.0s]
+✔ | 14 | Error model fitting [5.1s]
✔ | 5 | Time step normalisation
✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s]
✔ | 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.4s]
✔ | 1 | Fitting the logistic model [0.2s]
-✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [24.4s]
-✔ | 1 12 | Nonlinear mixed-effects models [0.3s]
+✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [24.6s]
+✔ | 1 13 | Nonlinear mixed-effects models [0.4s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_mixed.R:79'): saemix results are reproducible for biphasic fits
Reason: Fitting with saemix takes around 10 minutes when using deSolve
@@ -26,32 +26,32 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
✔ | 10 | Special cases of mkinfit calls [0.5s]
✔ | 3 | mkinfit features [0.7s]
✔ | 8 | mkinmod model generation and printing [0.2s]
-✔ | 3 | Model predictions with mkinpredict [0.3s]
-✔ | 9 | Multistart method for saem.mmkin models [36.2s]
+✔ | 3 | Model predictions with mkinpredict [0.4s]
+✔ | 9 | Multistart method for saem.mmkin models [37.3s]
✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.5s]
-✔ | 9 | Nonlinear mixed-effects models with nlme [8.8s]
-✔ | 16 | Plotting [10.0s]
+✔ | 9 | Nonlinear mixed-effects models with nlme [9.2s]
+✔ | 16 | Plotting [10.1s]
✔ | 4 | Residuals extracted from mkinfit models
-✔ | 1 36 | saemix parent models [65.8s]
+✔ | 1 36 | saemix parent models [71.4s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_saemix_parent.R:143'): We can also use mkin solution methods for saem
Reason: This still takes almost 2.5 minutes although we do not solve ODEs
────────────────────────────────────────────────────────────────────────────────
✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
✔ | 11 | Processing of residue series
-✔ | 10 | Fitting the SFORB model [3.7s]
+✔ | 10 | Fitting the SFORB model [3.9s]
✔ | 1 | Summaries of old mkinfit objects
✔ | 5 | Summary [0.2s]
-✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
-✔ | 9 | Hypothesis tests [8.0s]
+✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3s]
+✔ | 9 | Hypothesis tests [8.3s]
✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 210.4 s
+Duration: 220.0 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)
• This still takes almost 2.5 minutes although we do not solve ODEs (1)
-[ FAIL 0 | WARN 0 | SKIP 3 | PASS 269 ]
+[ FAIL 0 | WARN 0 | SKIP 3 | PASS 270 ]
diff --git a/man/ds_mixed.Rd b/man/ds_mixed.Rd
new file mode 100644
index 00000000..227b8e7f
--- /dev/null
+++ b/man/ds_mixed.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ds_mixed.R
+\name{ds_mixed}
+\alias{ds_mixed}
+\alias{ds_sfo}
+\alias{ds_fomc}
+\alias{ds_dfop}
+\alias{ds_hs}
+\alias{ds_dfop_sfo}
+\title{Synthetic data for hierarchical kinetic degradation models}
+\description{
+The R code used to create this data object is installed with this package in
+the 'dataset_generation' directory.
+}
+\examples{
+\dontrun{
+ sfo_mmkin <- mmkin("SFO", ds_sfo, quiet = TRUE, error_model = "tc", cores = 15)
+ sfo_saem <- saem(sfo_mmkin, no_random_effect = "parent_0")
+ plot(sfo_saem)
+}
+
+# This is the code used to generate the datasets
+cat(readLines(system.file("dataset_generation/ds_mixed.R", package = "mkin")), sep = "\n")
+}
diff --git a/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg
deleted file mode 100644
index 6015aed8..00000000
--- a/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg
+++ /dev/null
@@ -1,62 +0,0 @@
-
-
diff --git a/tests/testthat/_snaps/multistart/llhist-for-dfop-sfo-fit.svg b/tests/testthat/_snaps/multistart/llhist-for-dfop-sfo-fit.svg
new file mode 100644
index 00000000..6015aed8
--- /dev/null
+++ b/tests/testthat/_snaps/multistart/llhist-for-dfop-sfo-fit.svg
@@ -0,0 +1,62 @@
+
+
diff --git a/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg
index 98513d06..028c69de 100644
--- a/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg
+++ b/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg
@@ -25,20 +25,22 @@
--649.836
--649.834
--649.832
--649.830
--649.828
+-646.124
+-646.123
+-646.122
+-646.121
+-646.120
-
-
+
+
+0
-1
-2
-3
+1
+2
+3
+4
@@ -46,11 +48,11 @@
-
+
-
-
-
+
+
+original fit
diff --git a/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg
deleted file mode 100644
index 7017908e..00000000
--- a/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg
+++ /dev/null
@@ -1,201 +0,0 @@
-
-
diff --git a/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg
new file mode 100644
index 00000000..7017908e
--- /dev/null
+++ b/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg
@@ -0,0 +1,201 @@
+
+
diff --git a/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
index 18eb7fcc..a47a585a 100644
--- a/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
+++ b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
@@ -25,42 +25,44 @@
-
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -87,15 +89,15 @@
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
diff --git a/tests/testthat/anova_sfo_saem.txt b/tests/testthat/anova_sfo_saem.txt
index 9e4bf71f..0ccd6d5c 100644
--- a/tests/testthat/anova_sfo_saem.txt
+++ b/tests/testthat/anova_sfo_saem.txt
@@ -1,7 +1,7 @@
-Data: 262 observations of 1 variable(s) grouped in 15 datasets
+Data: 263 observations of 1 variable(s) grouped in 15 datasets
npar AIC BIC Lik Chisq Df Pr(>Chisq)
-sfo_saem_1_reduced 5 1310 1313 -650
-sfo_saem_1_reduced_mkin 5 1310 1313 -650 0 0
-sfo_saem_1 6 1312 1316 -650 0 1 1
-sfo_saem_1_mkin 6 1312 1316 -650 0 0
+sfo_saem_1_reduced 5 1302 1306 -646
+sfo_saem_1_reduced_mkin 5 1302 1306 -646 0 0
+sfo_saem_1 6 1304 1308 -646 0 1 1
+sfo_saem_1_mkin 6 1303 1308 -646 1 0
diff --git a/tests/testthat/print_dfop_saem_1.txt b/tests/testthat/print_dfop_saem_1.txt
new file mode 100644
index 00000000..bdc40065
--- /dev/null
+++ b/tests/testthat/print_dfop_saem_1.txt
@@ -0,0 +1,23 @@
+Kinetic nonlinear mixed-effects model fit by SAEM
+Structural model:
+d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+ time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+ * parent
+
+Data:
+270 observations of 1 variable(s) grouped in 15 datasets
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 1409 1415 -696
+
+Fitted parameters:
+ estimate lower upper
+parent_0 99.92 98.77 101.06
+log_k1 -2.72 -2.95 -2.50
+log_k2 -4.14 -4.27 -4.01
+g_qlogis -0.35 -0.53 -0.16
+a.1 0.92 0.68 1.16
+b.1 0.05 0.04 0.06
+SD.log_k1 0.37 0.23 0.51
+SD.log_k2 0.23 0.14 0.31
diff --git a/tests/testthat/print_dfop_saemix_1.txt b/tests/testthat/print_dfop_saemix_1.txt
deleted file mode 100644
index 1d399a52..00000000
--- a/tests/testthat/print_dfop_saemix_1.txt
+++ /dev/null
@@ -1,23 +0,0 @@
-Kinetic nonlinear mixed-effects model fit by SAEM
-Structural model:
-d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
- time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
- * parent
-
-Data:
-270 observations of 1 variable(s) grouped in 15 datasets
-
-Likelihood computed by importance sampling
- AIC BIC logLik
- 1409 1415 -696
-
-Fitted parameters:
- estimate lower upper
-parent_0 100.00 99.00 100.00
-log_k1 -2.70 -3.00 -2.50
-log_k2 -4.10 -4.30 -4.00
-g_qlogis -0.35 -0.53 -0.16
-a.1 0.92 0.68 1.20
-b.1 0.05 0.04 0.06
-SD.log_k1 0.37 0.23 0.51
-SD.log_k2 0.23 0.14 0.31
diff --git a/tests/testthat/print_fits_synth_const.txt b/tests/testthat/print_fits_synth_const.txt
index 2ea1f133..b4bbe6ca 100644
--- a/tests/testthat/print_fits_synth_const.txt
+++ b/tests/testthat/print_fits_synth_const.txt
@@ -4,7 +4,7 @@ Status of individual fits:
dataset
model 1 2 3 4 5 6
SFO OK OK OK OK OK OK
- FOMC C C OK OK OK OK
+ FOMC C OK OK OK OK C
C: Optimisation did not converge:
false convergence (8)
diff --git a/tests/testthat/print_mmkin_sfo_1_mixed.txt b/tests/testthat/print_mmkin_sfo_1_mixed.txt
index 33e5bf5c..c12cfe2b 100644
--- a/tests/testthat/print_mmkin_sfo_1_mixed.txt
+++ b/tests/testthat/print_mmkin_sfo_1_mixed.txt
@@ -3,7 +3,7 @@ Structural model:
d_parent/dt = - k_parent * parent
Data:
-262 observations of 1 variable(s) grouped in 15 datasets
+263 observations of 1 variable(s) grouped in 15 datasets
object
Status of individual fits:
@@ -16,4 +16,4 @@ OK: No warnings
Mean fitted parameters:
parent_0 log_k_parent
- 99.9 -3.3
+ 100.0 -3.4
diff --git a/tests/testthat/print_multistart_biphasic.txt b/tests/testthat/print_multistart_biphasic.txt
deleted file mode 100644
index b4344f22..00000000
--- a/tests/testthat/print_multistart_biphasic.txt
+++ /dev/null
@@ -1,4 +0,0 @@
- object with 8 fits:
-OK
- 8
-OK: Fit terminated successfully
diff --git a/tests/testthat/print_multistart_dfop_sfo.txt b/tests/testthat/print_multistart_dfop_sfo.txt
new file mode 100644
index 00000000..b4344f22
--- /dev/null
+++ b/tests/testthat/print_multistart_dfop_sfo.txt
@@ -0,0 +1,4 @@
+ object with 8 fits:
+OK
+ 8
+OK: Fit terminated successfully
diff --git a/tests/testthat/print_sfo_saem_1_reduced.txt b/tests/testthat/print_sfo_saem_1_reduced.txt
index bac8848e..1c7fb588 100644
--- a/tests/testthat/print_sfo_saem_1_reduced.txt
+++ b/tests/testthat/print_sfo_saem_1_reduced.txt
@@ -3,16 +3,16 @@ Structural model:
d_parent/dt = - k_parent * parent
Data:
-262 observations of 1 variable(s) grouped in 15 datasets
+263 observations of 1 variable(s) grouped in 15 datasets
Likelihood computed by importance sampling
AIC BIC logLik
- 1310 1313 -650
+ 1302 1306 -646
Fitted parameters:
estimate lower upper
-parent_0 1e+02 99.08 1e+02
-k_parent 4e-02 0.03 4e-02
-a.1 9e-01 0.75 1e+00
+parent_0 1e+02 99.03 1e+02
+k_parent 3e-02 0.03 4e-02
+a.1 9e-01 0.71 1e+00
b.1 5e-02 0.04 5e-02
-SD.k_parent 3e-01 0.20 4e-01
+SD.k_parent 2e-01 0.14 3e-01
diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R
index 362038c3..c554800d 100644
--- a/tests/testthat/setup_script.R
+++ b/tests/testthat/setup_script.R
@@ -81,112 +81,27 @@ fit_obs_1 <- mkinfit(m_synth_SFO_lin, SFO_lin_a, error_model = "obs", quiet = TR
fit_tc_1 <- mkinfit(m_synth_SFO_lin, SFO_lin_a, error_model = "tc", quiet = TRUE,
error_model_algorithm = "threestep")
-# Mixed models data and fits
-sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
-n <- n_biphasic <- 15
-log_sd <- 0.3
-err_1 = list(const = 1, prop = 0.05)
-tc <- function(value) sigma_twocomp(value, err_1$const, err_1$prop)
-const <- function(value) 2
-
-set.seed(123456)
-SFO <- mkinmod(parent = mkinsub("SFO"))
-k_parent = rlnorm(n, log(0.03), log_sd)
-set.seed(123456)
-ds_sfo <- lapply(1:n, function(i) {
- ds_mean <- mkinpredict(SFO, c(k_parent = k_parent[i]),
- c(parent = 100), sampling_times)
- add_err(ds_mean, tc, n = 1)[[1]]
-})
-
-set.seed(123456)
-FOMC <- mkinmod(parent = mkinsub("FOMC"))
-fomc_pop <- list(parent_0 = 100, alpha = 2, beta = 8)
-fomc_parms <- as.matrix(data.frame(
- alpha = rlnorm(n, log(fomc_pop$alpha), 0.4),
- beta = rlnorm(n, log(fomc_pop$beta), 0.2)))
-set.seed(123456)
-ds_fomc <- lapply(1:3, function(i) {
- ds_mean <- mkinpredict(FOMC, fomc_parms[i, ],
- c(parent = 100), sampling_times)
- add_err(ds_mean, tc, n = 1)[[1]]
-})
-
-set.seed(123456)
-DFOP <- mkinmod(parent = mkinsub("DFOP"))
-dfop_pop <- list(parent_0 = 100, k1 = 0.06, k2 = 0.015, g = 0.4)
-dfop_parms <- as.matrix(data.frame(
- k1 = rlnorm(n, log(dfop_pop$k1), log_sd),
- k2 = rlnorm(n, log(dfop_pop$k2), log_sd),
- g = plogis(rnorm(n, qlogis(dfop_pop$g), log_sd))))
-set.seed(123456)
-ds_dfop <- lapply(1:n, function(i) {
- ds_mean <- mkinpredict(DFOP, dfop_parms[i, ],
- c(parent = dfop_pop$parent_0), sampling_times)
- add_err(ds_mean, tc, n = 1)[[1]]
-})
-
-set.seed(123456)
-HS <- mkinmod(parent = mkinsub("HS"))
-hs_pop <- list(parent_0 = 100, k1 = 0.08, k2 = 0.01, tb = 15)
-hs_parms <- as.matrix(data.frame(
- k1 = rlnorm(n, log(hs_pop$k1), log_sd),
- k2 = rlnorm(n, log(hs_pop$k2), log_sd),
- tb = rlnorm(n, log(hs_pop$tb), 0.1)))
-set.seed(123456)
-ds_hs <- lapply(1:10, function(i) {
- ds_mean <- mkinpredict(HS, hs_parms[i, ],
- c(parent = hs_pop$parent_0), sampling_times)
- add_err(ds_mean, const, n = 1)[[1]]
-})
-
-set.seed(123456)
-DFOP_SFO <- mkinmod(
- parent = mkinsub("DFOP", "m1"),
- m1 = mkinsub("SFO"),
- quiet = TRUE)
-dfop_sfo_pop <- list(parent_0 = 100,
- k_m1 = 0.007, f_parent_to_m1 = 0.5,
- k1 = 0.1, k2 = 0.02, g = 0.5)
-syn_biphasic_parms <- as.matrix(data.frame(
- k1 = rlnorm(n_biphasic, log(dfop_sfo_pop$k1), log_sd),
- k2 = rlnorm(n_biphasic, log(dfop_sfo_pop$k2), log_sd),
- g = plogis(rnorm(n_biphasic, qlogis(dfop_sfo_pop$g), log_sd)),
- f_parent_to_m1 = plogis(rnorm(n_biphasic,
- qlogis(dfop_sfo_pop$f_parent_to_m1), log_sd)),
- k_m1 = rlnorm(n_biphasic, log(dfop_sfo_pop$k_m1), log_sd)))
-ds_biphasic_mean <- lapply(1:n_biphasic,
- function(i) {
- mkinpredict(DFOP_SFO, syn_biphasic_parms[i, ],
- c(parent = 100, m1 = 0), sampling_times)
- }
-)
-set.seed(123456)
-ds_biphasic <- lapply(ds_biphasic_mean, function(ds) {
- add_err(ds,
- sdfunc = function(value) sqrt(err_1$const^2 + value^2 * err_1$prop^2),
- n = 1, secondary = "m1")[[1]]
-})
-
# Mixed model fits
mmkin_sfo_1 <- mmkin("SFO", ds_sfo, quiet = TRUE, error_model = "tc", cores = n_cores)
mmkin_dfop_1 <- mmkin("DFOP", ds_dfop, quiet = TRUE, cores = n_cores,
error_model = "tc")
-mmkin_biphasic <- mmkin(list("DFOP-SFO" = DFOP_SFO), ds_biphasic, quiet = TRUE, cores = n_cores,
+DFOP_SFO <- mkinmod(parent = mkinsub("DFOP", "m1"),
+ m1 = mkinsub("SFO"), quiet = TRUE)
+mmkin_dfop_sfo <- mmkin(list("DFOP-SFO" = DFOP_SFO), ds_dfop_sfo, quiet = TRUE, cores = n_cores,
control = list(eval.max = 500, iter.max = 400),
error_model = "tc")
# nlme
dfop_nlme_1 <- suppressWarnings(nlme(mmkin_dfop_1))
-nlme_biphasic <- suppressWarnings(nlme(mmkin_biphasic))
+nlme_dfop_sfo <- suppressWarnings(nlme(mmkin_dfop_sfo))
# saemix
sfo_saem_1 <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "saemix")
sfo_saem_1_reduced <- update(sfo_saem_1, no_random_effect = "parent_0")
-dfop_saemix_1 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "mkin",
+dfop_saem_1 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "mkin",
no_random_effect = c("parent_0", "g_qlogis"))
-saem_biphasic_m <- saem(mmkin_biphasic, transformations = "mkin", quiet = TRUE)
-saem_biphasic_s <- saem(mmkin_biphasic, transformations = "saemix", quiet = TRUE)
+saem_dfop_sfo_m <- saem(mmkin_dfop_sfo, transformations = "mkin", quiet = TRUE)
+saem_dfop_sfo_s <- saem(mmkin_dfop_sfo, transformations = "saemix", quiet = TRUE)
diff --git a/tests/testthat/summary_hfit_sfo_tc.txt b/tests/testthat/summary_hfit_sfo_tc.txt
index 41743091..bb5bf6fb 100644
--- a/tests/testthat/summary_hfit_sfo_tc.txt
+++ b/tests/testthat/summary_hfit_sfo_tc.txt
@@ -8,7 +8,7 @@ Equations:
d_parent/dt = - k_parent * parent
Data:
-106 observations of 1 variable(s) grouped in 6 datasets
+104 observations of 1 variable(s) grouped in 6 datasets
Model predictions using solution type analytical
@@ -28,15 +28,15 @@ Results:
Likelihood computed by importance sampling
AIC BIC logLik
- 533 531 -261
+ 524 523 -257
Optimised parameters:
- est. lower upper
-parent_0 101.02 99.58 102.46
-log_k_parent -3.32 -3.53 -3.11
-a.1 0.91 0.64 1.17
-b.1 0.05 0.04 0.06
-SD.log_k_parent 0.27 0.11 0.42
+ est. lower upper
+parent_0 100.68 99.27 102.08
+log_k_parent -3.38 -3.55 -3.21
+a.1 0.87 0.59 1.14
+b.1 0.05 0.04 0.06
+SD.log_k_parent 0.21 0.09 0.33
Correlation:
pr_0
@@ -44,18 +44,18 @@ log_k_parent 0.1
Random effects:
est. lower upper
-SD.log_k_parent 0.3 0.1 0.4
+SD.log_k_parent 0.2 0.09 0.3
Variance model:
est. lower upper
-a.1 0.91 0.64 1.17
+a.1 0.87 0.59 1.14
b.1 0.05 0.04 0.06
Backtransformed parameters:
est. lower upper
-parent_0 1e+02 1e+02 1e+02
-k_parent 4e-02 3e-02 4e-02
+parent_0 1e+02 99.27 1e+02
+k_parent 3e-02 0.03 4e-02
Estimated disappearance times:
DT50 DT90
-parent 19 64
+parent 20 68
diff --git a/tests/testthat/summary_saem_biphasic_s.txt b/tests/testthat/summary_saem_biphasic_s.txt
deleted file mode 100644
index 7c337843..00000000
--- a/tests/testthat/summary_saem_biphasic_s.txt
+++ /dev/null
@@ -1,87 +0,0 @@
-saemix version used for fitting: Dummy 0.0 for testing
-mkin version used for pre-fitting: Dummy 0.0 for testing
-R version used for fitting: Dummy R version for testing
-Date of fit: Dummy date for testing
-Date of summary: Dummy date for testing
-
-Equations:
-d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
- time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
- * parent
-d_m1/dt = + f_parent_to_m1 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g)
- * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) *
- exp(-k2 * time))) * parent - k_m1 * m1
-
-Data:
-510 observations of 2 variable(s) grouped in 15 datasets
-
-Model predictions using solution type analytical
-
-Fitted in test time 0 s
-Using 300, 100 iterations and 4 chains
-
-Variance model: Two-component variance function
-
-Mean of starting values for individual parameters:
- parent_0 k_m1 f_parent_to_m1 k1 k2
- 1e+02 7e-03 5e-01 1e-01 2e-02
- g
- 5e-01
-
-Fixed degradation parameter values:
-None
-
-Results:
-
-Likelihood computed by importance sampling
- AIC BIC logLik
- 2334 2344 -1153
-
-Optimised parameters:
- est. lower upper
-parent_0 1e+02 1e+02 1e+02
-k_m1 7e-03 6e-03 7e-03
-f_parent_to_m1 5e-01 4e-01 5e-01
-k1 1e-01 9e-02 1e-01
-k2 2e-02 2e-02 3e-02
-g 5e-01 5e-01 5e-01
-a.1 9e-01 8e-01 1e+00
-b.1 5e-02 5e-02 6e-02
-SD.parent_0 3e-02 -5e+01 5e+01
-SD.k_m1 2e-01 1e-01 3e-01
-SD.f_parent_to_m1 3e-01 2e-01 4e-01
-SD.k1 4e-01 2e-01 5e-01
-SD.k2 3e-01 2e-01 5e-01
-SD.g 2e-01 6e-02 4e-01
-
-Correlation:
- pr_0 k_m1 f___ k1 k2
-k_m1 -0.2
-f_parent_to_m1 -0.3 0.1
-k1 0.1 0.0 0.0
-k2 0.0 0.0 0.0 0.1
-g 0.1 -0.1 0.0 -0.2 -0.2
-
-Random effects:
- est. lower upper
-SD.parent_0 0.03 -49.24 49.3
-SD.k_m1 0.23 0.13 0.3
-SD.f_parent_to_m1 0.30 0.19 0.4
-SD.k1 0.40 0.25 0.5
-SD.k2 0.34 0.21 0.5
-SD.g 0.21 0.06 0.4
-
-Variance model:
- est. lower upper
-a.1 0.93 0.79 1.06
-b.1 0.05 0.05 0.06
-
-Resulting formation fractions:
- ff
-parent_m1 0.5
-parent_sink 0.5
-
-Estimated disappearance times:
- DT50 DT90 DT50back DT50_k1 DT50_k2
-parent 13 73 22 6 32
-m1 105 348 NA NA NA
diff --git a/tests/testthat/summary_saem_dfop_sfo_s.txt b/tests/testthat/summary_saem_dfop_sfo_s.txt
new file mode 100644
index 00000000..7c337843
--- /dev/null
+++ b/tests/testthat/summary_saem_dfop_sfo_s.txt
@@ -0,0 +1,87 @@
+saemix version used for fitting: Dummy 0.0 for testing
+mkin version used for pre-fitting: Dummy 0.0 for testing
+R version used for fitting: Dummy R version for testing
+Date of fit: Dummy date for testing
+Date of summary: Dummy date for testing
+
+Equations:
+d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+ time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+ * parent
+d_m1/dt = + f_parent_to_m1 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g)
+ * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) *
+ exp(-k2 * time))) * parent - k_m1 * m1
+
+Data:
+510 observations of 2 variable(s) grouped in 15 datasets
+
+Model predictions using solution type analytical
+
+Fitted in test time 0 s
+Using 300, 100 iterations and 4 chains
+
+Variance model: Two-component variance function
+
+Mean of starting values for individual parameters:
+ parent_0 k_m1 f_parent_to_m1 k1 k2
+ 1e+02 7e-03 5e-01 1e-01 2e-02
+ g
+ 5e-01
+
+Fixed degradation parameter values:
+None
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 2334 2344 -1153
+
+Optimised parameters:
+ est. lower upper
+parent_0 1e+02 1e+02 1e+02
+k_m1 7e-03 6e-03 7e-03
+f_parent_to_m1 5e-01 4e-01 5e-01
+k1 1e-01 9e-02 1e-01
+k2 2e-02 2e-02 3e-02
+g 5e-01 5e-01 5e-01
+a.1 9e-01 8e-01 1e+00
+b.1 5e-02 5e-02 6e-02
+SD.parent_0 3e-02 -5e+01 5e+01
+SD.k_m1 2e-01 1e-01 3e-01
+SD.f_parent_to_m1 3e-01 2e-01 4e-01
+SD.k1 4e-01 2e-01 5e-01
+SD.k2 3e-01 2e-01 5e-01
+SD.g 2e-01 6e-02 4e-01
+
+Correlation:
+ pr_0 k_m1 f___ k1 k2
+k_m1 -0.2
+f_parent_to_m1 -0.3 0.1
+k1 0.1 0.0 0.0
+k2 0.0 0.0 0.0 0.1
+g 0.1 -0.1 0.0 -0.2 -0.2
+
+Random effects:
+ est. lower upper
+SD.parent_0 0.03 -49.24 49.3
+SD.k_m1 0.23 0.13 0.3
+SD.f_parent_to_m1 0.30 0.19 0.4
+SD.k1 0.40 0.25 0.5
+SD.k2 0.34 0.21 0.5
+SD.g 0.21 0.06 0.4
+
+Variance model:
+ est. lower upper
+a.1 0.93 0.79 1.06
+b.1 0.05 0.05 0.06
+
+Resulting formation fractions:
+ ff
+parent_m1 0.5
+parent_sink 0.5
+
+Estimated disappearance times:
+ DT50 DT90 DT50back DT50_k1 DT50_k2
+parent 13 73 22 6 32
+m1 105 348 NA NA NA
diff --git a/tests/testthat/test_mixed.R b/tests/testthat/test_mixed.R
index 2d53c6dd..ab8dfc27 100644
--- a/tests/testthat/test_mixed.R
+++ b/tests/testthat/test_mixed.R
@@ -1,22 +1,23 @@
context("Nonlinear mixed-effects models")
+
# Round error model parameters as they are not rounded in print methods
dfop_nlme_1$modelStruct$varStruct$const <-
signif(dfop_nlme_1$modelStruct$varStruct$const, 3)
dfop_nlme_1$modelStruct$varStruct$prop <-
signif(dfop_nlme_1$modelStruct$varStruct$prop, 4)
+dfop_sfo_pop <- attr(ds_dfop_sfo, "pop")
+
test_that("Print methods work", {
expect_known_output(print(fits[, 2:3], digits = 2), "print_mmkin_parent.txt")
expect_known_output(print(mixed(mmkin_sfo_1), digits = 2), "print_mmkin_sfo_1_mixed.txt")
expect_known_output(print(dfop_nlme_1, digits = 1), "print_dfop_nlme_1.txt")
+ expect_known_output(print(sfo_saem_1_reduced, digits = 1), "print_sfo_saem_1_reduced.txt")
- # In order to address the platform dependence of the results, we round to two
- # significant digits before printing
- dfop_saemix_1_print <- dfop_saemix_1
- dfop_saemix_1_print$so@results@conf.int[c("estimate", "lower", "upper")] <-
- signif(dfop_saemix_1_print$so@results@conf.int[c("estimate", "lower", "upper")], 2)
- expect_known_output(print(dfop_saemix_1_print, digits = 1), "print_dfop_saemix_1.txt")
+ skip_on_cran() # The following test is platform dependent and fails on
+ # win-builder with current (18 Nov 2022) R-devel and on the Fedora CRAN check systems
+ expect_known_output(print(dfop_saem_1, digits = 1), "print_dfop_saem_1.txt")
})
test_that("nlme results are reproducible to some degree", {
@@ -36,17 +37,16 @@ test_that("nlme results are reproducible to some degree", {
# k1 and k2 just fail the first test (lower bound of the ci), so we need to exclude it
dfop_no_k1_k2 <- c("parent_0", "k_m1", "f_parent_to_m1", "g")
dfop_sfo_pop_no_k1_k2 <- as.numeric(dfop_sfo_pop[dfop_no_k1_k2])
- dfop_sfo_pop <- as.numeric(dfop_sfo_pop) # to remove names
- ci_dfop_sfo_n <- summary(nlme_biphasic)$confint_back
+ ci_dfop_sfo_n <- summary(nlme_dfop_sfo)$confint_back
expect_true(all(ci_dfop_sfo_n[dfop_no_k1_k2, "lower"] < dfop_sfo_pop_no_k1_k2))
- expect_true(all(ci_dfop_sfo_n[, "upper"] > dfop_sfo_pop))
+ expect_true(all(ci_dfop_sfo_n[, "upper"] > as.numeric(dfop_sfo_pop)))
})
test_that("saemix results are reproducible for biphasic fits", {
- test_summary <- summary(saem_biphasic_s)
+ test_summary <- summary(saem_dfop_sfo_s)
test_summary$saemixversion <- "Dummy 0.0 for testing"
test_summary$mkinversion <- "Dummy 0.0 for testing"
test_summary$Rversion <- "Dummy R version for testing"
@@ -54,33 +54,33 @@ test_that("saemix results are reproducible for biphasic fits", {
test_summary$date.summary <- "Dummy date for testing"
test_summary$time <- c(elapsed = "test time 0")
- expect_known_output(print(test_summary, digits = 1), "summary_saem_biphasic_s.txt")
+ expect_known_output(print(test_summary, digits = 1), "summary_saem_dfop_sfo_s.txt")
dfop_sfo_pop <- as.numeric(dfop_sfo_pop)
no_k1 <- c(1, 2, 3, 5, 6)
no_k2 <- c(1, 2, 3, 4, 6)
no_k1_k2 <- c(1, 2, 3, 6)
- ci_dfop_sfo_s_s <- summary(saem_biphasic_s)$confint_back
+ ci_dfop_sfo_s_s <- summary(saem_dfop_sfo_s)$confint_back
expect_true(all(ci_dfop_sfo_s_s[, "lower"] < dfop_sfo_pop))
expect_true(all(ci_dfop_sfo_s_s[, "upper"] > dfop_sfo_pop))
# k2 is not fitted well
- ci_dfop_sfo_s_m <- summary(saem_biphasic_m)$confint_back
+ ci_dfop_sfo_s_m <- summary(saem_dfop_sfo_m)$confint_back
expect_true(all(ci_dfop_sfo_s_m[no_k2, "lower"] < dfop_sfo_pop[no_k2]))
expect_true(all(ci_dfop_sfo_s_m[no_k1, "upper"] > dfop_sfo_pop[no_k1]))
# I tried to only do few iterations in routine tests as this is so slow
# but then deSolve fails at some point (presumably at the switch between
# the two types of iterations)
- #saem_biphasic_2 <- saem(mmkin_biphasic, solution_type = "deSolve",
+ #saem_dfop_sfo_2 <- saem(mmkin_biphasic, solution_type = "deSolve",
# control = list(nbiter.saemix = c(10, 5), nbiter.burn = 5), quiet = TRUE)
skip("Fitting with saemix takes around 10 minutes when using deSolve")
- saem_biphasic_2 <- saem(mmkin_biphasic, solution_type = "deSolve", quiet = TRUE)
+ saem_dfop_sfo_2 <- saem(mmkin_dfop_sfo, solution_type = "deSolve", quiet = TRUE)
# As with the analytical solution, k1 and k2 are not fitted well
- ci_dfop_sfo_s_d <- summary(saem_biphasic_2)$confint_back
+ ci_dfop_sfo_s_d <- summary(saem_dfop_sfo_2)$confint_back
expect_true(all(ci_dfop_sfo_s_d[no_k2, "lower"] < dfop_sfo_pop[no_k2]))
expect_true(all(ci_dfop_sfo_s_d[no_k1, "upper"] > dfop_sfo_pop[no_k1]))
})
diff --git a/tests/testthat/test_multistart.R b/tests/testthat/test_multistart.R
index 502cee98..3a511e06 100644
--- a/tests/testthat/test_multistart.R
+++ b/tests/testthat/test_multistart.R
@@ -9,26 +9,22 @@ test_that("multistart works for saem.mmkin models", {
best(saem_sfo_s_multi),
test = TRUE
)
- # On winbuilder, sfo_saem_1 gives an AIC of 1310.8, while we get 1311.7
- # locally (using saemix 3.2, which likely makes the difference due to the
- # error parameter patch) on Linux and Windows. The other, well-determined
- # fits both give 1309.7.
- expect_equal(round(anova_sfo, 1)["sfo_saem_1_reduced", "AIC"], 1309.7)
- expect_equal(round(anova_sfo, 1)["best(saem_sfo_s_multi)", "AIC"], 1309.7)
- expect_true(anova_sfo[3, "Pr(>Chisq)"] > 0.2) # Local: 1, CRAN: 0.4
+ expect_equal(round(anova_sfo, 1)["sfo_saem_1_reduced", "AIC"], 1302.2)
+ expect_equal(round(anova_sfo, 1)["best(saem_sfo_s_multi)", "AIC"], 1302.2)
+ expect_true(anova_sfo[3, "Pr(>Chisq)"] > 0.2) # Local: 1, win-builder: 0.4
set.seed(123456)
- saem_biphasic_m_multi <- multistart(saem_biphasic_m, n = 8,
+ saem_dfop_sfo_m_multi <- multistart(saem_dfop_sfo_m, n = 8,
cores = n_cores)
- expect_known_output(print(saem_biphasic_m_multi),
- file = "print_multistart_biphasic.txt")
+ expect_known_output(print(saem_dfop_sfo_m_multi),
+ file = "print_multistart_dfop_sfo.txt")
- anova_biphasic <- anova(saem_biphasic_m,
- best(saem_biphasic_m_multi))
+ anova_dfop_sfo <- anova(saem_dfop_sfo_m,
+ best(saem_dfop_sfo_m_multi))
# With the new starting parameters we do not improve
# with multistart any more
- expect_equal(anova_biphasic[2, "AIC"], anova_biphasic[1, "AIC"],
+ expect_equal(anova_dfop_sfo[2, "AIC"], anova_dfop_sfo[1, "AIC"],
tolerance = 1e-4)
skip_on_travis() # Plots are platform dependent
@@ -37,10 +33,10 @@ test_that("multistart works for saem.mmkin models", {
vdiffr::expect_doppelganger("llhist for sfo fit", llhist_sfo)
vdiffr::expect_doppelganger("parplot for sfo fit", parplot_sfo)
- llhist_biphasic <- function() llhist(saem_biphasic_m_multi)
- parplot_biphasic <- function() parplot(saem_biphasic_m_multi,
+ llhist_dfop_sfo <- function() llhist(saem_dfop_sfo_m_multi)
+ parplot_dfop_sfo <- function() parplot(saem_dfop_sfo_m_multi,
ylim = c(0.5, 2))
- vdiffr::expect_doppelganger("llhist for biphasic saemix fit", llhist_biphasic)
- vdiffr::expect_doppelganger("parplot for biphasic saemix fit", parplot_biphasic)
+ vdiffr::expect_doppelganger("llhist for dfop sfo fit", llhist_dfop_sfo)
+ vdiffr::expect_doppelganger("parplot for dfop sfo fit", parplot_dfop_sfo)
})
diff --git a/tests/testthat/test_plot.R b/tests/testthat/test_plot.R
index 13058c00..01b0c1ee 100644
--- a/tests/testthat/test_plot.R
+++ b/tests/testthat/test_plot.R
@@ -44,24 +44,24 @@ test_that("Plotting mkinfit, mmkin and mixed model objects is reproducible", {
f_uba_dfop_sfo_saem <- saem(f_uba_mmkin["DFOP-SFO", ], quiet = TRUE, transformations = "saemix")
- plot_biphasic_mmkin <- function() plot(f_uba_dfop_sfo_mixed, pop_curve = TRUE)
- vdiffr::expect_doppelganger("mixed model fit for mmkin object", plot_biphasic_mmkin)
+ plot_dfop_sfo_mmkin <- function() plot(f_uba_dfop_sfo_mixed, pop_curve = TRUE)
+ vdiffr::expect_doppelganger("mixed model fit for mmkin object", plot_dfop_sfo_mmkin)
- plot_biphasic_saem_s <- function() plot(f_uba_dfop_sfo_saem)
- vdiffr::expect_doppelganger("mixed model fit for saem object with saemix transformations", plot_biphasic_saem_s)
+ plot_dfop_sfo_saem_s <- function() plot(f_uba_dfop_sfo_saem)
+ vdiffr::expect_doppelganger("mixed model fit for saem object with saemix transformations", plot_dfop_sfo_saem_s)
skip_on_travis()
- plot_biphasic_nlme <- function() plot(dfop_nlme_1)
- vdiffr::expect_doppelganger("mixed model fit for nlme object", plot_biphasic_nlme)
+ plot_dfop_sfo_nlme <- function() plot(dfop_nlme_1)
+ vdiffr::expect_doppelganger("mixed model fit for nlme object", plot_dfop_sfo_nlme)
- #plot_biphasic_mmkin <- function() plot(mixed(mmkin_biphasic))
+ #plot_dfop_sfo_mmkin <- function() plot(mixed(mmkin_dfop_sfo))
# Biphasic fits with lots of data and fits have lots of potential for differences
- plot_biphasic_nlme <- function() plot(nlme_biphasic)
- #plot_biphasic_saem_s <- function() plot(saem_biphasic_s)
- plot_biphasic_saem_m <- function() plot(saem_biphasic_m)
+ plot_dfop_sfo_nlme <- function() plot(nlme_dfop_sfo)
+ #plot_dfop_sfo_saem_s <- function() plot(saem_dfop_sfo_s)
+ plot_dfop_sfo_saem_m <- function() plot(saem_dfop_sfo_m)
- vdiffr::expect_doppelganger("mixed model fit for saem object with mkin transformations", plot_biphasic_saem_m)
+ vdiffr::expect_doppelganger("mixed model fit for saem object with mkin transformations", plot_dfop_sfo_saem_m)
# different results when working with eigenvalues
plot_errmod_fit_D_obs_eigen <- function() plot_err(fit_D_obs_eigen, sep_obs = FALSE)
diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R
index 20889c6c..31605931 100644
--- a/tests/testthat/test_saemix_parent.R
+++ b/tests/testthat/test_saemix_parent.R
@@ -38,11 +38,11 @@ test_that("Parent fits using saemix are correctly implemented", {
s_sfo_nlme_1 <- summary(sfo_nlme_1)
# Compare with input
- 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)
+ expect_equal(round(s_sfo_saem_1$confint_ranef["SD.k_parent", "est."], 1), 0.3, tol = 0.1)
+ expect_equal(round(s_sfo_saem_1_mkin$confint_ranef["SD.log_k_parent", "est."], 1), 0.3, tol = 0.1)
# k_parent is a bit different from input 0.03 here
- 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)
+ expect_equal(round(s_sfo_saem_1$confint_back["k_parent", "est."], 3), 0.033)
+ expect_equal(round(s_sfo_saem_1_mkin$confint_back["k_parent", "est."], 3), 0.033)
# But the result is pretty unanimous between methods
expect_equal(round(s_sfo_saem_1_reduced$confint_back["k_parent", "est."], 3),
@@ -74,7 +74,7 @@ test_that("Parent fits using saemix are correctly implemented", {
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", no_random_effect = "parent_0")
- fomc_pop <- as.numeric(fomc_pop)
+ fomc_pop <- as.numeric(attr(ds_fomc, "pop"))
ci_fomc_s1 <- summary(fomc_saem_1)$confint_back
expect_true(all(ci_fomc_s1[, "lower"] < fomc_pop))
expect_true(all(ci_fomc_s1[, "upper"] > fomc_pop))
@@ -87,14 +87,14 @@ test_that("Parent fits using saemix are correctly implemented", {
expect_equal(endpoints(fomc_saem_1), endpoints(fomc_saem_2), tol = 0.01)
# DFOP
- dfop_saemix_2 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "saemix",
+ dfop_saem_2 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "saemix",
no_random_effect = "parent_0")
- s_dfop_s1 <- summary(dfop_saemix_1) # mkin transformations
- s_dfop_s2 <- summary(dfop_saemix_2) # saemix transformations
+ s_dfop_s1 <- summary(dfop_saem_1) # mkin transformations
+ s_dfop_s2 <- summary(dfop_saem_2) # saemix transformations
s_dfop_n <- summary(dfop_nlme_1)
- dfop_pop <- as.numeric(dfop_pop)
+ dfop_pop <- as.numeric(attr(ds_dfop, "pop"))
expect_true(all(s_dfop_s1$confint_back[, "lower"] < dfop_pop))
expect_true(all(s_dfop_s1$confint_back[, "upper"] > dfop_pop))
@@ -111,18 +111,18 @@ test_that("Parent fits using saemix are correctly implemented", {
# SFORB
mmkin_sforb_1 <- mmkin("SFORB", ds_dfop, quiet = TRUE, cores = n_cores)
- sforb_saemix_1 <- saem(mmkin_sforb_1, quiet = TRUE,
+ sforb_saem_1 <- saem(mmkin_sforb_1, quiet = TRUE,
no_random_effect = c("parent_free_0"),
transformations = "mkin")
- sforb_saemix_2 <- saem(mmkin_sforb_1, quiet = TRUE,
+ sforb_saem_2 <- saem(mmkin_sforb_1, quiet = TRUE,
no_random_effect = c("parent_free_0"),
transformations = "saemix")
expect_equal(
- log(endpoints(dfop_saemix_1)$distimes[1:2]),
- log(endpoints(sforb_saemix_1)$distimes[1:2]), tolerance = 0.01)
+ log(endpoints(dfop_saem_1)$distimes[1:2]),
+ log(endpoints(sforb_saem_1)$distimes[1:2]), tolerance = 0.01)
expect_equal(
- log(endpoints(sforb_saemix_1)$distimes[1:2]),
- log(endpoints(sforb_saemix_2)$distimes[1:2]), tolerance = 0.01)
+ log(endpoints(sforb_saem_1)$distimes[1:2]),
+ log(endpoints(sforb_saem_2)$distimes[1:2]), tolerance = 0.01)
mmkin_hs_1 <- mmkin("HS", ds_hs, quiet = TRUE, error_model = "const", cores = n_cores)
hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE, no_random_effect = "parent_0")
@@ -131,7 +131,7 @@ test_that("Parent fits using saemix are correctly implemented", {
expect_equal(endpoints(hs_saem_1), endpoints(hs_saem_2), tol = 0.01)
ci_hs_s1 <- summary(hs_saem_1)$confint_back
- hs_pop <- as.numeric(hs_pop)
+ hs_pop <- as.numeric(attr(ds_hs, "pop"))
#expect_true(all(ci_hs_s1[, "lower"] < hs_pop)) # k1 is overestimated
expect_true(all(ci_hs_s1[, "upper"] > hs_pop))
})
@@ -141,10 +141,10 @@ test_that("We can also use mkin solution methods for saem", {
"saemix transformations is only supported if an analytical solution is implemented"
)
skip("This still takes almost 2.5 minutes although we do not solve ODEs")
- dfop_saemix_3 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "mkin",
- solution_type = "analytical", no_random_effect = "parent_0")
- distimes_dfop <- endpoints(dfop_saemix_1)$distimes
- distimes_dfop_analytical <- endpoints(dfop_saemix_3)$distimes
+ dfop_saem_3 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "mkin",
+ solution_type = "analytical", no_random_effect = c("parent_0", "g_qlogis"))
+ distimes_dfop <- endpoints(dfop_saem_1)$distimes
+ distimes_dfop_analytical <- endpoints(dfop_saem_3)$distimes
rel_diff <- abs(distimes_dfop_analytical - distimes_dfop) / distimes_dfop
expect_true(all(rel_diff < 0.01))
})
--
cgit v1.2.3
From 61b9a4046582da5cf88bd3c04d0d6ca77adc3405 Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Fri, 25 Nov 2022 19:10:28 +0100
Subject: mhmkin: Easy specification of ill-defined parms
The argument 'no_random_effect' now accepts an illparms.mhmkin object
---
NEWS.md | 2 +-
R/mhmkin.R | 38 ++++++++++++++++++++------------------
man/mhmkin.Rd | 16 +++++++---------
3 files changed, 28 insertions(+), 28 deletions(-)
(limited to 'R')
diff --git a/NEWS.md b/NEWS.md
index 7f14e06b..f07c6983 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,6 @@
# mkin 1.2.2
--
+- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models.
# mkin 1.2.1 (2022-11-19)
diff --git a/R/mhmkin.R b/R/mhmkin.R
index 1f29dc40..6a61e8bb 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -12,13 +12,13 @@
#' degradation models to the same data
#' @param backend The backend to be used for fitting. Currently, only saemix is
#' supported
-#' @param no_random_effect Default is NULL and will be passed to [saem]. If
-#' you specify "auto", random effects are only included if the number
-#' of datasets in which the parameter passed the t-test is at least 'auto_ranef_threshold'.
-#' Beware that while this may make for convenient model reduction or even
-#' numerical stability of the algorithm, it will likely lead to
-#' underparameterised models.
-#' @param auto_ranef_threshold See 'no_random_effect.
+#' @param no_random_effect Default is NULL and will be passed to [saem]. If a
+#' character vector is supplied, it will be passed to all calls to [saem],
+#' regardless if the corresponding parameter is in the model. Alternatively,
+#' an object of class [illparms.mhmkin] can be specified. This has to have
+#' the same dimensions as the return object of the current call. In this way,
+#' ill-defined parameters found in corresponding simpler versions of the
+#' degradation model can be specified.
#' @param algorithm The algorithm to be used for fitting (currently not used)
#' @param \dots Further arguments that will be passed to the nonlinear mixed-effects
#' model fitting function.
@@ -51,7 +51,7 @@ mhmkin.mmkin <- function(objects, ...) {
#' @export
#' @rdname mhmkin
mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
- no_random_effect = NULL, auto_ranef_threshold = 3,
+ no_random_effect = NULL,
...,
cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), cluster = NULL)
{
@@ -102,22 +102,24 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
deg_index <- w[1]
error_index <- w[2]
mmkin_row <- objects[[error_index]][deg_index, ]
- if (identical(no_random_effect, "auto")) {
- ip <- illparms(mmkin_row)
- ipt <- table(unlist(lapply(as.list(ip), strsplit, ", ")))
- no_ranef <- names(ipt)[which((length(ds) - ipt) <= auto_ranef_threshold)]
- transparms <- rownames(mmkin_row[[1]]$start_transformed)
- bparms <- rownames(mmkin_row[[1]]$start)
- names(transparms) <- bparms
- no_ranef_trans <- transparms[no_ranef]
+ if (is(no_random_effect, "illparms.mhmkin")) {
+ if (identical(dim(no_random_effect), dim(fit_indices))) {
+ no_ranef_split <- strsplit(no_random_effect[[fit_index]], ", ")
+ no_ranef <- sapply(no_ranef_split, function(x) {
+ gsub("sd\\((.*)\\)", "\\1", x)
+ })
+ } else {
+ stop("Dimensions of illparms.mhmkin object given as 'no_random_effect' are not suitable")
+ }
} else {
- no_ranef_trans <- no_random_effect
+ no_ranef <- no_random_effect
}
res <- try(do.call(backend_function,
- args = c(list(mmkin_row), dot_args, list(no_random_effect = no_ranef_trans))))
+ args = c(list(mmkin_row), dot_args, list(no_random_effect = no_ranef))))
return(res)
}
+
fit_time <- system.time({
if (is.null(cluster)) {
results <- parallel::mclapply(as.list(1:n.fits), fit_function,
diff --git a/man/mhmkin.Rd b/man/mhmkin.Rd
index 597fa8ac..4230e44f 100644
--- a/man/mhmkin.Rd
+++ b/man/mhmkin.Rd
@@ -18,7 +18,6 @@ mhmkin(objects, ...)
backend = "saemix",
algorithm = "saem",
no_random_effect = NULL,
- auto_ranef_threshold = 3,
...,
cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(),
cluster = NULL
@@ -42,14 +41,13 @@ supported}
\item{algorithm}{The algorithm to be used for fitting (currently not used)}
-\item{no_random_effect}{Default is NULL and will be passed to \link{saem}. If
-you specify "auto", random effects are only included if the number
-of datasets in which the parameter passed the t-test is at least 'auto_ranef_threshold'.
-Beware that while this may make for convenient model reduction or even
-numerical stability of the algorithm, it will likely lead to
-underparameterised models.}
-
-\item{auto_ranef_threshold}{See 'no_random_effect.}
+\item{no_random_effect}{Default is NULL and will be passed to \link{saem}. If a
+character vector is supplied, it will be passed to all calls to \link{saem},
+regardless if the corresponding parameter is in the model. Alternatively,
+an object of class \link{illparms.mhmkin} can be specified. This has to have
+the same dimensions as the return object of the current call. In this way,
+ill-defined parameters found in corresponding simpler versions of the
+degradation model can be specified.}
\item{cores}{The number of cores to be used for multicore processing. This
is only used when the \code{cluster} argument is \code{NULL}. On Windows
--
cgit v1.2.3
From aaa4cab7e0c7212f91147a9789af54b97fe342ca Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Tue, 29 Nov 2022 20:23:17 +0100
Subject: Complete starting values in summary for saem.mmkin fits
Also update tests to the changes in mhmkin (see NEWS)
---
NEWS.md | 4 ++-
R/summary.saem.mmkin.R | 11 ++++++-
log/build.log | 2 +-
log/test.log | 34 +++++++++++-----------
.../illparms_hfits_synth_no_ranef_auto.txt | 4 ---
tests/testthat/print_hfits_synth_no_ranef_auto.txt | 9 ------
tests/testthat/summary_hfit_sfo_tc.txt | 11 ++++++-
tests/testthat/summary_saem_dfop_sfo_s.txt | 15 +++++++++-
tests/testthat/test_mhmkin.R | 10 -------
9 files changed, 55 insertions(+), 45 deletions(-)
delete mode 100644 tests/testthat/illparms_hfits_synth_no_ranef_auto.txt
delete mode 100644 tests/testthat/print_hfits_synth_no_ranef_auto.txt
(limited to 'R')
diff --git a/NEWS.md b/NEWS.md
index f07c6983..35f499d5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,8 @@
# mkin 1.2.2
-- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models.
+- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well.
+
+- 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters
# mkin 1.2.1 (2022-11-19)
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
index 2754e9f0..7337b0f3 100644
--- a/R/summary.saem.mmkin.R
+++ b/R/summary.saem.mmkin.R
@@ -225,13 +225,22 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3)
obs = "Variance unique to each observed variable",
tc = "Two-component variance function"), "\n")
- cat("\nMean of starting values for individual parameters:\n")
+ cat("\nStarting values for degradation parameters:\n")
print(x$mean_dp_start, digits = digits)
cat("\nFixed degradation parameter values:\n")
if(length(x$fixed$value) == 0) cat("None\n")
else print(x$fixed, digits = digits)
+ cat("\nStarting values for random effects (square root of initial entries in omega):\n")
+ print(sqrt(x$so@model@omega.init), digits = digits)
+
+ cat("\nStarting values for error model parameters:\n")
+ errparms <- x$so@model@error.init
+ names(errparms) <- x$so@model@name.sigma
+ errparms <- errparms[x$so@model@indx.res]
+ print(errparms, digits = digits)
+
cat("\nResults:\n\n")
cat("Likelihood computed by importance sampling\n")
print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
diff --git a/log/build.log b/log/build.log
index 6be01938..dbe0cd5b 100644
--- a/log/build.log
+++ b/log/build.log
@@ -5,5 +5,5 @@
* creating vignettes ... OK
* checking for LF line-endings in source and make files and shell scripts
* checking for empty or unneeded directories
-* building ‘mkin_1.2.1.tar.gz’
+* building ‘mkin_1.2.2.tar.gz’
diff --git a/log/test.log b/log/test.log
index 0b9aa7eb..e17ecc1f 100644
--- a/log/test.log
+++ b/log/test.log
@@ -5,54 +5,54 @@
✔ | 5 | Calculation of Akaike weights
✔ | 3 | Export dataset for reading into CAKE
✔ | 12 | Confidence intervals and p-values [1.1s]
-✔ | 1 12 | Dimethenamid data from 2018 [33.3s]
+✔ | 1 12 | Dimethenamid data from 2018 [32.2s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_dmta.R:98'): Different backends get consistent results for SFO-SFO3+, dimethenamid data
Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
────────────────────────────────────────────────────────────────────────────────
-✔ | 14 | Error model fitting [5.0s]
+✔ | 14 | Error model fitting [4.9s]
✔ | 5 | Time step normalisation
✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s]
✔ | 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.4s]
-✔ | 1 | Fitting the logistic model [0.3s]
-✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [25.5s]
-✔ | 1 11 | Nonlinear mixed-effects models [13.2s]
+✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3s]
+✔ | 1 | Fitting the logistic model [0.2s]
+✔ | 8 | Batch fitting and diagnosing hierarchical kinetic models [14.5s]
+✔ | 1 11 | Nonlinear mixed-effects models [13.1s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_mixed.R:78'): saemix results are reproducible for biphasic fits
Reason: Fitting with saemix takes around 10 minutes when using deSolve
────────────────────────────────────────────────────────────────────────────────
✔ | 3 | Test dataset classes mkinds and mkindsg
-✔ | 10 | Special cases of mkinfit calls [0.5s]
+✔ | 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.3s]
-✔ | 12 | Multistart method for saem.mmkin models [51.1s]
-✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.3s]
-✔ | 9 | Nonlinear mixed-effects models with nlme [9.2s]
-✔ | 15 | Plotting [10.3s]
+✔ | 12 | Multistart method for saem.mmkin models [50.1s]
+✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.2s]
+✔ | 9 | Nonlinear mixed-effects models with nlme [8.7s]
+✔ | 15 | Plotting [10.2s]
✔ | 4 | Residuals extracted from mkinfit models
-✔ | 1 36 | saemix parent models [73.7s]
+✔ | 1 36 | saemix parent models [103.8s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_saemix_parent.R:143'): We can also use mkin solution methods for saem
Reason: This still takes almost 2.5 minutes although we do not solve ODEs
────────────────────────────────────────────────────────────────────────────────
-✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5s]
+✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
✔ | 11 | Processing of residue series
✔ | 10 | 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.4s]
-✔ | 9 | Hypothesis tests [8.8s]
+✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
+✔ | 9 | Hypothesis tests [8.1s]
✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 251.2 s
+Duration: 266.0 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)
• This still takes almost 2.5 minutes although we do not solve ODEs (1)
-[ FAIL 0 | WARN 0 | SKIP 3 | PASS 270 ]
+[ FAIL 0 | WARN 0 | SKIP 3 | PASS 268 ]
Error while shutting down parallel: unable to terminate some child processes
diff --git a/tests/testthat/illparms_hfits_synth_no_ranef_auto.txt b/tests/testthat/illparms_hfits_synth_no_ranef_auto.txt
deleted file mode 100644
index a64ed222..00000000
--- a/tests/testthat/illparms_hfits_synth_no_ranef_auto.txt
+++ /dev/null
@@ -1,4 +0,0 @@
- error
-degradation const tc
- SFO sd(parent_0)
- FOMC b.1
diff --git a/tests/testthat/print_hfits_synth_no_ranef_auto.txt b/tests/testthat/print_hfits_synth_no_ranef_auto.txt
deleted file mode 100644
index 9af1cbcd..00000000
--- a/tests/testthat/print_hfits_synth_no_ranef_auto.txt
+++ /dev/null
@@ -1,9 +0,0 @@
- object
-Status of individual fits:
-
- error
-degradation const tc
- SFO OK OK
- FOMC OK OK
-
-OK: Fit terminated successfully
diff --git a/tests/testthat/summary_hfit_sfo_tc.txt b/tests/testthat/summary_hfit_sfo_tc.txt
index bb5bf6fb..0a61f75f 100644
--- a/tests/testthat/summary_hfit_sfo_tc.txt
+++ b/tests/testthat/summary_hfit_sfo_tc.txt
@@ -17,13 +17,22 @@ Using 300, 100 iterations and 9 chains
Variance model: Two-component variance function
-Mean of starting values for individual parameters:
+Starting values for degradation parameters:
parent_0 log_k_parent
101 -3
Fixed degradation parameter values:
None
+Starting values for random effects (square root of initial entries in omega):
+ parent_0 log_k_parent
+parent_0 4 0.0
+log_k_parent 0 0.4
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
Results:
Likelihood computed by importance sampling
diff --git a/tests/testthat/summary_saem_dfop_sfo_s.txt b/tests/testthat/summary_saem_dfop_sfo_s.txt
index 7c337843..6468ff17 100644
--- a/tests/testthat/summary_saem_dfop_sfo_s.txt
+++ b/tests/testthat/summary_saem_dfop_sfo_s.txt
@@ -22,7 +22,7 @@ Using 300, 100 iterations and 4 chains
Variance model: Two-component variance function
-Mean of starting values for individual parameters:
+Starting values for degradation parameters:
parent_0 k_m1 f_parent_to_m1 k1 k2
1e+02 7e-03 5e-01 1e-01 2e-02
g
@@ -31,6 +31,19 @@ Mean of starting values for individual parameters:
Fixed degradation parameter values:
None
+Starting values for random effects (square root of initial entries in omega):
+ parent_0 k_m1 f_parent_to_m1 k1 k2 g
+parent_0 101 0 0 0 0 0
+k_m1 0 1 0 0 0 0
+f_parent_to_m1 0 0 1 0 0 0
+k1 0 0 0 1 0 0
+k2 0 0 0 0 1 0
+g 0 0 0 0 0 1
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
Results:
Likelihood computed by importance sampling
diff --git a/tests/testthat/test_mhmkin.R b/tests/testthat/test_mhmkin.R
index 93333ac1..e2339f28 100644
--- a/tests/testthat/test_mhmkin.R
+++ b/tests/testthat/test_mhmkin.R
@@ -46,14 +46,4 @@ test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", {
print(fits_synth_const),
"print_fits_synth_const.txt")
- hfits_no_ranef_auto <- update(hfits, no_random_effect = "auto", auto_ranef_threshold = 2)
-
- expect_known_output(
- print(hfits_no_ranef_auto),
- "print_hfits_synth_no_ranef_auto.txt")
-
- expect_known_output(
- print(illparms(hfits_no_ranef_auto)),
- "illparms_hfits_synth_no_ranef_auto.txt")
-
})
--
cgit v1.2.3
From 74e44dfed5af6e6fd421abe82d3e3f190771f85a Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Thu, 1 Dec 2022 11:20:00 +0100
Subject: Possibility to manually specify no_random_effects in mhmkin
---
NEWS.md | 2 +-
R/mhmkin.R | 91 ++++++++++++++++++++++++-------
log/check.log | 8 +--
log/test.log | 48 ++++++++--------
man/mhmkin.Rd | 50 +++++++++++++++--
tests/testthat/illparms_hfits_synth.txt | 10 +---
tests/testthat/print_fits_synth_const.txt | 4 +-
tests/testthat/summary_hfit_sfo_tc.txt | 34 ++++++------
tests/testthat/test_mhmkin.R | 34 ++++++++----
9 files changed, 187 insertions(+), 94 deletions(-)
(limited to 'R')
diff --git a/NEWS.md b/NEWS.md
index 35f499d5..65df4191 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,6 @@
# mkin 1.2.2
-- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well.
+- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object or a list with suitable dimensions as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well.
- 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters
diff --git a/R/mhmkin.R b/R/mhmkin.R
index 6a61e8bb..a9798ff4 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -14,11 +14,12 @@
#' supported
#' @param no_random_effect Default is NULL and will be passed to [saem]. If a
#' character vector is supplied, it will be passed to all calls to [saem],
-#' regardless if the corresponding parameter is in the model. Alternatively,
-#' an object of class [illparms.mhmkin] can be specified. This has to have
-#' the same dimensions as the return object of the current call. In this way,
-#' ill-defined parameters found in corresponding simpler versions of the
-#' degradation model can be specified.
+#' which will exclude random effects for all matching parameters. Alternatively,
+#' a list of character vectors or an object of class [illparms.mhmkin] can be
+#' specified. They have to have the same dimensions that the return object of
+#' the current call will have, i.e. the number of rows must match the number
+#' of degradation models in the mmkin object(s), and the number of columns must
+#' match the number of error models used in the mmkin object(s).
#' @param algorithm The algorithm to be used for fitting (currently not used)
#' @param \dots Further arguments that will be passed to the nonlinear mixed-effects
#' model fitting function.
@@ -50,6 +51,42 @@ mhmkin.mmkin <- function(objects, ...) {
#' @export
#' @rdname mhmkin
+#' @examples
+#' \dontrun{
+#' # We start with separate evaluations of all the first six datasets with two
+#' # degradation models and two error models
+#' f_sep_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = 2, quiet = TRUE)
+#' f_sep_tc <- update(f_sep_const, error_model = "tc")
+#' # The mhmkin function sets up hierarchical degradation models aka
+#' # nonlinear mixed-effects models for all four combinations, specifying
+#' # uncorrelated random effects for all degradation parameters
+#' f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cores = 2)
+#' status(f_saem_1)
+#' # The 'illparms' function shows that in all hierarchical fits, at least
+#' # one random effect is ill-defined (the confidence interval for the
+#' # random effect expressed as standard deviation includes zero)
+#' illparms(f_saem_1)
+#' # Therefore we repeat the fits, excluding the ill-defined random effects
+#' f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1))
+#' status(f_saem_2)
+#' illparms(f_saem_2)
+#' # Model comparisons show that FOMC with two-component error is preferable,
+#' # and confirms our reduction of the default parameter model
+#' anova(f_saem_1)
+#' anova(f_saem_2)
+#' # The convergence plot for the selected model looks fine
+#' saemix::plot(f_saem_2[["FOMC", "tc"]]$so, plot.type = "convergence")
+#' # The plot of predictions versus data shows that we have a pretty data-rich
+#' # situation with homogeneous distribution of residuals, because we used the
+#' # same degradation model, error model and parameter distribution model that
+#' # was used in the data generation.
+#' plot(f_saem_2[["FOMC", "tc"]])
+#' # We can specify the same parameter model reductions manually
+#' no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta"))
+#' dim(no_ranef) <- c(2, 2)
+#' f_saem_2m <- update(f_saem_1, no_random_effect = no_ranef)
+#' anova(f_saem_2m)
+#' }
mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
no_random_effect = NULL,
...,
@@ -97,25 +134,38 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
dimnames(fit_indices) <- list(degradation = names(deg_models),
error = error_models)
- fit_function <- function(fit_index) {
- w <- which(fit_indices == fit_index, arr.ind = TRUE)
- deg_index <- w[1]
- error_index <- w[2]
- mmkin_row <- objects[[error_index]][deg_index, ]
+ if (is.null(no_random_effect) || length(dim(no_random_effect)) == 1) {
+ no_ranef <- rep(list(no_random_effect), n.fits)
+ dim(no_ranef) <- dim(fit_indices)
+ } else {
+ if (!identical(dim(no_random_effect), dim(fit_indices))) {
+ stop("Dimensions of argument 'no_random_effect' are not suitable")
+ }
if (is(no_random_effect, "illparms.mhmkin")) {
- if (identical(dim(no_random_effect), dim(fit_indices))) {
- no_ranef_split <- strsplit(no_random_effect[[fit_index]], ", ")
- no_ranef <- sapply(no_ranef_split, function(x) {
- gsub("sd\\((.*)\\)", "\\1", x)
+ no_ranef_dim <- dim(no_random_effect)
+ no_ranef <- lapply(no_random_effect, function(x) {
+ no_ranef_split <- strsplit(x, ", ")
+ ret <- sapply(no_ranef_split, function(y) {
+ gsub("sd\\((.*)\\)", "\\1", y)
})
- } else {
- stop("Dimensions of illparms.mhmkin object given as 'no_random_effect' are not suitable")
- }
+ return(ret)
+ })
+ dim(no_ranef) <- no_ranef_dim
} else {
no_ranef <- no_random_effect
}
+ }
+
+ fit_function <- function(fit_index) {
+ w <- which(fit_indices == fit_index, arr.ind = TRUE)
+ deg_index <- w[1]
+ error_index <- w[2]
+ mmkin_row <- objects[[error_index]][deg_index, ]
res <- try(do.call(backend_function,
- args = c(list(mmkin_row), dot_args, list(no_random_effect = no_ranef))))
+ args = c(
+ list(mmkin_row),
+ dot_args,
+ list(no_random_effect = no_ranef[[deg_index, error_index]]))))
return(res)
}
@@ -145,15 +195,16 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
#' @param j Column index selecting the fits to specific datasets
#' @param drop If FALSE, the method always returns an mhmkin object, otherwise
#' either a list of fit objects or a single fit object.
-#' @return An object of class \code{\link{mhmkin}}.
+#' @return An object inheriting from \code{\link{mhmkin}}.
#' @rdname mhmkin
#' @export
`[.mhmkin` <- function(x, i, j, ..., drop = FALSE) {
+ original_class <- class(x)
class(x) <- NULL
x_sub <- x[i, j, drop = drop]
if (!drop) {
- class(x_sub) <- "mhmkin"
+ class(x_sub) <- original_class
}
return(x_sub)
}
diff --git a/log/check.log b/log/check.log
index 31fc31eb..42365918 100644
--- a/log/check.log
+++ b/log/check.log
@@ -5,7 +5,7 @@
* using options ‘--no-tests --as-cran’
* checking for file ‘mkin/DESCRIPTION’ ... OK
* checking extension type ... Package
-* this is package ‘mkin’ version ‘1.2.1’
+* this is package ‘mkin’ version ‘1.2.2’
* package encoding: UTF-8
* checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
Maintainer: ‘Johannes Ranke ’
@@ -18,7 +18,7 @@ Maintainer: ‘Johannes Ranke ’
* checking for portable file names ... OK
* checking for sufficient/correct file permissions ... OK
* checking serialization versions ... OK
-* checking whether package ‘mkin’ can be installed ... OK
+* checking whether package ‘mkin’ can be installed ... [11s/11s] OK
* checking installed package size ... OK
* checking package directory ... OK
* checking for future file timestamps ... OK
@@ -41,7 +41,7 @@ Maintainer: ‘Johannes Ranke ’
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
-* checking R code for possible problems ... [17s/17s] OK
+* checking R code for possible problems ... [19s/19s] OK
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd line widths ... OK
@@ -57,7 +57,7 @@ Maintainer: ‘Johannes Ranke ’
* checking data for ASCII and uncompressed saves ... OK
* checking installed files from ‘inst/doc’ ... OK
* checking files in ‘vignettes’ ... OK
-* checking examples ... [20s/20s] OK
+* checking examples ... [24s/24s] OK
* checking for unstated dependencies in ‘tests’ ... OK
* checking tests ... SKIPPED
* checking for unstated dependencies in vignettes ... OK
diff --git a/log/test.log b/log/test.log
index e17ecc1f..84fa49b9 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,58 +1,58 @@
ℹ Testing mkin
✔ | F W S OK | Context
✔ | 5 | AIC calculation
-✔ | 5 | Analytical solutions for coupled models [3.3s]
+✔ | 5 | Analytical solutions for coupled models [4.2s]
✔ | 5 | Calculation of Akaike weights
✔ | 3 | Export dataset for reading into CAKE
-✔ | 12 | Confidence intervals and p-values [1.1s]
-✔ | 1 12 | Dimethenamid data from 2018 [32.2s]
+✔ | 12 | Confidence intervals and p-values [1.2s]
+✔ | 1 12 | Dimethenamid data from 2018 [42.0s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_dmta.R:98'): Different backends get consistent results for SFO-SFO3+, dimethenamid data
Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
────────────────────────────────────────────────────────────────────────────────
-✔ | 14 | Error model fitting [4.9s]
+✔ | 14 | Error model fitting [6.5s]
✔ | 5 | Time step normalisation
-✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s]
-✔ | 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]
-✔ | 8 | Batch fitting and diagnosing hierarchical kinetic models [14.5s]
-✔ | 1 11 | Nonlinear mixed-effects models [13.1s]
+✔ | 4 | Calculation of FOCUS chi2 error levels [0.7s]
+✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.1s]
+✔ | 4 | Test fitting the decline of metabolites from their maximum [0.5s]
+✔ | 1 | Fitting the logistic model [0.3s]
+✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [54.1s]
+✔ | 1 11 | Nonlinear mixed-effects models [14.3s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_mixed.R:78'): saemix results are reproducible for biphasic fits
Reason: Fitting with saemix takes around 10 minutes when using deSolve
────────────────────────────────────────────────────────────────────────────────
✔ | 3 | Test dataset classes mkinds and mkindsg
-✔ | 10 | Special cases of mkinfit calls [0.4s]
-✔ | 3 | mkinfit features [0.7s]
-✔ | 8 | mkinmod model generation and printing [0.2s]
+✔ | 10 | Special cases of mkinfit calls [0.8s]
+✔ | 3 | mkinfit features [0.9s]
+✔ | 8 | mkinmod model generation and printing [0.3s]
✔ | 3 | Model predictions with mkinpredict [0.3s]
-✔ | 12 | Multistart method for saem.mmkin models [50.1s]
-✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.2s]
-✔ | 9 | Nonlinear mixed-effects models with nlme [8.7s]
-✔ | 15 | Plotting [10.2s]
+✔ | 12 | Multistart method for saem.mmkin models [80.1s]
+✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.8s]
+✔ | 9 | Nonlinear mixed-effects models with nlme [11.4s]
+✔ | 15 | Plotting [12.1s]
✔ | 4 | Residuals extracted from mkinfit models
-✔ | 1 36 | saemix parent models [103.8s]
+✔ | 1 36 | saemix parent models [85.9s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_saemix_parent.R:143'): We can also use mkin solution methods for saem
Reason: This still takes almost 2.5 minutes although we do not solve ODEs
────────────────────────────────────────────────────────────────────────────────
-✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
+✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.9s]
✔ | 11 | Processing of residue series
-✔ | 10 | Fitting the SFORB model [3.8s]
+✔ | 10 | Fitting the SFORB model [4.6s]
✔ | 1 | Summaries of old mkinfit objects
✔ | 5 | Summary [0.2s]
-✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
-✔ | 9 | Hypothesis tests [8.1s]
+✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.9s]
+✔ | 9 | Hypothesis tests [11.0s]
✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 266.0 s
+Duration: 342.6 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)
• This still takes almost 2.5 minutes although we do not solve ODEs (1)
-[ FAIL 0 | WARN 0 | SKIP 3 | PASS 268 ]
+[ FAIL 0 | WARN 0 | SKIP 3 | PASS 270 ]
Error while shutting down parallel: unable to terminate some child processes
diff --git a/man/mhmkin.Rd b/man/mhmkin.Rd
index 4230e44f..c77f4289 100644
--- a/man/mhmkin.Rd
+++ b/man/mhmkin.Rd
@@ -43,11 +43,12 @@ supported}
\item{no_random_effect}{Default is NULL and will be passed to \link{saem}. If a
character vector is supplied, it will be passed to all calls to \link{saem},
-regardless if the corresponding parameter is in the model. Alternatively,
-an object of class \link{illparms.mhmkin} can be specified. This has to have
-the same dimensions as the return object of the current call. In this way,
-ill-defined parameters found in corresponding simpler versions of the
-degradation model can be specified.}
+which will exclude random effects for all matching parameters. Alternatively,
+a list of character vectors or an object of class \link{illparms.mhmkin} can be
+specified. They have to have the same dimensions that the return object of
+the current call will have, i.e. the number of rows must match the number
+of degradation models in the mmkin object(s), and the number of columns must
+match the number of error models used in the mmkin object(s).}
\item{cores}{The number of cores to be used for multicore processing. This
is only used when the \code{cluster} argument is \code{NULL}. On Windows
@@ -74,7 +75,7 @@ be indexed using the degradation model names for the first index (row index)
and the error model names for the second index (column index), with class
attribute 'mhmkin'.
-An object of class \code{\link{mhmkin}}.
+An object inheriting from \code{\link{mhmkin}}.
}
\description{
The name of the methods expresses that (\strong{m}ultiple) \strong{h}ierarchichal
@@ -82,6 +83,43 @@ The name of the methods expresses that (\strong{m}ultiple) \strong{h}ierarchicha
fitted. Our kinetic models are nonlinear, so we can use various nonlinear
mixed-effects model fitting functions.
}
+\examples{
+\dontrun{
+# We start with separate evaluations of all the first six datasets with two
+# degradation models and two error models
+f_sep_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = 2, quiet = TRUE)
+f_sep_tc <- update(f_sep_const, error_model = "tc")
+# The mhmkin function sets up hierarchical degradation models aka
+# nonlinear mixed-effects models for all four combinations, specifying
+# uncorrelated random effects for all degradation parameters
+f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cores = 2)
+status(f_saem_1)
+# The 'illparms' function shows that in all hierarchical fits, at least
+# one random effect is ill-defined (the confidence interval for the
+# random effect expressed as standard deviation includes zero)
+illparms(f_saem_1)
+# Therefore we repeat the fits, excluding the ill-defined random effects
+f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1))
+status(f_saem_2)
+illparms(f_saem_2)
+# Model comparisons show that FOMC with two-component error is preferable,
+# and confirms our reduction of the default parameter model
+anova(f_saem_1)
+anova(f_saem_2)
+# The convergence plot for the selected model looks fine
+saemix::plot(f_saem_2[["FOMC", "tc"]]$so, plot.type = "convergence")
+# The plot of predictions versus data shows that we have a pretty data-rich
+# situation with homogeneous distribution of residuals, because we used the
+# same degradation model, error model and parameter distribution model that
+# was used in the data generation.
+plot(f_saem_2[["FOMC", "tc"]])
+# We can specify the same parameter model reductions manually
+no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta"))
+dim(no_ranef) <- c(2, 2)
+f_saem_2m <- update(f_saem_1, no_random_effect = no_ranef)
+anova(f_saem_2m)
+}
+}
\seealso{
\code{\link{[.mhmkin}} for subsetting \link{mhmkin} objects
}
diff --git a/tests/testthat/illparms_hfits_synth.txt b/tests/testthat/illparms_hfits_synth.txt
index affd1318..7a69645b 100644
--- a/tests/testthat/illparms_hfits_synth.txt
+++ b/tests/testthat/illparms_hfits_synth.txt
@@ -1,8 +1,4 @@
error
-degradation const
- SFO
- FOMC sd(log_alpha), sd(log_beta)
- error
-degradation tc
- SFO sd(parent_0)
- FOMC sd(parent_0), sd(log_alpha), sd(log_beta)
+degradation const tc
+ SFO sd(parent_0) sd(parent_0)
+ FOMC sd(log_beta) sd(parent_0), sd(log_beta)
diff --git a/tests/testthat/print_fits_synth_const.txt b/tests/testthat/print_fits_synth_const.txt
index b4bbe6ca..5d076d3d 100644
--- a/tests/testthat/print_fits_synth_const.txt
+++ b/tests/testthat/print_fits_synth_const.txt
@@ -4,8 +4,6 @@ Status of individual fits:
dataset
model 1 2 3 4 5 6
SFO OK OK OK OK OK OK
- FOMC C OK OK OK OK C
+ FOMC OK OK OK OK OK OK
-C: Optimisation did not converge:
-false convergence (8)
OK: No warnings
diff --git a/tests/testthat/summary_hfit_sfo_tc.txt b/tests/testthat/summary_hfit_sfo_tc.txt
index 0a61f75f..0618c715 100644
--- a/tests/testthat/summary_hfit_sfo_tc.txt
+++ b/tests/testthat/summary_hfit_sfo_tc.txt
@@ -8,7 +8,7 @@ Equations:
d_parent/dt = - k_parent * parent
Data:
-104 observations of 1 variable(s) grouped in 6 datasets
+95 observations of 1 variable(s) grouped in 6 datasets
Model predictions using solution type analytical
@@ -19,7 +19,7 @@ Variance model: Two-component variance function
Starting values for degradation parameters:
parent_0 log_k_parent
- 101 -3
+ 94 -2
Fixed degradation parameter values:
None
@@ -27,7 +27,7 @@ None
Starting values for random effects (square root of initial entries in omega):
parent_0 log_k_parent
parent_0 4 0.0
-log_k_parent 0 0.4
+log_k_parent 0 0.7
Starting values for error model parameters:
a.1 b.1
@@ -37,15 +37,15 @@ Results:
Likelihood computed by importance sampling
AIC BIC logLik
- 524 523 -257
+ 542 541 -266
Optimised parameters:
- est. lower upper
-parent_0 100.68 99.27 102.08
-log_k_parent -3.38 -3.55 -3.21
-a.1 0.87 0.59 1.14
-b.1 0.05 0.04 0.06
-SD.log_k_parent 0.21 0.09 0.33
+ est. lower upper
+parent_0 92.52 89.11 95.9
+log_k_parent -1.66 -2.07 -1.3
+a.1 2.03 1.60 2.5
+b.1 0.09 0.07 0.1
+SD.log_k_parent 0.51 0.22 0.8
Correlation:
pr_0
@@ -53,18 +53,18 @@ log_k_parent 0.1
Random effects:
est. lower upper
-SD.log_k_parent 0.2 0.09 0.3
+SD.log_k_parent 0.5 0.2 0.8
Variance model:
est. lower upper
-a.1 0.87 0.59 1.14
-b.1 0.05 0.04 0.06
+a.1 2.03 1.60 2.5
+b.1 0.09 0.07 0.1
Backtransformed parameters:
- est. lower upper
-parent_0 1e+02 99.27 1e+02
-k_parent 3e-02 0.03 4e-02
+ est. lower upper
+parent_0 92.5 89.1 95.9
+k_parent 0.2 0.1 0.3
Estimated disappearance times:
DT50 DT90
-parent 20 68
+parent 4 12
diff --git a/tests/testthat/test_mhmkin.R b/tests/testthat/test_mhmkin.R
index e2339f28..da063326 100644
--- a/tests/testthat/test_mhmkin.R
+++ b/tests/testthat/test_mhmkin.R
@@ -3,8 +3,11 @@ context("Batch fitting and diagnosing hierarchical kinetic models")
test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", {
skip_on_cran()
- fits_synth_const <- suppressWarnings(
- mmkin(c("SFO", "FOMC"), ds_sfo[1:6], cores = n_cores, quiet = TRUE))
+ fits_synth_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = n_cores, quiet = TRUE)
+
+ expect_known_output(
+ print(fits_synth_const),
+ "print_fits_synth_const.txt")
fits_synth_tc <- suppressWarnings(
update(fits_synth_const, error_model = "tc"))
@@ -19,8 +22,8 @@ test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", {
print(illparms(hfits)),
"illparms_hfits_synth.txt")
- expect_equal(which.min(AIC(hfits)), 3)
- expect_equal(which.min(BIC(hfits)), 3)
+ expect_equal(which.min(AIC(hfits)), 4)
+ expect_equal(which.min(BIC(hfits)), 4)
hfit_sfo_tc <- update(hfits[["SFO", "tc"]],
covariance.model = diag(c(0, 1)))
@@ -38,12 +41,19 @@ test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", {
expect_known_output(print(test_summary, digits = 1),
"summary_hfit_sfo_tc.txt")
- # It depends on the platform exactly which of the datasets fail to converge
- # with FOMC, because they were generated to be SFO
- skip_on_travis()
-
- expect_known_output(
- print(fits_synth_const),
- "print_fits_synth_const.txt")
-
+ hfits_sfo_reduced <- update(hfits,
+ no_random_effect = illparms(hfits))
+ expect_equal(
+ as.character(illparms(hfits_sfo_reduced)),
+ rep("", 4))
+
+ # We can also manually set up an object specifying random effects to be
+ # excluded. Entries in the inital list have to be by column
+ no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta"))
+ dim(no_ranef) <- c(2, 2)
+
+ hfits_sfo_reduced_2 <- update(hfits,
+ no_random_effect = no_ranef)
+ expect_equivalent(round(anova(hfits_sfo_reduced), 0),
+ round(anova(hfits_sfo_reduced_2), 0))
})
--
cgit v1.2.3
From cf73354f985be17352a4d538c6f9ea69f6be9aa9 Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Fri, 2 Dec 2022 13:28:26 +0100
Subject: Avoid redundant warnings in summaries
---
NEWS.md | 2 +-
R/summary.saem.mmkin.R | 5 +++--
2 files changed, 4 insertions(+), 3 deletions(-)
(limited to 'R')
diff --git a/NEWS.md b/NEWS.md
index 65df4191..15c45cc9 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object or a list with suitable dimensions as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well.
-- 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters
+- 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters. Avoid redundant warnings that occurred in the calculation of correlations of the fixed effects in the case that the Fisher information matrix could not be inverted.
# mkin 1.2.1 (2022-11-19)
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
index 7337b0f3..46ab548b 100644
--- a/R/summary.saem.mmkin.R
+++ b/R/summary.saem.mmkin.R
@@ -136,10 +136,11 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
}
# Correlation of fixed effects (inspired by summary.nlme)
- varFix <- try(vcov(object$so)[1:n_fixed, 1:n_fixed])
- if (inherits(varFix, "try-error")) {
+ cov_so <- try(solve(object$so@results@fim), silent = TRUE)
+ if (inherits(cov_so, "try-error")) {
object$corFixed <- NA
} else {
+ varFix <- cov_so[1:n_fixed, 1:n_fixed]
stdFix <- sqrt(diag(varFix))
object$corFixed <- array(
t(varFix/stdFix)/stdFix,
--
cgit v1.2.3
From 478c6d5eec4c84b22b43adcbdf36888b302ead00 Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Tue, 6 Dec 2022 10:33:24 +0100
Subject: Some parplot improvements
llquant argument, improved legend text, tests
---
NEWS.md | 2 +
R/parplot.R | 16 +-
log/test.log | 41 +++--
man/parplot.Rd | 7 +-
.../_snaps/multistart/parplot-for-dfop-sfo-fit.svg | 169 ++++++++++-----------
.../_snaps/multistart/parplot-for-sfo-fit.svg | 4 +-
tests/testthat/test_multistart.R | 2 +-
vignettes/FOCUS_D.html | 37 ++---
vignettes/FOCUS_L.html | 83 +++++-----
9 files changed, 183 insertions(+), 178 deletions(-)
(limited to 'R')
diff --git a/NEWS.md b/NEWS.md
index 15c45cc9..c7929fbe 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,8 @@
- 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters. Avoid redundant warnings that occurred in the calculation of correlations of the fixed effects in the case that the Fisher information matrix could not be inverted.
+- 'R/parplot.R': Possibility to select the top 'llquant' fraction of the fits for the parameter plots, and improved legend text.
+
# mkin 1.2.1 (2022-11-19)
- '{data,R}/ds_mixed.rda': Include the test data in the package instead of generating it in 'tests/testthat/setup_script.R'. Refactor the generating code to make it consistent and update tests.
diff --git a/R/parplot.R b/R/parplot.R
index 98579779..63306ac2 100644
--- a/R/parplot.R
+++ b/R/parplot.R
@@ -10,7 +10,10 @@
#'
#' @param object The [multistart] object
#' @param llmin The minimum likelihood of objects to be shown
-#' @param scale By default, scale parameters using the best available fit.
+#' @param llquant Fractional value for selecting only the fits with higher
+#' likelihoods. Overrides 'llmin'.
+#' @param scale By default, scale parameters using the best
+#' available fit.
#' If 'median', parameters are scaled using the median parameters from all fits.
#' @param main Title of the plot
#' @param lpos Positioning of the legend.
@@ -28,7 +31,8 @@ parplot <- function(object, ...) {
#' @rdname parplot
#' @export
-parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best", "median"),
+parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
+ scale = c("best", "median"),
lpos = "bottomleft", main = "", ...)
{
oldpar <- par(no.readonly = TRUE)
@@ -48,6 +52,10 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
stop("parplot is only implemented for multistart.saem.mmkin objects")
}
ll <- sapply(object, llfunc)
+ if (!is.na(llquant[1])) {
+ if (llmin != -Inf) warning("Overriding 'llmin' because 'llquant' was specified")
+ llmin <- quantile(ll, 1 - llquant)
+ }
selected <- which(ll > llmin)
selected_parms <- all_parms[selected, ]
@@ -110,7 +118,7 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
legend(lpos, inset = c(0.05, 0.05), bty = "n",
pch = 1, col = 3:1, lty = c(NA, NA, 1),
legend = c(
- "Starting parameters",
- "Original run",
+ "Original start",
+ "Original results",
"Multistart runs"))
}
diff --git a/log/test.log b/log/test.log
index 5764f209..7614b136 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,22 +1,22 @@
ℹ Testing mkin
✔ | F W S OK | Context
✔ | 5 | AIC calculation
-✔ | 5 | Analytical solutions for coupled models [3.3s]
+✔ | 5 | Analytical solutions for coupled models [3.2s]
✔ | 5 | Calculation of Akaike weights
✔ | 3 | Export dataset for reading into CAKE
✔ | 12 | Confidence intervals and p-values [1.1s]
-✔ | 1 12 | Dimethenamid data from 2018 [31.9s]
+✔ | 1 12 | Dimethenamid data from 2018 [32.0s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_dmta.R:98'): Different backends get consistent results for SFO-SFO3+, dimethenamid data
Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
────────────────────────────────────────────────────────────────────────────────
-✔ | 14 | Error model fitting [4.9s]
+✔ | 14 | Error model fitting [4.8s]
✔ | 5 | Time step normalisation
✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s]
✔ | 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]
-✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [41.4s]
+✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [42.2s]
✔ | 1 11 | Nonlinear mixed-effects models [13.3s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_mixed.R:78'): saemix results are reproducible for biphasic fits
@@ -26,32 +26,41 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
✔ | 10 | Special cases of mkinfit calls [0.6s]
✔ | 3 | mkinfit features [0.7s]
✔ | 8 | mkinmod model generation and printing [0.2s]
-✔ | 3 | Model predictions with mkinpredict [0.3s]
-✔ | 12 | Multistart method for saem.mmkin models [47.5s]
-✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.3s]
-✔ | 9 | Nonlinear mixed-effects models with nlme [9.5s]
-✔ | 15 | Plotting [10.2s]
+✔ | 3 | Model predictions with mkinpredict [0.4s]
+✖ | 1 11 | Multistart method for saem.mmkin models [46.7s]
+────────────────────────────────────────────────────────────────────────────────
+Failure ('test_multistart.R:56'): multistart works for saem.mmkin models
+Snapshot of `testcase` to 'multistart/parplot-for-dfop-sfo-fit.svg' has changed
+Run `testthat::snapshot_review('multistart/')` to review changes
+Backtrace:
+ 1. vdiffr::expect_doppelganger("parplot for dfop sfo fit", parplot_dfop_sfo)
+ at test_multistart.R:56:2
+ 3. testthat::expect_snapshot_file(...)
+────────────────────────────────────────────────────────────────────────────────
+✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.2s]
+✔ | 9 | Nonlinear mixed-effects models with nlme [9.4s]
+✔ | 15 | Plotting [10.1s]
✔ | 4 | Residuals extracted from mkinfit models
-✔ | 1 36 | saemix parent models [73.0s]
+✔ | 1 36 | saemix parent models [73.6s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_saemix_parent.R:143'): We can also use mkin solution methods for saem
Reason: This still takes almost 2.5 minutes although we do not solve ODEs
────────────────────────────────────────────────────────────────────────────────
-✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5s]
+✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
✔ | 11 | Processing of residue series
-✔ | 10 | Fitting the SFORB model [3.9s]
+✔ | 10 | Fitting the SFORB model [3.7s]
✔ | 1 | Summaries of old mkinfit objects
✔ | 5 | Summary [0.2s]
-✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
-✔ | 9 | Hypothesis tests [8.4s]
+✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3s]
+✔ | 9 | Hypothesis tests [8.1s]
✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 261.1 s
+Duration: 260.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)
• This still takes almost 2.5 minutes although we do not solve ODEs (1)
-[ FAIL 0 | WARN 0 | SKIP 3 | PASS 270 ]
+[ FAIL 1 | WARN 0 | SKIP 3 | PASS 269 ]
diff --git a/man/parplot.Rd b/man/parplot.Rd
index ac9e02cf..67bf0cc1 100644
--- a/man/parplot.Rd
+++ b/man/parplot.Rd
@@ -10,6 +10,7 @@ parplot(object, ...)
\method{parplot}{multistart.saem.mmkin}(
object,
llmin = -Inf,
+ llquant = NA,
scale = c("best", "median"),
lpos = "bottomleft",
main = "",
@@ -23,7 +24,11 @@ parplot(object, ...)
\item{llmin}{The minimum likelihood of objects to be shown}
-\item{scale}{By default, scale parameters using the best available fit.
+\item{llquant}{Fractional value for selecting only the fits with higher
+likelihoods. Overrides 'llmin'.}
+
+\item{scale}{By default, scale parameters using the best
+available fit.
If 'median', parameters are scaled using the median parameters from all fits.}
\item{lpos}{Positioning of the legend.}
diff --git a/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg
index 7017908e..b01dac74 100644
--- a/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg
+++ b/tests/testthat/_snaps/multistart/parplot-for-dfop-sfo-fit.svg
@@ -25,109 +25,104 @@
-
-
-
-
+
+
+
+
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-
-
-
-
-
+
+
+
+
+
-
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
-
+
+
+
+
+
-
-
-
-
-
-
+
+
+
+
+
+
-
-
-
-
-
+
+
+
+
+
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
-
-
-
-
-
+
+
+
+
+
-
-
-
-
-
-
+
+
+
+
+
+
-
-
-
-
-
-
+
+
+
+
+
+
-
-
-
-
-
-
+
+
+
+
+
+
-
-
-
-
-
+
+
+
+
+
-
+
@@ -194,8 +189,8 @@
-Starting parameters
-Original run
+Original start
+Original resultsMultistart runs
diff --git a/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
index a47a585a..c8d4970b 100644
--- a/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
+++ b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
@@ -103,8 +103,8 @@
-Starting parameters
-Original run
+Original start
+Original resultsMultistart runs
diff --git a/tests/testthat/test_multistart.R b/tests/testthat/test_multistart.R
index 91ef71f0..dda0ea23 100644
--- a/tests/testthat/test_multistart.R
+++ b/tests/testthat/test_multistart.R
@@ -50,7 +50,7 @@ test_that("multistart works for saem.mmkin models", {
llhist_dfop_sfo <- function() llhist(saem_dfop_sfo_m_multi)
parplot_dfop_sfo <- function() parplot(saem_dfop_sfo_m_multi,
- ylim = c(0.5, 2))
+ ylim = c(0.5, 2), llquant = 0.5)
vdiffr::expect_doppelganger("llhist for dfop sfo fit", llhist_dfop_sfo)
vdiffr::expect_doppelganger("parplot for dfop sfo fit", parplot_dfop_sfo)
diff --git a/vignettes/FOCUS_D.html b/vignettes/FOCUS_D.html
index 9f768b06..b8a63a7b 100644
--- a/vignettes/FOCUS_D.html
+++ b/vignettes/FOCUS_D.html
@@ -299,8 +299,8 @@ pre code {
border-radius: 4px;
}
-.tabset-dropdown > .nav-tabs > li.active:before {
- content: "";
+.tabset-dropdown > .nav-tabs > li.active:before, .tabset-dropdown > .nav-tabs.nav-tabs-open:before {
+ content: "\e259";
font-family: 'Glyphicons Halflings';
display: inline-block;
padding: 10px;
@@ -308,16 +308,9 @@ pre code {
}
.tabset-dropdown > .nav-tabs.nav-tabs-open > li.active:before {
- content: "";
- border: none;
-}
-
-.tabset-dropdown > .nav-tabs.nav-tabs-open:before {
- content: "";
+ content: "\e258";
font-family: 'Glyphicons Halflings';
- display: inline-block;
- padding: 10px;
- border-right: 1px solid #ddd;
+ border: none;
}
.tabset-dropdown > .nav-tabs > li.active {
@@ -364,7 +357,7 @@ pre code {
Example evaluation of FOCUS Example Dataset D
Johannes Ranke
-
Last change 31 January 2019 (rebuilt 2022-07-08)
+
Last change 31 January 2019 (rebuilt 2022-12-06)
@@ -438,10 +431,10 @@ print(FOCUS_2006_D)
A comprehensive report of the results is obtained using the summary method for mkinfit objects.
summary(fit)
-
## mkin version used for fitting: 1.1.0
-## R version used for fitting: 4.2.1
-## Date of fit: Fri Jul 8 15:44:37 2022
-## Date of summary: Fri Jul 8 15:44:38 2022
+
## mkin version used for fitting: 1.2.2
+## R version used for fitting: 4.2.2
+## Date of fit: Tue Dec 6 09:39:42 2022
+## Date of summary: Tue Dec 6 09:39:42 2022
##
## Equations:
## d_parent/dt = - k_parent * parent
@@ -449,7 +442,7 @@ print(FOCUS_2006_D)
##
## Model predictions using solution type analytical
##
-## Fitted using 401 model solutions performed in 0.13 s
+## Fitted using 401 model solutions performed in 0.158 s
##
## Error model: Constant variance
##
@@ -492,11 +485,11 @@ print(FOCUS_2006_D)
Since mkin version 0.9-32 (July 2014), we can use shorthand notation like "SFO" for parent only degradation models. The following two lines fit the model and produce the summary report of the model fit. This covers the numerical analysis given in the FOCUS report.
## mkin version used for fitting: 1.1.2
-## R version used for fitting: 4.2.1
-## Date of fit: Wed Sep 14 22:28:35 2022
-## Date of summary: Wed Sep 14 22:28:35 2022
+
## mkin version used for fitting: 1.2.2
+## R version used for fitting: 4.2.2
+## Date of fit: Tue Dec 6 09:39:45 2022
+## Date of summary: Tue Dec 6 09:39:45 2022
##
## Equations:
## d_parent/dt = - k_parent * parent
##
## Model predictions using solution type analytical
##
-## Fitted using 133 model solutions performed in 0.032 s
+## Fitted using 133 model solutions performed in 0.033 s
##
## Error model: Constant variance
##
@@ -1637,10 +1630,10 @@ summary(m.L1.SFO)
## Warning in sqrt(1/diag(V)): NaNs produced
## Warning in cov2cor(ans$covar): diag(.) had 0 or NA entries; non-finite result is
## doubtful
-
## mkin version used for fitting: 1.1.2
-## R version used for fitting: 4.2.1
-## Date of fit: Wed Sep 14 22:28:35 2022
-## Date of summary: Wed Sep 14 22:28:35 2022
+
## mkin version used for fitting: 1.2.2
+## R version used for fitting: 4.2.2
+## Date of fit: Tue Dec 6 09:39:45 2022
+## Date of summary: Tue Dec 6 09:39:45 2022
##
## Equations:
## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
@@ -1742,17 +1735,17 @@ plot(m.L2.FOMC, show_residuals = TRUE,
main = "FOCUS L2 - FOMC")
summary(m.L2.FOMC, data = FALSE)
-
## mkin version used for fitting: 1.1.2
-## R version used for fitting: 4.2.1
-## Date of fit: Wed Sep 14 22:28:35 2022
-## Date of summary: Wed Sep 14 22:28:35 2022
+
## mkin version used for fitting: 1.2.2
+## R version used for fitting: 4.2.2
+## Date of fit: Tue Dec 6 09:39:45 2022
+## Date of summary: Tue Dec 6 09:39:45 2022
##
## Equations:
## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
##
## Model predictions using solution type analytical
##
-## Fitted using 239 model solutions performed in 0.049 s
+## Fitted using 239 model solutions performed in 0.048 s
##
## Error model: Constant variance
##
@@ -1820,10 +1813,10 @@ plot(m.L2.DFOP, show_residuals = TRUE, show_errmin = TRUE,
main = "FOCUS L2 - DFOP")
summary(m.L2.DFOP, data = FALSE)
-
## mkin version used for fitting: 1.1.2
-## R version used for fitting: 4.2.1
-## Date of fit: Wed Sep 14 22:28:36 2022
-## Date of summary: Wed Sep 14 22:28:36 2022
+
## mkin version used for fitting: 1.2.2
+## R version used for fitting: 4.2.2
+## Date of fit: Tue Dec 6 09:39:46 2022
+## Date of summary: Tue Dec 6 09:39:46 2022
##
## Equations:
## d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -1832,7 +1825,7 @@ plot(m.L2.DFOP, show_residuals = TRUE, show_errmin = TRUE,
##
## Model predictions using solution type analytical
##
-## Fitted using 581 model solutions performed in 0.135 s
+## Fitted using 581 model solutions performed in 0.131 s
##
## Error model: Constant variance
##
@@ -1920,10 +1913,10 @@ plot(mm.L3)
The objects returned by mmkin are arranged like a matrix, with models as a row index and datasets as a column index.
We can extract the summary and plot for e.g. the DFOP fit, using square brackets for indexing which will result in the use of the summary and plot functions working on mkinfit objects.
summary(mm.L3[["DFOP", 1]])
-
## mkin version used for fitting: 1.1.2
-## R version used for fitting: 4.2.1
-## Date of fit: Wed Sep 14 22:28:36 2022
-## Date of summary: Wed Sep 14 22:28:36 2022
+
## mkin version used for fitting: 1.2.2
+## R version used for fitting: 4.2.2
+## Date of fit: Tue Dec 6 09:39:46 2022
+## Date of summary: Tue Dec 6 09:39:46 2022
##
## Equations:
## d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -1932,7 +1925,7 @@ plot(mm.L3)
##
## Model predictions using solution type analytical
##
-## Fitted using 376 model solutions performed in 0.081 s
+## Fitted using 376 model solutions performed in 0.078 s
##
## Error model: Constant variance
##
@@ -2028,17 +2021,17 @@ plot(mm.L4)
The χ2 error level of 3.3% as well as the plot suggest that the SFO model fits very well. The error level at which the χ2 test passes is slightly lower for the FOMC model. However, the difference appears negligible.
summary(mm.L4[["SFO", 1]], data = FALSE)
-
## mkin version used for fitting: 1.1.2
-## R version used for fitting: 4.2.1
-## Date of fit: Wed Sep 14 22:28:36 2022
-## Date of summary: Wed Sep 14 22:28:37 2022
+
## mkin version used for fitting: 1.2.2
+## R version used for fitting: 4.2.2
+## Date of fit: Tue Dec 6 09:39:47 2022
+## Date of summary: Tue Dec 6 09:39:47 2022
##
## Equations:
## d_parent/dt = - k_parent * parent
##
## Model predictions using solution type analytical
##
-## Fitted using 142 model solutions performed in 0.034 s
+## Fitted using 142 model solutions performed in 0.03 s
##
## Error model: Constant variance
##
@@ -2092,10 +2085,10 @@ plot(mm.L4)
## DT50 DT90
## parent 106 352
summary(mm.L4[["FOMC", 1]], data = FALSE)
-
## mkin version used for fitting: 1.1.2
-## R version used for fitting: 4.2.1
-## Date of fit: Wed Sep 14 22:28:37 2022
-## Date of summary: Wed Sep 14 22:28:37 2022
+
## mkin version used for fitting: 1.2.2
+## R version used for fitting: 4.2.2
+## Date of fit: Tue Dec 6 09:39:47 2022
+## Date of summary: Tue Dec 6 09:39:47 2022
##
## Equations:
## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
--
cgit v1.2.3
From 97f71fc3d086bd447ab3e4d19abf32bb3114085b Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Wed, 7 Dec 2022 14:41:52 +0100
Subject: Check slopes in saemix covariate models
---
NEWS.md | 2 ++
R/illparms.R | 14 +++++++++++++-
log/test.log | 41 ++++++++++++++++-------------------------
man/illparms.Rd | 13 ++++++++++++-
4 files changed, 43 insertions(+), 27 deletions(-)
(limited to 'R')
diff --git a/NEWS.md b/NEWS.md
index c7929fbe..1931142c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -6,6 +6,8 @@
- 'R/parplot.R': Possibility to select the top 'llquant' fraction of the fits for the parameter plots, and improved legend text.
+- 'R/illparms.R': Also check if confidence intervals for slope parameters in covariate models include zero. Only implemented for fits obtained with the saemix backend.
+
# mkin 1.2.1 (2022-11-19)
- '{data,R}/ds_mixed.rda': Include the test data in the package instead of generating it in 'tests/testthat/setup_script.R'. Refactor the generating code to make it consistent and update tests.
diff --git a/R/illparms.R b/R/illparms.R
index 01e75cf1..eef4bd33 100644
--- a/R/illparms.R
+++ b/R/illparms.R
@@ -20,6 +20,9 @@
#' @param random For hierarchical fits, should random effects be tested?
#' @param errmod For hierarchical fits, should error model parameters be
#' tested?
+#' @param slopes For hierarchical [saem] fits using saemix as backend,
+#' should slope parameters in the covariate model(starting with 'beta_') be
+#' tested?
#' @return For [mkinfit] or [saem] objects, a character vector of parameter
#' names. For [mmkin] or [mhmkin] objects, a matrix like object of class
#' 'illparms.mmkin' or 'illparms.mhmkin'.
@@ -92,7 +95,7 @@ print.illparms.mmkin <- function(x, ...) {
#' @rdname illparms
#' @export
-illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) {
+illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, slopes = TRUE, ...) {
if (inherits(object$so, "try-error")) {
ill_parms <- NA
} else {
@@ -106,6 +109,15 @@ illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod
ill_parms_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "est."] > 0
ill_parms <- c(ill_parms, names(which(ill_parms_errmod)))
}
+ if (slopes) {
+ if (is.null(object$so)) stop("Slope testing is only implemented for the saemix backend")
+ slope_names <- grep("^beta_", object$so@model@name.fixed, value = TRUE)
+ ci <- object$so@results@conf.int
+ rownames(ci) <- ci$name
+ slope_ci <- ci[slope_names, ]
+ ill_parms_slopes <- slope_ci[, "lower"] < 0 & slope_ci[, "estimate"] > 0
+ ill_parms <- c(ill_parms, slope_names[ill_parms_slopes])
+ }
}
class(ill_parms) <- "illparms.saem.mmkin"
return(ill_parms)
diff --git a/log/test.log b/log/test.log
index 7614b136..af2ffc41 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,23 +1,23 @@
ℹ Testing mkin
✔ | F W S OK | Context
✔ | 5 | AIC calculation
-✔ | 5 | Analytical solutions for coupled models [3.2s]
+✔ | 5 | Analytical solutions for coupled models [3.4s]
✔ | 5 | Calculation of Akaike weights
✔ | 3 | Export dataset for reading into CAKE
✔ | 12 | Confidence intervals and p-values [1.1s]
-✔ | 1 12 | Dimethenamid data from 2018 [32.0s]
+✔ | 1 12 | Dimethenamid data from 2018 [32.7s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_dmta.R:98'): Different backends get consistent results for SFO-SFO3+, dimethenamid data
Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
────────────────────────────────────────────────────────────────────────────────
-✔ | 14 | Error model fitting [4.8s]
+✔ | 14 | Error model fitting [5.0s]
✔ | 5 | Time step normalisation
✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s]
✔ | 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]
+✔ | 4 | Test fitting the decline of metabolites from their maximum [0.4s]
✔ | 1 | Fitting the logistic model [0.2s]
-✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [42.2s]
-✔ | 1 11 | Nonlinear mixed-effects models [13.3s]
+✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [43.4s]
+✔ | 1 11 | Nonlinear mixed-effects models [13.5s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_mixed.R:78'): saemix results are reproducible for biphasic fits
Reason: Fitting with saemix takes around 10 minutes when using deSolve
@@ -27,40 +27,31 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
✔ | 3 | mkinfit features [0.7s]
✔ | 8 | mkinmod model generation and printing [0.2s]
✔ | 3 | Model predictions with mkinpredict [0.4s]
-✖ | 1 11 | Multistart method for saem.mmkin models [46.7s]
-────────────────────────────────────────────────────────────────────────────────
-Failure ('test_multistart.R:56'): multistart works for saem.mmkin models
-Snapshot of `testcase` to 'multistart/parplot-for-dfop-sfo-fit.svg' has changed
-Run `testthat::snapshot_review('multistart/')` to review changes
-Backtrace:
- 1. vdiffr::expect_doppelganger("parplot for dfop sfo fit", parplot_dfop_sfo)
- at test_multistart.R:56:2
- 3. testthat::expect_snapshot_file(...)
-────────────────────────────────────────────────────────────────────────────────
+✔ | 12 | Multistart method for saem.mmkin models [50.1s]
✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.2s]
-✔ | 9 | Nonlinear mixed-effects models with nlme [9.4s]
-✔ | 15 | Plotting [10.1s]
+✔ | 9 | Nonlinear mixed-effects models with nlme [9.8s]
+✔ | 15 | Plotting [10.5s]
✔ | 4 | Residuals extracted from mkinfit models
-✔ | 1 36 | saemix parent models [73.6s]
+✔ | 1 36 | saemix parent models [85.7s]
────────────────────────────────────────────────────────────────────────────────
Skip ('test_saemix_parent.R:143'): We can also use mkin solution methods for saem
Reason: This still takes almost 2.5 minutes although we do not solve ODEs
────────────────────────────────────────────────────────────────────────────────
-✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
+✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5s]
✔ | 11 | Processing of residue series
-✔ | 10 | Fitting the SFORB model [3.7s]
+✔ | 10 | Fitting the SFORB model [3.9s]
✔ | 1 | Summaries of old mkinfit objects
✔ | 5 | Summary [0.2s]
-✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3s]
-✔ | 9 | Hypothesis tests [8.1s]
+✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
+✔ | 9 | Hypothesis tests [8.2s]
✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 260.9 s
+Duration: 280.1 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)
• This still takes almost 2.5 minutes although we do not solve ODEs (1)
-[ FAIL 1 | WARN 0 | SKIP 3 | PASS 269 ]
+[ FAIL 0 | WARN 0 | SKIP 3 | PASS 270 ]
diff --git a/man/illparms.Rd b/man/illparms.Rd
index 14be9c35..75eb18f0 100644
--- a/man/illparms.Rd
+++ b/man/illparms.Rd
@@ -22,7 +22,14 @@ illparms(object, ...)
\method{print}{illparms.mmkin}(x, ...)
-\method{illparms}{saem.mmkin}(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...)
+\method{illparms}{saem.mmkin}(
+ object,
+ conf.level = 0.95,
+ random = TRUE,
+ errmod = TRUE,
+ slopes = TRUE,
+ ...
+)
\method{print}{illparms.saem.mmkin}(x, ...)
@@ -43,6 +50,10 @@ illparms(object, ...)
\item{errmod}{For hierarchical fits, should error model parameters be
tested?}
+
+\item{slopes}{For hierarchical \link{saem} fits using saemix as backend,
+should slope parameters in the covariate model(starting with 'beta_') be
+tested?}
}
\value{
For \link{mkinfit} or \link{saem} objects, a character vector of parameter
--
cgit v1.2.3
From 904ba9668eb76eaae4960e2188134e8c88da07ee Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Wed, 7 Dec 2022 16:19:54 +0100
Subject: Fix parplot for the case of failed multistart runs
---
NEWS.md | 2 ++
R/parms.R | 2 +-
R/parplot.R | 2 +-
3 files changed, 4 insertions(+), 2 deletions(-)
(limited to 'R')
diff --git a/NEWS.md b/NEWS.md
index 1931142c..4540b517 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -8,6 +8,8 @@
- 'R/illparms.R': Also check if confidence intervals for slope parameters in covariate models include zero. Only implemented for fits obtained with the saemix backend.
+- 'R/parplot.R': Make the function work also in the case that some of the multistart runs failed.
+
# mkin 1.2.1 (2022-11-19)
- '{data,R}/ds_mixed.rda': Include the test data in the package instead of generating it in 'tests/testthat/setup_script.R'. Refactor the generating code to make it consistent and update tests.
diff --git a/R/parms.R b/R/parms.R
index bd4e479b..bb04a570 100644
--- a/R/parms.R
+++ b/R/parms.R
@@ -77,6 +77,6 @@ parms.multistart <- function(object, exclude_failed = TRUE, ...) {
successful <- which(!is.na(res[, 1]))
first_success <- successful[1]
colnames(res) <- names(parms(object[[first_success]]))
- if (exclude_failed) res <- res[successful, ]
+ if (exclude_failed[1]) res <- res[successful, ]
return(res)
}
diff --git a/R/parplot.R b/R/parplot.R
index 63306ac2..e9c18947 100644
--- a/R/parplot.R
+++ b/R/parplot.R
@@ -41,7 +41,7 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
orig <- attr(object, "orig")
orig_parms <- parms(orig)
start_degparms <- orig$mean_dp_start
- all_parms <- parms(object)
+ all_parms <- parms(object, exclude_failed = FALSE)
if (inherits(object, "multistart.saem.mmkin")) {
llfunc <- function(object) {
--
cgit v1.2.3
From a54bd290bc3884d0000c52c1b29bc557825d9eae Mon Sep 17 00:00:00 2001
From: Johannes Ranke
Date: Thu, 15 Dec 2022 14:50:28 +0100
Subject: List random effects correlations in output if any
Update docs
---
NEWS.md | 4 +-
R/intervals.R | 8 +-
R/summary.saem.mmkin.R | 16 +-
docs/dev/index.html | 21 +-
docs/dev/news/index.html | 57 +++-
docs/dev/pkgdown.yml | 4 +-
docs/dev/reference/Rplot001.png | Bin 18113 -> 1011 bytes
docs/dev/reference/Rplot002.png | Bin 38732 -> 16953 bytes
docs/dev/reference/illparms.html | 15 +-
docs/dev/reference/parplot.html | 9 +-
docs/dev/reference/saem.html | 21 +-
docs/dev/reference/summary.saem.mmkin.html | 508 ++++++++++++++++-------------
log/test.log | 20 +-
man/summary.saem.mmkin.Rd | 13 +-
14 files changed, 422 insertions(+), 274 deletions(-)
(limited to 'R')
diff --git a/NEWS.md b/NEWS.md
index 4540b517..7e65204f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object or a list with suitable dimensions as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well.
-- 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters. Avoid redundant warnings that occurred in the calculation of correlations of the fixed effects in the case that the Fisher information matrix could not be inverted.
+- 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters. Avoid redundant warnings that occurred in the calculation of correlations of the fixed effects in the case that the Fisher information matrix could not be inverted. List correlations of random effects if specified by the user in the covariance model.
- 'R/parplot.R': Possibility to select the top 'llquant' fraction of the fits for the parameter plots, and improved legend text.
@@ -10,6 +10,8 @@
- 'R/parplot.R': Make the function work also in the case that some of the multistart runs failed.
+- 'R/intervals.R': Include correlations of random effects in the model in case there are any.
+
# mkin 1.2.1 (2022-11-19)
- '{data,R}/ds_mixed.rda': Include the test data in the package instead of generating it in 'tests/testthat/setup_script.R'. Refactor the generating code to make it consistent and update tests.
diff --git a/R/intervals.R b/R/intervals.R
index 705ef6eb..fcdbaea9 100644
--- a/R/intervals.R
+++ b/R/intervals.R
@@ -78,8 +78,12 @@ intervals.saem.mmkin <- function(object, level = 0.95, backtransform = TRUE, ...
# Random effects
sdnames <- intersect(rownames(conf.int), paste("SD", pnames, sep = "."))
- ranef_ret <- as.matrix(conf.int[sdnames, c("lower", "est.", "upper")])
- rownames(ranef_ret) <- paste0(gsub("SD\\.", "sd(", sdnames), ")")
+ corrnames <- grep("^Corr.", rownames(conf.int), value = TRUE)
+ ranef_ret <- as.matrix(conf.int[c(sdnames, corrnames), c("lower", "est.", "upper")])
+ sdnames_ret <- paste0(gsub("SD\\.", "sd(", sdnames), ")")
+ corrnames_ret <- gsub("Corr\\.(.*)\\.(.*)", "corr(\\1,\\2)", corrnames)
+ rownames(ranef_ret) <- c(sdnames_ret, corrnames_ret)
+
attr(ranef_ret, "label") <- "Random effects:"
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
index 46ab548b..49b02a50 100644
--- a/R/summary.saem.mmkin.R
+++ b/R/summary.saem.mmkin.R
@@ -75,10 +75,21 @@
#' f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo)
#' print(f_saem_dfop_sfo)
#' illparms(f_saem_dfop_sfo)
-#' f_saem_dfop_sfo_2 <- update(f_saem_dfop_sfo, covariance.model = diag(c(0, 0, 1, 1, 1, 0)))
+#' f_saem_dfop_sfo_2 <- update(f_saem_dfop_sfo,
+#' no_random_effect = c("parent_0", "log_k_m1"))
#' illparms(f_saem_dfop_sfo_2)
#' intervals(f_saem_dfop_sfo_2)
#' summary(f_saem_dfop_sfo_2, data = TRUE)
+#' # Add a correlation between random effects of g and k2
+#' cov_model_3 <- f_saem_dfop_sfo_2$so@model@covariance.model
+#' cov_model_3["log_k2", "g_qlogis"] <- 1
+#' cov_model_3["g_qlogis", "log_k2"] <- 1
+#' f_saem_dfop_sfo_3 <- update(f_saem_dfop_sfo,
+#' covariance.model = cov_model_3)
+#' intervals(f_saem_dfop_sfo_3)
+#' # The correlation does not improve the fit judged by AIC and BIC, although
+#' # the likelihood is higher with the additional parameter
+#' anova(f_saem_dfop_sfo, f_saem_dfop_sfo_2, f_saem_dfop_sfo_3)
#' }
#'
#' @export
@@ -150,7 +161,8 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
# Random effects
sdnames <- intersect(rownames(conf.int), paste0("SD.", pnames))
- confint_ranef <- as.matrix(conf.int[sdnames, c("estimate", "lower", "upper")])
+ corrnames <- grep("^Corr.", rownames(conf.int), value = TRUE)
+ confint_ranef <- as.matrix(conf.int[c(sdnames, corrnames), c("estimate", "lower", "upper")])
colnames(confint_ranef)[1] <- "est."
# Error model
diff --git a/docs/dev/index.html b/docs/dev/index.html
index 2615d389..4723879e 100644
--- a/docs/dev/index.html
+++ b/docs/dev/index.html
@@ -222,12 +222,21 @@
References
-
Ranke J, Wöltjen J, Schmidt J, and Comets E (2021) Taking kinetic evaluations of degradation data to the next level with nonlinear mixed-effects models. Environments8 (8) 71 doi:10.3390/environments8080071
-
-
Ranke J, Meinecke S (2019) Error Models for the Kinetic Evaluation of Chemical Degradation Data Environments6 (12) 124 doi:10.3390/environments6120124
-
-
Ranke J, Wöltjen J, Meinecke S (2018) Comparison of software tools for kinetic evaluation of chemical degradation data Environmental Sciences Europe30 17 doi:10.1186/s12302-018-0145-1
-
+
+
+Ranke J, Wöltjen J, Schmidt J, and Comets E (2021) Taking kinetic evaluations of degradation data to the next level with nonlinear mixed-effects models. Environments8 (8) 71 doi:10.3390/environments8080071
+
+
+
+
+Ranke J, Meinecke S (2019) Error Models for the Kinetic Evaluation of Chemical Degradation Data Environments6 (12) 124 doi:10.3390/environments6120124
+
+
+
+
+Ranke J, Wöltjen J, Meinecke S (2018) Comparison of software tools for kinetic evaluation of chemical degradation data Environmental Sciences Europe30 17 doi:10.1186/s12302-018-0145-1
+
‘R/mhmkin.R’: Allow an ‘illparms.mhmkin’ object or a list with suitable dimensions as value of the argument ‘no_random_effects’, making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well.
-
‘R/summary.saem.mmkin.R’: List all initial parameter values in the summary, including random effects and error model parameters. Avoid redundant warnings that occurred in the calculation of correlations of the fixed effects in the case that the Fisher information matrix could not be inverted.
+
‘R/summary.saem.mmkin.R’: List all initial parameter values in the summary, including random effects and error model parameters. Avoid redundant warnings that occurred in the calculation of correlations of the fixed effects in the case that the Fisher information matrix could not be inverted. List correlations of random effects if specified by the user in the covariance model.
+
‘R/parplot.R’: Possibility to select the top ‘llquant’ fraction of the fits for the parameter plots, and improved legend text.
+
‘R/illparms.R’: Also check if confidence intervals for slope parameters in covariate models include zero. Only implemented for fits obtained with the saemix backend.
+
‘R/parplot.R’: Make the function work also in the case that some of the multistart runs failed.
+
‘R/intervals.R’: Include correlations of random effects in the model in case there are any.
mkin 1.2.1 (2022-11-19)
@@ -140,7 +144,8 @@
mkin 1.0.5 (2021-09-15)
-
‘dimethenamid_2018’: Correct the data for the Borstel soil. The five observations from Staudenmaier (2013) that were previously stored as “Borstel 2” are actually just a subset of the 16 observations in “Borstel 1” which is now simply “Borstel”
+
‘dimethenamid_2018’: Correct the data for the Borstel soil. The five observations from Staudenmaier (2013) that were previously stored as “Borstel 2” are actually just a subset of the 16 observations in “Borstel 1” which is now simply “Borstel”
+
mkin 1.0.4 (2021-04-20)
All plotting functions setting graphical parameters: Use on.exit() for resetting graphical parameters
@@ -149,10 +154,12 @@
mkin 1.0.3 (2021-02-15)
-
Review and update README, the ‘Introduction to mkin’ vignette and some of the help pages
+
Review and update README, the ‘Introduction to mkin’ vignette and some of the help pages
+
mkin 1.0.2 (Unreleased)
-
‘mkinfit’: Keep model names stored in ‘mkinmod’ objects, avoiding their loss in ‘gmkin’
+
‘mkinfit’: Keep model names stored in ‘mkinmod’ objects, avoiding their loss in ‘gmkin’
+
mkin 1.0.1 (2021-02-10)
‘confint.mmkin’, ‘nlme.mmkin’, ‘transform_odeparms’: Fix example code in dontrun sections that failed with current defaults
@@ -207,7 +214,8 @@
mkin 0.9.49.11 (2020-04-20)
-
Increase a test tolerance to make it pass on all CRAN check machines
+
Increase a test tolerance to make it pass on all CRAN check machines
+
mkin 0.9.49.10 (2020-04-18)
‘nlme.mmkin’: An nlme method for mmkin row objects and an associated S3 class with print, plot, anova and endpoint methods
@@ -322,7 +330,8 @@
mkin 0.9.46 (2017-07-24)
-
Remove test_FOMC_ill-defined.R as it is too platform dependent
+
Remove test_FOMC_ill-defined.R as it is too platform dependent
+
mkin 0.9.45.2 (2017-07-24)
Rename twa to max_twa_parent to avoid conflict with twa from my pfm package
@@ -334,7 +343,8 @@
mkin 0.9.45.1 (2016-12-20)
New features
-
A twa function, calculating maximum time weighted average concentrations for the parent (SFO, FOMC and DFOP).
+
A twa function, calculating maximum time weighted average concentrations for the parent (SFO, FOMC and DFOP).
+
mkin 0.9.45 (2016-12-08)
@@ -349,7 +359,8 @@
mkin 0.9.44 (2016-06-29)
Bug fixes
-
The test test_FOMC_ill-defined failed on several architectures, so the test is now skipped
+
The test test_FOMC_ill-defined failed on several architectures, so the test is now skipped
+
mkin 0.9.43 (2016-06-28)
@@ -383,7 +394,8 @@
mkin 0.9.42 (2016-03-25)
Major changes
-
Add the argument from_max_mean to mkinfit, for fitting only the decline from the maximum observed value for models with a single observed variable
+
Add the argument from_max_mean to mkinfit, for fitting only the decline from the maximum observed value for models with a single observed variable
+
Minor changes
Add plots to compiled_models vignette
@@ -403,18 +415,21 @@
Bug fixes
-print.summary.mkinfit(): Avoid an error that occurred when printing summaries generated with mkin versions before 0.9-36
+print.summary.mkinfit(): Avoid an error that occurred when printing summaries generated with mkin versions before 0.9-36
+
mkin 0.9-40 (2015-07-21)
Bug fixes
-endpoints(): For DFOP and SFORB models, where optimize() is used, make use of the fact that the DT50 must be between DT50_k1 and DT50_k2 (DFOP) or DT50_b1 and DT50_b2 (SFORB), as optimize() sometimes did not find the minimum. Likewise for finding DT90 values. Also fit on the log scale to make the function more efficient.
+endpoints(): For DFOP and SFORB models, where optimize() is used, make use of the fact that the DT50 must be between DT50_k1 and DT50_k2 (DFOP) or DT50_b1 and DT50_b2 (SFORB), as optimize() sometimes did not find the minimum. Likewise for finding DT90 values. Also fit on the log scale to make the function more efficient.
+
Internal changes
-DESCRIPTION, NAMESPACE, R/*.R: Import (from) stats, graphics and methods packages, and qualify some function calls for non-base packages installed with R to avoid NOTES made by R CMD check –as-cran with upcoming R versions.
+DESCRIPTION, NAMESPACE, R/*.R: Import (from) stats, graphics and methods packages, and qualify some function calls for non-base packages installed with R to avoid NOTES made by R CMD check –as-cran with upcoming R versions.
+
mkin 0.9-39 (2015-06-26)
@@ -426,7 +441,8 @@
Bug fixes
-mkinparplot(): Fix the x axis scaling for rate constants and formation fractions that got confused by the introduction of the t-values of transformed parameters.
+mkinparplot(): Fix the x axis scaling for rate constants and formation fractions that got confused by the introduction of the t-values of transformed parameters.
+
mkin 0.9-38 (2015-06-24)
@@ -438,7 +454,8 @@
Bug fixes
-mkinmod(): When generating the C code for the derivatives, only declare the time variable when it is needed and remove the ‘-W-no-unused-variable’ compiler flag as the C compiler used in the CRAN checks on Solaris does not know it.
+mkinmod(): When generating the C code for the derivatives, only declare the time variable when it is needed and remove the ‘-W-no-unused-variable’ compiler flag as the C compiler used in the CRAN checks on Solaris does not know it.
+
mkin 0.9-36 (2015-06-21)
@@ -451,13 +468,15 @@
Minor changes
-
Added a simple showcase vignette with an evaluation of FOCUS example dataset D
+
Added a simple showcase vignette with an evaluation of FOCUS example dataset D
+
mkin 0.9-35 (2015-05-15)
Major changes
-
Switch from RUnit to testthat for testing
+
Switch from RUnit to testthat for testing
+
Bug fixes
mkinparplot(): Avoid warnings that occurred when not all confidence intervals were available in the summary of the fit
@@ -539,13 +558,15 @@
mkin 0.9-31 (2014-07-14)
Bug fixes
-
The internal renaming of optimised parameters in Version 0.9-30 led to errors in the determination of the degrees of freedom for the chi2 error level calulations in mkinerrmin() used by the summary function.
+
The internal renaming of optimised parameters in Version 0.9-30 led to errors in the determination of the degrees of freedom for the chi2 error level calulations in mkinerrmin() used by the summary function.
+
mkin 0.9-30 (2014-07-11)
New features
-
It is now possible to use formation fractions in combination with turning off the sink in mkinmod().
+
It is now possible to use formation fractions in combination with turning off the sink in mkinmod().
+
Major changes
The original and the transformed parameters now have different names (e.g. k_parent and log_k_parent. They also differ in how many they are when we have formation fractions but no pathway to sink.
diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml
index 858f6c59..642efcde 100644
--- a/docs/dev/pkgdown.yml
+++ b/docs/dev/pkgdown.yml
@@ -1,4 +1,4 @@
-pandoc: 2.9.2.1
+pandoc: 2.17.1.1
pkgdown: 2.0.6
pkgdown_sha: ~
articles:
@@ -13,7 +13,7 @@ articles:
dimethenamid_2018: web_only/dimethenamid_2018.html
multistart: web_only/multistart.html
saem_benchmarks: web_only/saem_benchmarks.html
-last_built: 2022-12-02T13:08Z
+last_built: 2022-12-15T13:46Z
urls:
reference: https://pkgdown.jrwb.de/mkin/reference
article: https://pkgdown.jrwb.de/mkin/articles
diff --git a/docs/dev/reference/Rplot001.png b/docs/dev/reference/Rplot001.png
index 8a77fc7f..17a35806 100644
Binary files a/docs/dev/reference/Rplot001.png and b/docs/dev/reference/Rplot001.png differ
diff --git a/docs/dev/reference/Rplot002.png b/docs/dev/reference/Rplot002.png
index c1621707..27feab09 100644
Binary files a/docs/dev/reference/Rplot002.png and b/docs/dev/reference/Rplot002.png differ
diff --git a/docs/dev/reference/illparms.html b/docs/dev/reference/illparms.html
index 8fe71568..9c498e1c 100644
--- a/docs/dev/reference/illparms.html
+++ b/docs/dev/reference/illparms.html
@@ -116,7 +116,14 @@ without parameter transformations is used.
print(x, ...)# S3 method for saem.mmkin
-illparms(object, conf.level =0.95, random =TRUE, errmod =TRUE, ...)
+illparms(
+object,
+ conf.level =0.95,
+ random =TRUE,
+ errmod =TRUE,
+ slopes =TRUE,
+...
+)# S3 method for illparms.saem.mmkinprint(x, ...)
@@ -154,6 +161,12 @@ without parameter transformations is used.
For hierarchical fits, should error model parameters be
tested?
+
+
slopes
+
For hierarchical saem fits using saemix as backend,
+should slope parameters in the covariate model(starting with 'beta_') be
+tested?
+
Value
diff --git a/docs/dev/reference/parplot.html b/docs/dev/reference/parplot.html
index 9852b694..720c0b2a 100644
--- a/docs/dev/reference/parplot.html
+++ b/docs/dev/reference/parplot.html
@@ -103,6 +103,7 @@ or by their medians as proposed in the paper by Duchesne et al. (2021).
parplot(object, llmin =-Inf,
+ llquant =NA, scale =c("best", "median"), lpos ="bottomleft", main ="",
@@ -124,8 +125,14 @@ or by their medians as proposed in the paper by Duchesne et al. (2021).
The minimum likelihood of objects to be shown
+
llquant
+
Fractional value for selecting only the fits with higher
+likelihoods. Overrides 'llmin'.
+
+
scale
-
By default, scale parameters using the best available fit.
+
By default, scale parameters using the best
+available fit.
If 'median', parameters are scaled using the median parameters from all fits.
diff --git a/docs/dev/reference/saem.html b/docs/dev/reference/saem.html
index d18cb848..131b168b 100644
--- a/docs/dev/reference/saem.html
+++ b/docs/dev/reference/saem.html
@@ -432,8 +432,8 @@ using mmkin.
#> saemix version used for fitting: 3.2 #> mkin version used for pre-fitting: 1.2.2 #> R version used for fitting: 4.2.2
-#> Date of fit: Thu Nov 24 08:11:00 2022
-#> Date of summary: Thu Nov 24 08:11:01 2022
+#> Date of fit: Wed Dec 7 16:22:26 2022
+#> Date of summary: Wed Dec 7 16:22:26 2022 #>#> Equations:#> d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -448,12 +448,12 @@ using mmkin.
#>#> Model predictions using solution type analytical #>
-#> Fitted in 8.778 s
+#> Fitted in 8.508 s#> Using 300, 100 iterations and 10 chains#>#> Variance model: Constant variance #>
-#> Mean of starting values for individual parameters:
+#> Starting values for degradation parameters:#> parent_0 log_k_A1 f_parent_qlogis log_k1 log_k2 #> 93.8102 -5.3734 -0.9711 -1.8799 -4.2708 #> g_qlogis
@@ -462,6 +462,19 @@ using mmkin.
#> Fixed degradation parameter values:#> None#>
+#> Starting values for random effects (square root of initial entries in omega):
+#> parent_0 log_k_A1 f_parent_qlogis log_k1 log_k2 g_qlogis
+#> parent_0 4.941 0.000 0.0000 0.000 0.000 0.0000
+#> log_k_A1 0.000 2.551 0.0000 0.000 0.000 0.0000
+#> f_parent_qlogis 0.000 0.000 0.7251 0.000 0.000 0.0000
+#> log_k1 0.000 0.000 0.0000 1.449 0.000 0.0000
+#> log_k2 0.000 0.000 0.0000 0.000 2.228 0.0000
+#> g_qlogis 0.000 0.000 0.0000 0.000 0.000 0.7814
+#>
+#> Starting values for error model parameters:
+#> a.1
+#> 1
+#>#> Results:#>#> Likelihood computed by importance sampling
diff --git a/docs/dev/reference/summary.saem.mmkin.html b/docs/dev/reference/summary.saem.mmkin.html
index a4150959..3b5869f1 100644
--- a/docs/dev/reference/summary.saem.mmkin.html
+++ b/docs/dev/reference/summary.saem.mmkin.html
@@ -102,7 +102,7 @@ endpoints such as formation fractions and DT50 values. Optionally
# S3 method for saem.mmkin
-summary(object, data =FALSE, verbose =FALSE, distimes =TRUE, ...)
+summary(object, data =FALSE, verbose =FALSE, distimes =TRUE, ...)# S3 method for summary.saem.mmkinprint(x, digits =max(3, getOption("digits")-3), verbose =x$verbose, ...)
The purpose of this document is to demonstrate how nonlinear
+hierarchical models (NLHM) based on the parent degradation models SFO,
+FOMC, DFOP and HS can be fitted with the mkin package.
+
The mkin package is used in version 1.2.2. It contains the test data
+and the functions used in the evaluations. The saemix
+package is used as a backend for fitting the NLHM, but is also loaded to
+make the convergence plot function available.
+
This document is processed with the knitr package, which
+also provides the kable function that is used to improve
+the display of tabular data in R markdown documents. For parallel
+processing, the parallel package is used.
The test data are available in the mkin package as an object of class
+mkindsg (mkin dataset group) under the identifier
+dimethenamid_2018. The following preprocessing steps are
+still necessary:
+
+
The data available for the enantiomer dimethenamid-P (DMTAP) are
+renamed to have the same substance name as the data for the racemic
+mixture dimethenamid (DMTA). The reason for this is that no difference
+between their degradation behaviour was identified in the EU risk
+assessment.
+
The data for transformation products and unnecessary columns are
+discarded
+
The observation times of each dataset are multiplied with the
+corresponding normalisation factor also available in the dataset, in
+order to make it possible to describe all datasets with a single set of
+parameters that are independent of temperature
+
Finally, datasets observed in the same soil (Elliot 1
+and Elliot 2) are combined, resulting in dimethenamid
+(DMTA) data from six soils.
+
+
The following commented R code performs this preprocessing.
+
+# Apply a function to each of the seven datasets in the mkindsg object to create a list
+dmta_ds<-lapply(1:7, function(i){
+ds_i<-dimethenamid_2018$ds[[i]]$data# Get a dataset
+ds_i[ds_i$name=="DMTAP", "name"]<-"DMTA"# Rename DMTAP to DMTA
+ds_i<-subset(ds_i, name=="DMTA", c("name", "time", "value"))# Select data
+ds_i$time<-ds_i$time*dimethenamid_2018$f_time_norm[i]# Normalise time
+ds_i# Return the dataset
+})
+
+# Use dataset titles as names for the list elements
+names(dmta_ds)<-sapply(dimethenamid_2018$ds, function(ds)ds$title)
+
+# Combine data for Elliot soil to obtain a named list with six elements
+dmta_ds[["Elliot"]]<-rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])#
+dmta_ds[["Elliot 1"]]<-NULL
+dmta_ds[["Elliot 2"]]<-NULL
In order to obtain suitable starting parameters for the NLHM fits,
+separate fits of the four models to the data for each soil are generated
+using the mmkin function from the mkin
+package. In a first step, constant variance is assumed. Convergence is
+checked with the status function.
In the table above, OK indicates convergence, and C indicates failure
+to converge. All separate fits with constant variance converged, with
+the sole exception of the HS fit to the BBA 2.2 data. To prepare for
+fitting NLHM using the two-component error model, the separate fits are
+updated assuming two-component error.
Using the two-component error model, the one fit that did not
+converge with constant variance did converge, but other non-SFO fits
+failed to converge.
+
+
+
Hierarchichal model fits
+
+
The following code fits eight versions of hierarchical models to the
+data, using SFO, FOMC, DFOP and HS for the parent compound, and using
+either constant variance or two-component error for the error model. The
+default parameter distribution model in mkin allows for variation of all
+degradation parameters across the assumed population of soils. In other
+words, each degradation parameter is associated with a random effect as
+a first step. The mhmkin function makes it possible to fit
+all eight versions in parallel (given a sufficient number of computing
+cores being available) to save execution time.
+
Convergence plots and summaries for these fits are shown in the
+appendix.
The DFOP model is preferred here, as it has a better mechanistic
+basis for batch experiments with constant incubation conditions. Also,
+it shows the lowest AIC and BIC values in the first set of fits when
+combined with the two-component error model. Therefore, the DFOP model
+was selected for further refinements of the fits with the aim to make
+the model fully identifiable.
+
+
Parameter identifiability based on the Fisher Information
+Matrix
+
+
Using the illparms function, ill-defined statistical
+model parameters such as standard deviations of the degradation
+parameters in the population and error model parameters can be
+found.
According to the illparms function, the fitted standard
+deviation of the second kinetic rate constant k2 is
+ill-defined in both DFOP fits. This suggests that different values would
+be obtained for this standard deviation when using different starting
+values.
+
The thus identified overparameterisation is addressed by removing the
+random effect for k2 from the parameter model.
which is not the case. Below, the refined model is compared with the
+previous best model. The model without random effect for k2
+is a reduced version of the previous model. Therefore, the models are
+nested and can be compared using the likelihood ratio test. This is
+achieved with the argument test = TRUE to the
+anova function.
+
+anova(f_saem[["DFOP", "tc"]], f_saem_dfop_tc_no_ranef_k2, test =TRUE)|>
+kable(format.args =list(digits =4))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
npar
+
AIC
+
BIC
+
Lik
+
Chisq
+
Df
+
Pr(>Chisq)
+
+
+
+
f_saem_dfop_tc_no_ranef_k2
+
9
+
663.8
+
661.9
+
-322.9
+
NA
+
NA
+
NA
+
+
+
f_saem[[“DFOP”, “tc”]]
+
10
+
665.5
+
663.4
+
-322.8
+
0.2809
+
1
+
0.5961
+
+
+
+
The AIC and BIC criteria are lower after removal of the ill-defined
+random effect for k2. The p value of the likelihood ratio
+test is much greater than 0.05, indicating that the model with the
+higher likelihood (here the model with random effects for all
+degradation parameters f_saem[["DFOP", "tc"]]) does not fit
+significantly better than the model with the lower likelihood (the
+reduced model f_saem_dfop_tc_no_ranef_k2).
+
Therefore, AIC, BIC and likelihood ratio test suggest the use of the
+reduced model.
+
The convergence of the fit is checked visually.
+
+
+Convergence plot for the NLHM DFOP fit with two-component error and
+without a random effect on ‘k2’
+
+
+
All parameters appear to have converged to a satisfactory degree. The
+final fit is plotted using the plot method from the mkin package.
saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:13 2023
+Date of summary: Thu Jan 5 08:19:13 2023
+
+Equations:
+d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+ time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+ * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 4.075 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Two-component variance function
+
+Starting values for degradation parameters:
+ DMTA_0 k1 k2 g
+98.759266 0.087034 0.009933 0.930827
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k1 k2 g
+DMTA_0 98.76 0 0 0
+k1 0.00 1 0 0
+k2 0.00 0 1 0
+g 0.00 0 0 1
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 663.8 661.9 -322.9
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.228939 96.285869 100.17201
+k1 0.064063 0.033477 0.09465
+k2 0.008297 0.005824 0.01077
+g 0.953821 0.914328 0.99331
+a.1 1.068479 0.869538 1.26742
+b.1 0.029424 0.022406 0.03644
+SD.DMTA_0 2.030437 0.404824 3.65605
+SD.k1 0.594692 0.256660 0.93272
+SD.g 1.006754 0.361327 1.65218
+
+Correlation:
+ DMTA_0 k1 k2
+k1 0.0218
+k2 0.0556 0.0355
+g -0.0516 -0.0284 -0.2800
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 2.0304 0.4048 3.6560
+SD.k1 0.5947 0.2567 0.9327
+SD.g 1.0068 0.3613 1.6522
+
+Variance model:
+ est. lower upper
+a.1 1.06848 0.86954 1.26742
+b.1 0.02942 0.02241 0.03644
+
+Estimated disappearance times:
+ DT50 DT90 DT50back DT50_k1 DT50_k2
+DMTA 11.45 41.4 12.46 10.82 83.54
+
+
+
Alternative check of parameter identifiability
+
+
The parameter check used in the illparms function is
+based on a quadratic approximation of the likelihood surface near its
+optimum, which is calculated using the Fisher Information Matrix (FIM).
+An alternative way to check parameter identifiability based on a
+multistart approach has recently been implemented in mkin.
+
The graph below shows boxplots of the parameters obtained in 50 runs
+of the saem algorithm with different parameter combinations, sampled
+from the range of the parameters obtained for the individual datasets
+fitted separately using nonlinear regression.
+
+f_saem_dfop_tc_multi<-multistart(f_saem[["DFOP", "tc"]], n =50, cores =15)
+Scaled parameters from the multistart runs, full model
+
+
+
The graph clearly confirms the lack of identifiability of the
+variance of k2 in the full model. The overparameterisation
+of the model also indicates a lack of identifiability of the variance of
+parameter g.
+
The parameter boxplots of the multistart runs with the reduced model
+shown below indicate that all runs give similar results, regardless of
+the starting parameters.
+
+f_saem_dfop_tc_no_ranef_k2_multi<-multistart(f_saem_dfop_tc_no_ranef_k2,
+ n =50, cores =15)
+Scaled parameters from the multistart runs, reduced model
+
+
+
When only the parameters of the top 25% of the fits are shown (based
+on a feature introduced in mkin 1.2.2 currently under development), the
+scatter is even less as shown below.
+Scaled parameters from the multistart runs, reduced model, fits with the
+top 25% likelihood values
+
+
+
+
+
+
Conclusions
+
+
Fitting the four parent degradation models SFO, FOMC, DFOP and HS as
+part of hierarchical model fits with two different error models and
+normal distributions of the transformed degradation parameters works
+without technical problems. The biphasic models DFOP and HS gave the
+best fit to the data, but the default parameter distribution model was
+not fully identifiable. Removing the random effect for the second
+kinetic rate constant of the DFOP model resulted in a reduced model that
+was fully identifiable and showed the lowest values for the model
+selection criteria AIC and BIC. The reliability of the identification of
+all model parameters was confirmed using multiple starting values.
+
+
+
Appendix
+
+
+
Hierarchical model fit listings
+
+
+Hierarchical mkin fit of the SFO model with error model const
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:06 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - k_DMTA * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 1.09 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Constant variance
+
+Starting values for degradation parameters:
+ DMTA_0 k_DMTA
+97.2953 0.0566
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k_DMTA
+DMTA_0 97.3 0
+k_DMTA 0.0 1
+
+Starting values for error model parameters:
+a.1
+ 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 796.3 795.3 -393.2
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 97.28130 95.71113 98.8515
+k_DMTA 0.05665 0.02909 0.0842
+a.1 2.66442 2.35579 2.9731
+SD.DMTA_0 1.54776 0.15447 2.9411
+SD.k_DMTA 0.60690 0.26248 0.9513
+
+Correlation:
+ DMTA_0
+k_DMTA 0.0168
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.5478 0.1545 2.9411
+SD.k_DMTA 0.6069 0.2625 0.9513
+
+Variance model:
+ est. lower upper
+a.1 2.664 2.356 2.973
+
+Estimated disappearance times:
+ DT50 DT90
+DMTA 12.24 40.65
+
+
+
+
+Hierarchical mkin fit of the SFO model with error model tc
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:07 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - k_DMTA * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 2.441 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Two-component variance function
+
+Starting values for degradation parameters:
+ DMTA_0 k_DMTA
+96.99175 0.05603
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k_DMTA
+DMTA_0 96.99 0
+k_DMTA 0.00 1
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 798.3 797.1 -393.2
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 97.271822 95.703157 98.84049
+k_DMTA 0.056638 0.029110 0.08417
+a.1 2.660081 2.230398 3.08976
+b.1 0.001665 -0.006911 0.01024
+SD.DMTA_0 1.545520 0.145035 2.94601
+SD.k_DMTA 0.606422 0.262274 0.95057
+
+Correlation:
+ DMTA_0
+k_DMTA 0.0169
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.5455 0.1450 2.9460
+SD.k_DMTA 0.6064 0.2623 0.9506
+
+Variance model:
+ est. lower upper
+a.1 2.660081 2.230398 3.08976
+b.1 0.001665 -0.006911 0.01024
+
+Estimated disappearance times:
+ DT50 DT90
+DMTA 12.24 40.65
+
+
+
+
+Hierarchical mkin fit of the FOMC model with error model const
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:06 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - (alpha/beta) * 1/((time/beta) + 1) * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 1.156 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Constant variance
+
+Starting values for degradation parameters:
+ DMTA_0 alpha beta
+ 98.292 9.909 156.341
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 alpha beta
+DMTA_0 98.29 0 0
+alpha 0.00 1 0
+beta 0.00 0 1
+
+Starting values for error model parameters:
+a.1
+ 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 734.2 732.7 -360.1
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.3435 96.9033 99.784
+alpha 7.2007 2.5889 11.812
+beta 112.8746 34.8816 190.868
+a.1 2.0459 1.8054 2.286
+SD.DMTA_0 1.4795 0.2717 2.687
+SD.alpha 0.6396 0.1509 1.128
+SD.beta 0.6874 0.1587 1.216
+
+Correlation:
+ DMTA_0 alpha
+alpha -0.1125
+beta -0.1227 0.3632
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.4795 0.2717 2.687
+SD.alpha 0.6396 0.1509 1.128
+SD.beta 0.6874 0.1587 1.216
+
+Variance model:
+ est. lower upper
+a.1 2.046 1.805 2.286
+
+Estimated disappearance times:
+ DT50 DT90 DT50back
+DMTA 11.41 42.53 12.8
+
+
+
+
+Hierarchical mkin fit of the FOMC model with error model tc
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:07 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - (alpha/beta) * 1/((time/beta) + 1) * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 2.729 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Two-component variance function
+
+Starting values for degradation parameters:
+DMTA_0 alpha beta
+98.772 4.663 92.597
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 alpha beta
+DMTA_0 98.77 0 0
+alpha 0.00 1 0
+beta 0.00 0 1
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 720.4 718.8 -352.2
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.99136 97.26011 100.72261
+alpha 5.86312 2.57485 9.15138
+beta 88.55571 29.20889 147.90254
+a.1 1.51063 1.24384 1.77741
+b.1 0.02824 0.02040 0.03609
+SD.DMTA_0 1.57436 -0.04867 3.19739
+SD.alpha 0.59871 0.17132 1.02611
+SD.beta 0.72994 0.22849 1.23139
+
+Correlation:
+ DMTA_0 alpha
+alpha -0.1363
+beta -0.1414 0.2542
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.5744 -0.04867 3.197
+SD.alpha 0.5987 0.17132 1.026
+SD.beta 0.7299 0.22849 1.231
+
+Variance model:
+ est. lower upper
+a.1 1.51063 1.2438 1.77741
+b.1 0.02824 0.0204 0.03609
+
+Estimated disappearance times:
+ DT50 DT90 DT50back
+DMTA 11.11 42.6 12.82
+
+
+
+
+Hierarchical mkin fit of the DFOP model with error model const
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:07 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+ time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+ * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 2.007 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Constant variance
+
+Starting values for degradation parameters:
+ DMTA_0 k1 k2 g
+98.64383 0.09211 0.02999 0.76814
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k1 k2 g
+DMTA_0 98.64 0 0 0
+k1 0.00 1 0 0
+k2 0.00 0 1 0
+g 0.00 0 0 1
+
+Starting values for error model parameters:
+a.1
+ 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 711.8 710 -346.9
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.092481 96.573898 99.61106
+k1 0.062499 0.030336 0.09466
+k2 0.009065 -0.005133 0.02326
+g 0.948967 0.862079 1.03586
+a.1 1.821671 1.604774 2.03857
+SD.DMTA_0 1.677785 0.472066 2.88350
+SD.k1 0.634962 0.270788 0.99914
+SD.k2 1.033498 -0.205994 2.27299
+SD.g 1.710046 0.428642 2.99145
+
+Correlation:
+ DMTA_0 k1 k2
+k1 0.0246
+k2 0.0491 0.0953
+g -0.0552 -0.0889 -0.4795
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.678 0.4721 2.8835
+SD.k1 0.635 0.2708 0.9991
+SD.k2 1.033 -0.2060 2.2730
+SD.g 1.710 0.4286 2.9914
+
+Variance model:
+ est. lower upper
+a.1 1.822 1.605 2.039
+
+Estimated disappearance times:
+ DT50 DT90 DT50back DT50_k1 DT50_k2
+DMTA 11.79 42.8 12.88 11.09 76.46
+
+
+
+
+Hierarchical mkin fit of the DFOP model with error model tc
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:08 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+ time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+ * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 3.033 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Two-component variance function
+
+Starting values for degradation parameters:
+ DMTA_0 k1 k2 g
+98.759266 0.087034 0.009933 0.930827
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k1 k2 g
+DMTA_0 98.76 0 0 0
+k1 0.00 1 0 0
+k2 0.00 0 1 0
+g 0.00 0 0 1
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 665.5 663.4 -322.8
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.377019 96.447952 100.30609
+k1 0.064843 0.034607 0.09508
+k2 0.008895 0.006368 0.01142
+g 0.949696 0.903815 0.99558
+a.1 1.065241 0.865754 1.26473
+b.1 0.029340 0.022336 0.03634
+SD.DMTA_0 2.007754 0.387982 3.62753
+SD.k1 0.580473 0.250286 0.91066
+SD.k2 0.006105 -4.920337 4.93255
+SD.g 1.097149 0.412779 1.78152
+
+Correlation:
+ DMTA_0 k1 k2
+k1 0.0235
+k2 0.0595 0.0424
+g -0.0470 -0.0278 -0.2731
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 2.007754 0.3880 3.6275
+SD.k1 0.580473 0.2503 0.9107
+SD.k2 0.006105 -4.9203 4.9325
+SD.g 1.097149 0.4128 1.7815
+
+Variance model:
+ est. lower upper
+a.1 1.06524 0.86575 1.26473
+b.1 0.02934 0.02234 0.03634
+
+Estimated disappearance times:
+ DT50 DT90 DT50back DT50_k1 DT50_k2
+DMTA 11.36 41.32 12.44 10.69 77.92
+
+
+
+
+Hierarchical mkin fit of the HS model with error model const
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:07 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - ifelse(time <= tb, k1, k2) * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 2.004 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Constant variance
+
+Starting values for degradation parameters:
+ DMTA_0 k1 k2 tb
+97.82176 0.06931 0.02997 11.13945
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k1 k2 tb
+DMTA_0 97.82 0 0 0
+k1 0.00 1 0 0
+k2 0.00 0 1 0
+tb 0.00 0 0 1
+
+Starting values for error model parameters:
+a.1
+ 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 714 712.1 -348
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 98.16102 96.47747 99.84456
+k1 0.07876 0.05261 0.10491
+k2 0.02227 0.01706 0.02747
+tb 13.99089 -7.40049 35.38228
+a.1 1.82305 1.60700 2.03910
+SD.DMTA_0 1.88413 0.56204 3.20622
+SD.k1 0.34292 0.10482 0.58102
+SD.k2 0.19851 0.01718 0.37985
+SD.tb 1.68168 0.58064 2.78272
+
+Correlation:
+ DMTA_0 k1 k2
+k1 0.0142
+k2 0.0001 -0.0025
+tb 0.0165 -0.1256 -0.0301
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 1.8841 0.56204 3.2062
+SD.k1 0.3429 0.10482 0.5810
+SD.k2 0.1985 0.01718 0.3798
+SD.tb 1.6817 0.58064 2.7827
+
+Variance model:
+ est. lower upper
+a.1 1.823 1.607 2.039
+
+Estimated disappearance times:
+ DT50 DT90 DT50back DT50_k1 DT50_k2
+DMTA 8.801 67.91 20.44 8.801 31.13
+
+
+
+
+Hierarchical mkin fit of the HS model with error model tc
+
+
+saemix version used for fitting: 3.2
+mkin version used for pre-fitting: 1.2.2
+R version used for fitting: 4.2.2
+Date of fit: Thu Jan 5 08:19:08 2023
+Date of summary: Thu Jan 5 08:20:11 2023
+
+Equations:
+d_DMTA/dt = - ifelse(time <= tb, k1, k2) * DMTA
+
+Data:
+155 observations of 1 variable(s) grouped in 6 datasets
+
+Model predictions using solution type analytical
+
+Fitted in 3.287 s
+Using 300, 100 iterations and 9 chains
+
+Variance model: Two-component variance function
+
+Starting values for degradation parameters:
+ DMTA_0 k1 k2 tb
+98.45190 0.07525 0.02576 19.19375
+
+Fixed degradation parameter values:
+None
+
+Starting values for random effects (square root of initial entries in omega):
+ DMTA_0 k1 k2 tb
+DMTA_0 98.45 0 0 0
+k1 0.00 1 0 0
+k2 0.00 0 1 0
+tb 0.00 0 0 1
+
+Starting values for error model parameters:
+a.1 b.1
+ 1 1
+
+Results:
+
+Likelihood computed by importance sampling
+ AIC BIC logLik
+ 667.1 665 -323.6
+
+Optimised parameters:
+ est. lower upper
+DMTA_0 97.76570 95.81350 99.71791
+k1 0.05855 0.03080 0.08630
+k2 0.02337 0.01664 0.03010
+tb 31.09638 29.38289 32.80987
+a.1 1.08835 0.88590 1.29080
+b.1 0.02964 0.02257 0.03671
+SD.DMTA_0 2.04877 0.42607 3.67147
+SD.k1 0.59166 0.25621 0.92711
+SD.k2 0.30698 0.09561 0.51835
+SD.tb 0.01274 -0.10914 0.13462
+
+Correlation:
+ DMTA_0 k1 k2
+k1 0.0160
+k2 -0.0070 -0.0024
+tb -0.0668 -0.0103 -0.2013
+
+Random effects:
+ est. lower upper
+SD.DMTA_0 2.04877 0.42607 3.6715
+SD.k1 0.59166 0.25621 0.9271
+SD.k2 0.30698 0.09561 0.5183
+SD.tb 0.01274 -0.10914 0.1346
+
+Variance model:
+ est. lower upper
+a.1 1.08835 0.88590 1.29080
+b.1 0.02964 0.02257 0.03671
+
+Estimated disappearance times:
+ DT50 DT90 DT50back DT50_k1 DT50_k2
+DMTA 11.84 51.71 15.57 11.84 29.66
+
+
+
+
+
+
Hierarchical model convergence plots
+
+
+
+Convergence plot for the NLHM SFO fit with constant variance
+
+
+
+
+Convergence plot for the NLHM SFO fit with two-component error
+
+
+
+
+Convergence plot for the NLHM FOMC fit with constant variance
+
+
+
+
+Convergence plot for the NLHM FOMC fit with two-component error
+
+
+
+
+Convergence plot for the NLHM DFOP fit with constant variance
+
+
+
+
+Convergence plot for the NLHM DFOP fit with two-component error
+
+
+
+
+Convergence plot for the NLHM HS fit with constant variance
+
+
+
+
+Convergence plot for the NLHM HS fit with two-component error
+
diff --git a/docs/dev/index.html b/docs/dev/index.html
index 4723879e..993b8eea 100644
--- a/docs/dev/index.html
+++ b/docs/dev/index.html
@@ -19,11 +19,11 @@
equation models are solved using automatically generated C functions.
Heteroscedasticity can be taken into account using variance by variable or
two-component error models as described by Ranke and Meinecke (2018)
- <doi:10.3390/environments6120124>. Interfaces to several nonlinear
- mixed-effects model packages are available, some of which are described by
- Ranke et al. (2021) <doi:10.3390/environments8080071>. Please note that no
- warranty is implied for correctness of results or fitness for a particular
- purpose.">
+ <doi:10.3390/environments6120124>. Hierarchical degradation models can
+ be fitted using nonlinear mixed-effects model packages as a backend as
+ described by Ranke et al. (2021) <doi:10.3390/environments8080071>. Please
+ note that no warranty is implied for correctness of results or fitness for a
+ particular purpose.">
Hierarchical kinetics template — hierarchical_kinetics • mkin
+
+
+
diff --git a/docs/dev/reference/mkinmod.html b/docs/dev/reference/mkinmod.html
index 251215a7..145dee83 100644
--- a/docs/dev/reference/mkinmod.html
+++ b/docs/dev/reference/mkinmod.html
@@ -132,7 +132,7 @@ the source compartment.
Additionally, mkinsub() has an argument to, specifying names of
variables to which a transfer is to be assumed in the model.
If the argument use_of_ff is set to "min"
-(default) and the model for the compartment is "SFO" or "SFORB", an
+and the model for the compartment is "SFO" or "SFORB", an
additional mkinsub() argument can be sink = FALSE, effectively
fixing the flux to sink to zero.
In print.mkinmod, this argument is currently not used.
@@ -247,7 +247,7 @@ in the FOCUS and NAFTA guidance documents are used.
For kinetic models with more than one observed variable, a symbolic solution
of the system of differential equations is included in the resulting
mkinmod object in some cases, speeding up the solution.
If a C compiler is found by pkgbuild::has_compiler() and there
is more than one observed variable in the specification, C code is generated
for evaluating the differential equations, compiled using
inline::cfunction() and added to the resulting mkinmod object.
@@ -310,7 +310,7 @@ Evaluating and Calculating Degradation Kinetics in Environmental Media
parent =mkinsub("SFO", "m1", full_name ="Test compound"), m1 =mkinsub("SFO", full_name ="Metabolite M1"), name ="SFO_SFO", dll_dir =DLL_dir, unload =TRUE, overwrite =TRUE)
-#> Copied DLL from /tmp/RtmpbZbZ8Y/file8c6a9f402f42.so to /home/jranke/.local/share/mkin/SFO_SFO.so
+#> Copied DLL from /tmp/RtmpelWAOB/fileb43c31a25a86.so to /home/jranke/.local/share/mkin/SFO_SFO.so# Now we can save the model and restore it in a new sessionsaveRDS(SFO_SFO.2, file ="~/SFO_SFO.rds")# Terminate the R session here if you would like to check, and then do
@@ -363,7 +363,7 @@ Evaluating and Calculating Degradation Kinetics in Environmental Media
#> })#> return(predicted)#> }
-#> <environment: 0x55556029f678>
+#> <environment: 0x55555f013820># If we have several parallel metabolites# (compare tests/testthat/test_synthetic_data_for_UBA_2014.R)
@@ -392,7 +392,7 @@ Evaluating and Calculating Degradation Kinetics in Environmental Media
diff --git a/docs/dev/reference/summary_listing.html b/docs/dev/reference/summary_listing.html
new file mode 100644
index 00000000..876412cc
--- /dev/null
+++ b/docs/dev/reference/summary_listing.html
@@ -0,0 +1,147 @@
+
+Display the output of a summary function according to the output format — summary_listing • mkin
+
+
+
Should a new page be started after the listing? Ignored in html output
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/dev/sitemap.xml b/docs/dev/sitemap.xml
index 06f56c5d..6e51371c 100644
--- a/docs/dev/sitemap.xml
+++ b/docs/dev/sitemap.xml
@@ -3,6 +3,9 @@
https://pkgdown.jrwb.de/mkin/404.html
+
+ https://pkgdown.jrwb.de/mkin/articles/2022_wp_1.1_dmta_parent.html
+ https://pkgdown.jrwb.de/mkin/articles/FOCUS_D.html
@@ -138,6 +141,9 @@
https://pkgdown.jrwb.de/mkin/reference/get_deg_func.html
+
+ https://pkgdown.jrwb.de/mkin/reference/hierarchical_kinetics.html
+ https://pkgdown.jrwb.de/mkin/reference/illparms.html
@@ -306,6 +312,9 @@
https://pkgdown.jrwb.de/mkin/reference/summary.saem.mmkin.html
+
+ https://pkgdown.jrwb.de/mkin/reference/summary_listing.html
+ https://pkgdown.jrwb.de/mkin/reference/synthetic_data_for_UBA_2014.html
diff --git a/inst/testdata/cyantraniliprole_soil_efsa_2014.xlsx b/inst/testdata/cyantraniliprole_soil_efsa_2014.xlsx
new file mode 100644
index 00000000..3252fdf1
Binary files /dev/null and b/inst/testdata/cyantraniliprole_soil_efsa_2014.xlsx differ
diff --git a/log/check.log b/log/check.log
index aec61e33..a81475d9 100644
--- a/log/check.log
+++ b/log/check.log
@@ -7,11 +7,28 @@
* checking extension type ... Package
* this is package ‘mkin’ version ‘1.2.2’
* package encoding: UTF-8
-* checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
+* checking CRAN incoming feasibility ... NOTE
Maintainer: ‘Johannes Ranke ’
+
+Size of tarball: 6636884 bytes
* checking package namespace information ... OK
* checking package dependencies ... OK
-* checking if this is a source package ... OK
+* checking if this is a source package ... NOTE
+Found the following apparent object files/libraries:
+ vignettes/2022_wp_1/cyan_dlls/dfop_path_1.so
+ vignettes/2022_wp_1/cyan_dlls/dfop_path_2.so
+ vignettes/2022_wp_1/cyan_dlls/fomc_path_1.so
+ vignettes/2022_wp_1/cyan_dlls/fomc_path_2.so
+ vignettes/2022_wp_1/cyan_dlls/hs_path_1.so
+ vignettes/2022_wp_1/cyan_dlls/sfo_path_1.so
+ vignettes/2022_wp_1/cyan_dlls/sforb_path_1.so
+ vignettes/2022_wp_1/cyan_dlls/sforb_path_2.so
+ vignettes/2022_wp_1/dmta_dlls/m_dfop_path.so
+ vignettes/2022_wp_1/dmta_dlls/m_fomc_path.so
+ vignettes/2022_wp_1/dmta_dlls/m_hs_path.so
+ vignettes/2022_wp_1/dmta_dlls/m_sfo_path.so
+ vignettes/2022_wp_1/dmta_dlls/m_sforb_path.so
+Object files/libraries should not be included in a source package.
* checking if there is a namespace ... OK
* checking for executable files ... OK
* checking for hidden files and directories ... OK
@@ -41,14 +58,7 @@ Maintainer: ‘Johannes Ranke ’
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
-* checking R code for possible problems ... NOTE
-parplot.multistart.saem.mmkin: no visible global function definition
- for ‘quantile’
-Undefined global functions or variables:
- quantile
-Consider adding
- importFrom("stats", "quantile")
-to your NAMESPACE file.
+* checking R code for possible problems ... OK
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd line widths ... OK
@@ -64,7 +74,7 @@ to your NAMESPACE file.
* checking data for ASCII and uncompressed saves ... OK
* checking installed files from ‘inst/doc’ ... OK
* checking files in ‘vignettes’ ... OK
-* checking examples ... [11s/11s] OK
+* checking examples ... [10s/10s] OK
* checking for unstated dependencies in ‘tests’ ... OK
* checking tests ... SKIPPED
* checking for unstated dependencies in vignettes ... OK
@@ -76,7 +86,7 @@ to your NAMESPACE file.
* checking for detritus in the temp directory ... OK
* DONE
-Status: 1 NOTE
+Status: 2 NOTEs
See
‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’
for details.
diff --git a/man/hierarchical_kinetics.Rd b/man/hierarchical_kinetics.Rd
index 4bb82a4c..2a8e211c 100644
--- a/man/hierarchical_kinetics.Rd
+++ b/man/hierarchical_kinetics.Rd
@@ -23,7 +23,7 @@ provided with the mkin package.
\dontrun{
library(rmarkdown)
-draft("New analysis.rmd", template = "hierarchical_kinetics", package = "mkin")
+draft("example_analysis.rmd", template = "hierarchical_kinetics", package = "mkin")
}
}
diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd
index 87ce9016..65b5de1a 100644
--- a/man/mkinmod.Rd
+++ b/man/mkinmod.Rd
@@ -33,7 +33,7 @@ the source compartment.
Additionally, \code{\link[=mkinsub]{mkinsub()}} has an argument \code{to}, specifying names of
variables to which a transfer is to be assumed in the model.
If the argument \code{use_of_ff} is set to "min"
-(default) and the model for the compartment is "SFO" or "SFORB", an
+and the model for the compartment is "SFO" or "SFORB", an
additional \code{\link[=mkinsub]{mkinsub()}} argument can be \code{sink = FALSE}, effectively
fixing the flux to sink to zero.
In print.mkinmod, this argument is currently not used.}
diff --git a/man/summary_listing.Rd b/man/summary_listing.Rd
new file mode 100644
index 00000000..995ebd8d
--- /dev/null
+++ b/man/summary_listing.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/summary_listing.R
+\name{summary_listing}
+\alias{summary_listing}
+\alias{tex_listing}
+\alias{html_listing}
+\title{Display the output of a summary function according to the output format}
+\usage{
+summary_listing(object, caption = NULL, label = NULL, clearpage = TRUE)
+
+tex_listing(object, caption = NULL, label = NULL, clearpage = TRUE)
+
+html_listing(object, caption = NULL)
+}
+\arguments{
+\item{object}{The object for which the summary is to be listed}
+
+\item{caption}{An optional caption}
+
+\item{label}{An optional label, ignored in html output}
+
+\item{clearpage}{Should a new page be started after the listing? Ignored in html output}
+}
+\description{
+This function is intended for use in a R markdown code chunk with the chunk
+option \code{results = "asis"}.
+}
diff --git a/man/tex_listing.Rd b/man/tex_listing.Rd
deleted file mode 100644
index 2f11d211..00000000
--- a/man/tex_listing.Rd
+++ /dev/null
@@ -1,21 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tex_listing.R
-\name{tex_listing}
-\alias{tex_listing}
-\title{Wrap the output of a summary function in tex listing environment}
-\usage{
-tex_listing(object, caption = NULL, label = NULL, clearpage = TRUE)
-}
-\arguments{
-\item{object}{The object for which the summary is to be listed}
-
-\item{caption}{An optional caption}
-
-\item{label}{An optional label}
-
-\item{clearpage}{Should a new page be started after the listing?}
-}
-\description{
-This function can be used in a R markdown code chunk with the chunk
-option \code{results = "asis"}.
-}
diff --git a/vignettes/FOCUS_D.html b/vignettes/FOCUS_D.html
index b8a63a7b..c729e3c2 100644
--- a/vignettes/FOCUS_D.html
+++ b/vignettes/FOCUS_D.html
@@ -31,7 +31,7 @@ document.addEventListener('DOMContentLoaded', function(e) {
!function(e,t){"use strict";"object"==typeof module&&"object"==typeof module.exports?module.exports=e.document?t(e,!0):function(e){if(!e.document)throw new Error("jQuery requires a window with a document");return t(e)}:t(e)}("undefined"!=typeof window?window:this,function(C,e){"use strict";var t=[],r=Object.getPrototypeOf,s=t.slice,g=t.flat?function(e){return t.flat.call(e)}:function(e){return t.concat.apply([],e)},u=t.push,i=t.indexOf,n={},o=n.toString,v=n.hasOwnProperty,a=v.toString,l=a.call(Object),y={},m=function(e){return"function"==typeof e&&"number"!=typeof e.nodeType&&"function"!=typeof e.item},x=function(e){return null!=e&&e===e.window},E=C.document,c={type:!0,src:!0,nonce:!0,noModule:!0};function b(e,t,n){var r,i,o=(n=n||E).createElement("script");if(o.text=e,t)for(r in c)(i=t[r]||t.getAttribute&&t.getAttribute(r))&&o.setAttribute(r,i);n.head.appendChild(o).parentNode.removeChild(o)}function w(e){return null==e?e+"":"object"==typeof e||"function"==typeof e?n[o.call(e)]||"object":typeof e}var f="3.6.0",S=function(e,t){return new S.fn.init(e,t)};function p(e){var t=!!e&&"length"in e&&e.length,n=w(e);return!m(e)&&!x(e)&&("array"===n||0===t||"number"==typeof t&&0+~]|"+M+")"+M+"*"),U=new RegExp(M+"|>"),X=new RegExp(F),V=new RegExp("^"+I+"$"),G={ID:new RegExp("^#("+I+")"),CLASS:new RegExp("^\\.("+I+")"),TAG:new RegExp("^("+I+"|[*])"),ATTR:new RegExp("^"+W),PSEUDO:new RegExp("^"+F),CHILD:new RegExp("^:(only|first|last|nth|nth-last)-(child|of-type)(?:\\("+M+"*(even|odd|(([+-]|)(\\d*)n|)"+M+"*(?:([+-]|)"+M+"*(\\d+)|))"+M+"*\\)|)","i"),bool:new RegExp("^(?:"+R+")$","i"),needsContext:new RegExp("^"+M+"*[>+~]|:(even|odd|eq|gt|lt|nth|first|last)(?:\\("+M+"*((?:-\\d)?\\d*)"+M+"*\\)|)(?=[^-]|$)","i")},Y=/HTML$/i,Q=/^(?:input|select|textarea|button)$/i,J=/^h\d$/i,K=/^[^{]+\{\s*\[native \w/,Z=/^(?:#([\w-]+)|(\w+)|\.([\w-]+))$/,ee=/[+~]/,te=new RegExp("\\\\[\\da-fA-F]{1,6}"+M+"?|\\\\([^\\r\\n\\f])","g"),ne=function(e,t){var n="0x"+e.slice(1)-65536;return t||(n<0?String.fromCharCode(n+65536):String.fromCharCode(n>>10|55296,1023&n|56320))},re=/([\0-\x1f\x7f]|^-?\d)|^-$|[^\0-\x1f\x7f-\uFFFF\w-]/g,ie=function(e,t){return t?"\0"===e?"\ufffd":e.slice(0,-1)+"\\"+e.charCodeAt(e.length-1).toString(16)+" ":"\\"+e},oe=function(){T()},ae=be(function(e){return!0===e.disabled&&"fieldset"===e.nodeName.toLowerCase()},{dir:"parentNode",next:"legend"});try{H.apply(t=O.call(p.childNodes),p.childNodes),t[p.childNodes.length].nodeType}catch(e){H={apply:t.length?function(e,t){L.apply(e,O.call(t))}:function(e,t){var n=e.length,r=0;while(e[n++]=t[r++]);e.length=n-1}}}function se(t,e,n,r){var i,o,a,s,u,l,c,f=e&&e.ownerDocument,p=e?e.nodeType:9;if(n=n||[],"string"!=typeof t||!t||1!==p&&9!==p&&11!==p)return n;if(!r&&(T(e),e=e||C,E)){if(11!==p&&(u=Z.exec(t)))if(i=u[1]){if(9===p){if(!(a=e.getElementById(i)))return n;if(a.id===i)return n.push(a),n}else if(f&&(a=f.getElementById(i))&&y(e,a)&&a.id===i)return n.push(a),n}else{if(u[2])return H.apply(n,e.getElementsByTagName(t)),n;if((i=u[3])&&d.getElementsByClassName&&e.getElementsByClassName)return H.apply(n,e.getElementsByClassName(i)),n}if(d.qsa&&!N[t+" "]&&(!v||!v.test(t))&&(1!==p||"object"!==e.nodeName.toLowerCase())){if(c=t,f=e,1===p&&(U.test(t)||z.test(t))){(f=ee.test(t)&&ye(e.parentNode)||e)===e&&d.scope||((s=e.getAttribute("id"))?s=s.replace(re,ie):e.setAttribute("id",s=S)),o=(l=h(t)).length;while(o--)l[o]=(s?"#"+s:":scope")+" "+xe(l[o]);c=l.join(",")}try{return H.apply(n,f.querySelectorAll(c)),n}catch(e){N(t,!0)}finally{s===S&&e.removeAttribute("id")}}}return g(t.replace($,"$1"),e,n,r)}function ue(){var r=[];return function e(t,n){return r.push(t+" ")>b.cacheLength&&delete e[r.shift()],e[t+" "]=n}}function le(e){return e[S]=!0,e}function ce(e){var t=C.createElement("fieldset");try{return!!e(t)}catch(e){return!1}finally{t.parentNode&&t.parentNode.removeChild(t),t=null}}function fe(e,t){var n=e.split("|"),r=n.length;while(r--)b.attrHandle[n[r]]=t}function pe(e,t){var n=t&&e,r=n&&1===e.nodeType&&1===t.nodeType&&e.sourceIndex-t.sourceIndex;if(r)return r;if(n)while(n=n.nextSibling)if(n===t)return-1;return e?1:-1}function de(t){return function(e){return"input"===e.nodeName.toLowerCase()&&e.type===t}}function he(n){return function(e){var t=e.nodeName.toLowerCase();return("input"===t||"button"===t)&&e.type===n}}function ge(t){return function(e){return"form"in e?e.parentNode&&!1===e.disabled?"label"in e?"label"in e.parentNode?e.parentNode.disabled===t:e.disabled===t:e.isDisabled===t||e.isDisabled!==!t&&ae(e)===t:e.disabled===t:"label"in e&&e.disabled===t}}function ve(a){return le(function(o){return o=+o,le(function(e,t){var n,r=a([],e.length,o),i=r.length;while(i--)e[n=r[i]]&&(e[n]=!(t[n]=e[n]))})})}function ye(e){return e&&"undefined"!=typeof e.getElementsByTagName&&e}for(e in d=se.support={},i=se.isXML=function(e){var t=e&&e.namespaceURI,n=e&&(e.ownerDocument||e).documentElement;return!Y.test(t||n&&n.nodeName||"HTML")},T=se.setDocument=function(e){var t,n,r=e?e.ownerDocument||e:p;return r!=C&&9===r.nodeType&&r.documentElement&&(a=(C=r).documentElement,E=!i(C),p!=C&&(n=C.defaultView)&&n.top!==n&&(n.addEventListener?n.addEventListener("unload",oe,!1):n.attachEvent&&n.attachEvent("onunload",oe)),d.scope=ce(function(e){return a.appendChild(e).appendChild(C.createElement("div")),"undefined"!=typeof e.querySelectorAll&&!e.querySelectorAll(":scope fieldset div").length}),d.attributes=ce(function(e){return e.className="i",!e.getAttribute("className")}),d.getElementsByTagName=ce(function(e){return e.appendChild(C.createComment("")),!e.getElementsByTagName("*").length}),d.getElementsByClassName=K.test(C.getElementsByClassName),d.getById=ce(function(e){return a.appendChild(e).id=S,!C.getElementsByName||!C.getElementsByName(S).length}),d.getById?(b.filter.ID=function(e){var t=e.replace(te,ne);return function(e){return e.getAttribute("id")===t}},b.find.ID=function(e,t){if("undefined"!=typeof t.getElementById&&E){var n=t.getElementById(e);return n?[n]:[]}}):(b.filter.ID=function(e){var n=e.replace(te,ne);return function(e){var t="undefined"!=typeof e.getAttributeNode&&e.getAttributeNode("id");return t&&t.value===n}},b.find.ID=function(e,t){if("undefined"!=typeof t.getElementById&&E){var n,r,i,o=t.getElementById(e);if(o){if((n=o.getAttributeNode("id"))&&n.value===e)return[o];i=t.getElementsByName(e),r=0;while(o=i[r++])if((n=o.getAttributeNode("id"))&&n.value===e)return[o]}return[]}}),b.find.TAG=d.getElementsByTagName?function(e,t){return"undefined"!=typeof t.getElementsByTagName?t.getElementsByTagName(e):d.qsa?t.querySelectorAll(e):void 0}:function(e,t){var n,r=[],i=0,o=t.getElementsByTagName(e);if("*"===e){while(n=o[i++])1===n.nodeType&&r.push(n);return r}return o},b.find.CLASS=d.getElementsByClassName&&function(e,t){if("undefined"!=typeof t.getElementsByClassName&&E)return t.getElementsByClassName(e)},s=[],v=[],(d.qsa=K.test(C.querySelectorAll))&&(ce(function(e){var t;a.appendChild(e).innerHTML="",e.querySelectorAll("[msallowcapture^='']").length&&v.push("[*^$]="+M+"*(?:''|\"\")"),e.querySelectorAll("[selected]").length||v.push("\\["+M+"*(?:value|"+R+")"),e.querySelectorAll("[id~="+S+"-]").length||v.push("~="),(t=C.createElement("input")).setAttribute("name",""),e.appendChild(t),e.querySelectorAll("[name='']").length||v.push("\\["+M+"*name"+M+"*="+M+"*(?:''|\"\")"),e.querySelectorAll(":checked").length||v.push(":checked"),e.querySelectorAll("a#"+S+"+*").length||v.push(".#.+[+~]"),e.querySelectorAll("\\\f"),v.push("[\\r\\n\\f]")}),ce(function(e){e.innerHTML="";var t=C.createElement("input");t.setAttribute("type","hidden"),e.appendChild(t).setAttribute("name","D"),e.querySelectorAll("[name=d]").length&&v.push("name"+M+"*[*^$|!~]?="),2!==e.querySelectorAll(":enabled").length&&v.push(":enabled",":disabled"),a.appendChild(e).disabled=!0,2!==e.querySelectorAll(":disabled").length&&v.push(":enabled",":disabled"),e.querySelectorAll("*,:x"),v.push(",.*:")})),(d.matchesSelector=K.test(c=a.matches||a.webkitMatchesSelector||a.mozMatchesSelector||a.oMatchesSelector||a.msMatchesSelector))&&ce(function(e){d.disconnectedMatch=c.call(e,"*"),c.call(e,"[s!='']:x"),s.push("!=",F)}),v=v.length&&new RegExp(v.join("|")),s=s.length&&new RegExp(s.join("|")),t=K.test(a.compareDocumentPosition),y=t||K.test(a.contains)?function(e,t){var n=9===e.nodeType?e.documentElement:e,r=t&&t.parentNode;return e===r||!(!r||1!==r.nodeType||!(n.contains?n.contains(r):e.compareDocumentPosition&&16&e.compareDocumentPosition(r)))}:function(e,t){if(t)while(t=t.parentNode)if(t===e)return!0;return!1},j=t?function(e,t){if(e===t)return l=!0,0;var n=!e.compareDocumentPosition-!t.compareDocumentPosition;return n||(1&(n=(e.ownerDocument||e)==(t.ownerDocument||t)?e.compareDocumentPosition(t):1)||!d.sortDetached&&t.compareDocumentPosition(e)===n?e==C||e.ownerDocument==p&&y(p,e)?-1:t==C||t.ownerDocument==p&&y(p,t)?1:u?P(u,e)-P(u,t):0:4&n?-1:1)}:function(e,t){if(e===t)return l=!0,0;var n,r=0,i=e.parentNode,o=t.parentNode,a=[e],s=[t];if(!i||!o)return e==C?-1:t==C?1:i?-1:o?1:u?P(u,e)-P(u,t):0;if(i===o)return pe(e,t);n=e;while(n=n.parentNode)a.unshift(n);n=t;while(n=n.parentNode)s.unshift(n);while(a[r]===s[r])r++;return r?pe(a[r],s[r]):a[r]==p?-1:s[r]==p?1:0}),C},se.matches=function(e,t){return se(e,null,null,t)},se.matchesSelector=function(e,t){if(T(e),d.matchesSelector&&E&&!N[t+" "]&&(!s||!s.test(t))&&(!v||!v.test(t)))try{var n=c.call(e,t);if(n||d.disconnectedMatch||e.document&&11!==e.document.nodeType)return n}catch(e){N(t,!0)}return 0":{dir:"parentNode",first:!0}," ":{dir:"parentNode"},"+":{dir:"previousSibling",first:!0},"~":{dir:"previousSibling"}},preFilter:{ATTR:function(e){return e[1]=e[1].replace(te,ne),e[3]=(e[3]||e[4]||e[5]||"").replace(te,ne),"~="===e[2]&&(e[3]=" "+e[3]+" "),e.slice(0,4)},CHILD:function(e){return e[1]=e[1].toLowerCase(),"nth"===e[1].slice(0,3)?(e[3]||se.error(e[0]),e[4]=+(e[4]?e[5]+(e[6]||1):2*("even"===e[3]||"odd"===e[3])),e[5]=+(e[7]+e[8]||"odd"===e[3])):e[3]&&se.error(e[0]),e},PSEUDO:function(e){var t,n=!e[6]&&e[2];return G.CHILD.test(e[0])?null:(e[3]?e[2]=e[4]||e[5]||"":n&&X.test(n)&&(t=h(n,!0))&&(t=n.indexOf(")",n.length-t)-n.length)&&(e[0]=e[0].slice(0,t),e[2]=n.slice(0,t)),e.slice(0,3))}},filter:{TAG:function(e){var t=e.replace(te,ne).toLowerCase();return"*"===e?function(){return!0}:function(e){return e.nodeName&&e.nodeName.toLowerCase()===t}},CLASS:function(e){var t=m[e+" "];return t||(t=new RegExp("(^|"+M+")"+e+"("+M+"|$)"))&&m(e,function(e){return t.test("string"==typeof e.className&&e.className||"undefined"!=typeof e.getAttribute&&e.getAttribute("class")||"")})},ATTR:function(n,r,i){return function(e){var t=se.attr(e,n);return null==t?"!="===r:!r||(t+="","="===r?t===i:"!="===r?t!==i:"^="===r?i&&0===t.indexOf(i):"*="===r?i&&-1:\x20\t\r\n\f]*)[\x20\t\r\n\f]*\/?>(?:<\/\1>|)$/i;function j(e,n,r){return m(n)?S.grep(e,function(e,t){return!!n.call(e,t,e)!==r}):n.nodeType?S.grep(e,function(e){return e===n!==r}):"string"!=typeof n?S.grep(e,function(e){return-1)[^>]*|#([\w-]+))$/;(S.fn.init=function(e,t,n){var r,i;if(!e)return this;if(n=n||D,"string"==typeof e){if(!(r="<"===e[0]&&">"===e[e.length-1]&&3<=e.length?[null,e,null]:q.exec(e))||!r[1]&&t)return!t||t.jquery?(t||n).find(e):this.constructor(t).find(e);if(r[1]){if(t=t instanceof S?t[0]:t,S.merge(this,S.parseHTML(r[1],t&&t.nodeType?t.ownerDocument||t:E,!0)),N.test(r[1])&&S.isPlainObject(t))for(r in t)m(this[r])?this[r](t[r]):this.attr(r,t[r]);return this}return(i=E.getElementById(r[2]))&&(this[0]=i,this.length=1),this}return e.nodeType?(this[0]=e,this.length=1,this):m(e)?void 0!==n.ready?n.ready(e):e(S):S.makeArray(e,this)}).prototype=S.fn,D=S(E);var L=/^(?:parents|prev(?:Until|All))/,H={children:!0,contents:!0,next:!0,prev:!0};function O(e,t){while((e=e[t])&&1!==e.nodeType);return e}S.fn.extend({has:function(e){var t=S(e,this),n=t.length;return this.filter(function(){for(var e=0;e\x20\t\r\n\f]*)/i,he=/^$|^module$|\/(?:java|ecma)script/i;ce=E.createDocumentFragment().appendChild(E.createElement("div")),(fe=E.createElement("input")).setAttribute("type","radio"),fe.setAttribute("checked","checked"),fe.setAttribute("name","t"),ce.appendChild(fe),y.checkClone=ce.cloneNode(!0).cloneNode(!0).lastChild.checked,ce.innerHTML="",y.noCloneChecked=!!ce.cloneNode(!0).lastChild.defaultValue,ce.innerHTML="",y.option=!!ce.lastChild;var ge={thead:[1,"