diff options
| -rw-r--r-- | R/saem.R | 37 | ||||
| -rw-r--r-- | docs/dev/pkgdown.yml | 2 | ||||
| -rw-r--r-- | log/test.log | 20 | ||||
| -rw-r--r-- | man/saem.Rd | 12 | ||||
| -rw-r--r-- | tests/testthat/test_saemix_parent.R | 52 |
5 files changed, 94 insertions, 29 deletions
@@ -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") |