aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--R/saem.R37
-rw-r--r--docs/dev/pkgdown.yml2
-rw-r--r--log/test.log20
-rw-r--r--man/saem.Rd12
-rw-r--r--tests/testthat/test_saemix_parent.R52
5 files changed, 94 insertions, 29 deletions
diff --git a/R/saem.R b/R/saem.R
index 8929cf6f..05cce682 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -33,6 +33,12 @@ utils::globalVariables(c("predicted", "std"))
#' for parameter that are tested if requested by 'test_log_parms'.
#' @param solution_type Possibility to specify the solution type in case the
#' automatic choice is not desired
+#' @param no_random_effect Character vector of degradation parameters for
+#' which there should be no variability over the groups. Only used
+#' if the covariance model is not explicitly specified.
+#' @param covariance.model Will be passed to [saemix::SaemixModel()]. Per
+#' default, uncorrelated random effects are specified for all degradation
+#' parameters.
#' @param fail_with_errors Should a failure to compute standard errors
#' from the inverse of the Fisher Information Matrix be a failure?
#' @param quiet Should we suppress the messages saemix prints at the beginning
@@ -121,6 +127,8 @@ saem.mmkin <- function(object,
test_log_parms = TRUE,
conf.level = 0.6,
solution_type = "auto",
+ covariance.model = "auto",
+ no_random_effect = NULL,
nbiter.saemix = c(300, 100),
control = list(displayProgress = FALSE, print = FALSE,
nbiter.saemix = nbiter.saemix,
@@ -134,7 +142,10 @@ saem.mmkin <- function(object,
degparms_start = degparms_start,
test_log_parms = test_log_parms, conf.level = conf.level,
solution_type = solution_type,
- transformations = transformations, ...)
+ transformations = transformations,
+ covariance.model = covariance.model,
+ no_random_effect = no_random_effect,
+ ...)
d_saemix <- saemix_data(object, verbose = verbose)
fit_failed <- FALSE
@@ -253,14 +264,22 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
#' @rdname saem
#' @return An [saemix::SaemixModel] object.
#' @export
-saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"),
- degparms_start = numeric(), test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...)
+saemix_model <- function(object, solution_type = "auto",
+ transformations = c("mkin", "saemix"), degparms_start = numeric(),
+ covariance.model = "auto", no_random_effect = NULL,
+ test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...)
{
if (nrow(object) > 1) stop("Only row objects allowed")
mkin_model <- object[[1]]$mkinmod
degparms_optim <- mean_degparms(object, test_log_parms = test_log_parms)
+ na_degparms <- names(which(is.na(degparms_optim)))
+ if (length(na_degparms) > 0) {
+ message("Did not find valid starting values for ", paste(na_degparms, collapse = ", "), "\n",
+ "Now trying with test_log_parms = FALSE")
+ degparms_optim <- mean_degparms(object, test_log_parms = FALSE)
+ }
if (transformations == "saemix") {
degparms_optim <- backtransform_odeparms(degparms_optim,
object[[1]]$mkinmod,
@@ -378,7 +397,7 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
}
} else {
model_function <- function(psi, id, xidep) {
- tb <- exp(psi[id, 4])
+ tb <- psi[id, 4]
t <- xidep[, "time"]
psi[id, 1] * ifelse(t <= tb,
exp(- psi[id, 2] * t),
@@ -560,12 +579,22 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
psi0_matrix <- matrix(degparms_psi0, nrow = 1)
colnames(psi0_matrix) <- names(degparms_psi0)
+ if (covariance.model[1] == "auto") {
+ covariance_diagonal <- rep(1, length(degparms_optim))
+ if (!is.null(no_random_effect)) {
+ degparms_no_random <- which(names(degparms_psi0) %in% no_random_effect)
+ covariance_diagonal[degparms_no_random] <- 0
+ }
+ covariance.model = diag(covariance_diagonal)
+ }
+
res <- saemix::saemixModel(model_function,
psi0 = psi0_matrix,
"Mixed model generated from mmkin object",
transform.par = transform.par,
error.model = error.model,
verbose = verbose,
+ covariance.model = covariance.model,
...
)
attr(res, "mean_dp_start") <- degparms_optim
diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml
index ee93ec22..22ed5aa7 100644
--- a/docs/dev/pkgdown.yml
+++ b/docs/dev/pkgdown.yml
@@ -12,7 +12,7 @@ articles:
compiled_models: web_only/compiled_models.html
dimethenamid_2018: web_only/dimethenamid_2018.html
multistart: web_only/multistart.html
-last_built: 2022-09-28T14:45Z
+last_built: 2022-09-30T10:35Z
urls:
reference: https://pkgdown.jrwb.de/mkin/reference
article: https://pkgdown.jrwb.de/mkin/articles
diff --git a/log/test.log b/log/test.log
index df663253..942ed50d 100644
--- a/log/test.log
+++ b/log/test.log
@@ -5,7 +5,7 @@
✔ | 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.9s]
+✔ | 1 12 | Dimethenamid data from 2018 [31.6s]
────────────────────────────────────────────────────────────────────────────────
Skip (test_dmta.R:98:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data
Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
@@ -16,7 +16,7 @@ Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
✔ | 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]
-✔ | 7 | Batch fitting and diagnosing hierarchical kinetic models [14.3s]
+✔ | 7 | Batch fitting and diagnosing hierarchical kinetic models [14.5s]
✔ | 1 12 | Nonlinear mixed-effects models [0.3s]
────────────────────────────────────────────────────────────────────────────────
Skip (test_mixed.R:68:3): saemix results are reproducible for biphasic fits
@@ -26,26 +26,26 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
✔ | 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.4s]
-✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.8s]
+✔ | 3 | Model predictions with mkinpredict [0.3s]
+✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.7s]
✔ | 9 | Nonlinear mixed-effects models with nlme [8.6s]
-✔ | 16 | Plotting [9.9s]
+✔ | 16 | Plotting [9.7s]
✔ | 4 | Residuals extracted from mkinfit models
-✔ | 28 | saemix parent models [181.0s]
+✔ | 35 | saemix parent models [191.8s]
✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
✔ | 11 | Processing of residue series
✔ | 7 | 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.1s]
-✔ | 9 | Hypothesis tests [8.0s]
+✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
+✔ | 9 | Hypothesis tests [7.8s]
✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 279.0 s
+Duration: 288.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)
-[ FAIL 0 | WARN 0 | SKIP 2 | PASS 246 ]
+[ FAIL 0 | WARN 0 | SKIP 2 | PASS 253 ]
diff --git a/man/saem.Rd b/man/saem.Rd
index dd787475..a9c4a1bb 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -18,6 +18,8 @@ saem(object, ...)
test_log_parms = TRUE,
conf.level = 0.6,
solution_type = "auto",
+ covariance.model = "auto",
+ no_random_effect = NULL,
nbiter.saemix = c(300, 100),
control = list(displayProgress = FALSE, print = FALSE, nbiter.saemix = nbiter.saemix,
save = FALSE, save.graphs = FALSE),
@@ -34,6 +36,8 @@ saemix_model(
solution_type = "auto",
transformations = c("mkin", "saemix"),
degparms_start = numeric(),
+ covariance.model = "auto",
+ no_random_effect = NULL,
test_log_parms = FALSE,
conf.level = 0.6,
verbose = FALSE,
@@ -70,6 +74,14 @@ for parameter that are tested if requested by 'test_log_parms'.}
\item{solution_type}{Possibility to specify the solution type in case the
automatic choice is not desired}
+\item{covariance.model}{Will be passed to \code{\link[saemix:SaemixModel-class]{saemix::SaemixModel()}}. Per
+default, uncorrelated random effects are specified for all degradation
+parameters.}
+
+\item{no_random_effect}{Character vector of degradation parameters for
+which there should be no variability over the groups. Only used
+if the covariance model is not explicitly specified.}
+
\item{nbiter.saemix}{Convenience option to increase the number of
iterations}
diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R
index 5d9a01de..ce776bf7 100644
--- a/tests/testthat/test_saemix_parent.R
+++ b/tests/testthat/test_saemix_parent.R
@@ -4,30 +4,54 @@ test_that("Parent fits using saemix are correctly implemented", {
skip_on_cran()
expect_error(saem(fits), "Only row objects")
- # Some fits were done in the setup script
+
+ # mmkin_sfo_1 was generated in the setup script
+ # We did not introduce variance of parent_0 in the data generation
+ # This is correctly detected
+ expect_equal(illparms(sfo_saem_1), "sd(parent_0)")
+ # So we remove this variance
+ sfo_saem_1_reduced <- update(sfo_saem_1, no_random_effect = "parent_0")
+ expect_equal(illparms(sfo_saem_1_reduced), character(0))
+
+ # We cannot currently do the fit with completely fixed initial values
mmkin_sfo_2 <- update(mmkin_sfo_1, fixed_initials = c(parent = 100))
+ sfo_saem_3 <- expect_error(saem(mmkin_sfo_2, quiet = TRUE), "at least two parameters")
+
+ # We get an error if we do not supply a suitable model specification
expect_error(update(mmkin_sfo_1, models = c("SFOOO")), "Please supply models.*")
- sfo_saem_2 <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "mkin")
- sfo_saem_3 <- expect_error(saem(mmkin_sfo_2, quiet = TRUE), "at least two parameters")
- expect_equal(endpoints(sfo_saem_1), endpoints(sfo_saem_2), tol = 0.01)
- s_sfo_s1 <- summary(sfo_saem_1)
- s_sfo_s2 <- summary(sfo_saem_2)
+ sfo_saem_1_mkin <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "mkin")
+ expect_equal(illparms(sfo_saem_1_mkin), "sd(parent_0)")
+ sfo_saem_1_reduced_mkin <- update(sfo_saem_1_mkin, no_random_effect = "parent_0")
+
+ # The endpoints obtained do not depend on the transformation
+ expect_equal(endpoints(sfo_saem_1), endpoints(sfo_saem_1_mkin), tol = 0.01)
+ expect_equal(endpoints(sfo_saem_1_reduced), endpoints(sfo_saem_1_reduced_mkin), tol = 0.01)
+
+ s_sfo_saem_1 <- summary(sfo_saem_1)
+ s_sfo_saem_1_reduced <- summary(sfo_saem_1_reduced)
+ s_sfo_saem_1_mkin <- summary(sfo_saem_1_mkin)
+ s_sfo_saem_1_reduced_mkin <- summary(sfo_saem_1_reduced_mkin)
sfo_nlme_1 <- expect_warning(nlme(mmkin_sfo_1), "not converge")
- s_sfo_n <- summary(sfo_nlme_1)
+ s_sfo_nlme_1 <- summary(sfo_nlme_1)
# Compare with input
- expect_equal(round(s_sfo_s2$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)
+ expect_equal(round(s_sfo_saem_1_mkin$confint_ranef["SD.log_k_parent", "est."], 1), 0.3)
# k_parent is a bit different from input 0.03 here
- expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), 0.035)
- expect_equal(round(s_sfo_s2$confint_back["k_parent", "est."], 3), 0.035)
+ 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)
# But the result is pretty unanimous between methods
- expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3),
- round(s_sfo_s2$confint_back["k_parent", "est."], 3))
- expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3),
- round(s_sfo_n$confint_back["k_parent", "est."], 3))
+ expect_equal(round(s_sfo_saem_1_reduced$confint_back["k_parent", "est."], 3),
+ round(s_sfo_saem_1$confint_back["k_parent", "est."], 3))
+ expect_equal(round(s_sfo_saem_1_mkin$confint_back["k_parent", "est."], 3),
+ round(s_sfo_saem_1$confint_back["k_parent", "est."], 3))
+ expect_equal(round(s_sfo_saem_1_reduced_mkin$confint_back["k_parent", "est."], 3),
+ round(s_sfo_saem_1$confint_back["k_parent", "est."], 3))
+ expect_equal(round(s_sfo_nlme_1$confint_back["k_parent", "est."], 3),
+ round(s_sfo_saem_1$confint_back["k_parent", "est."], 3))
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")

Contact - Imprint