diff options
| -rw-r--r-- | NAMESPACE | 2 | ||||
| -rw-r--r-- | R/mhmkin.R | 4 | ||||
| -rw-r--r-- | R/multistart.R | 29 | ||||
| -rw-r--r-- | R/saem.R | 37 | ||||
| -rw-r--r-- | log/test.log | 16 | ||||
| -rw-r--r-- | tests/testthat/test_saemix_parent.R | 15 |
6 files changed, 94 insertions, 9 deletions
@@ -68,8 +68,10 @@ S3method(summary,mkinfit) S3method(summary,mmkin) S3method(summary,nlme.mmkin) S3method(summary,saem.mmkin) +S3method(update,mhmkin) S3method(update,mkinfit) S3method(update,mmkin) +S3method(update,multistart) S3method(update,nlme.mmkin) S3method(update,saem.mmkin) S3method(which.best,default) @@ -207,6 +207,10 @@ BIC.mhmkin <- function(object, ...) { #' @export update.mhmkin <- function(object, ..., evaluate = TRUE) { call <- attr(object, "call") + # For some reason we get mhkin.list in call[[1]] when using mhmkin from the + # loaded package so we need to fix this so we do not have to export + # mhmkin.list in addition to the S3 method mhmkin + call[[1]] <- mhmkin update_arguments <- match.call(expand.dots = FALSE)$... diff --git a/R/multistart.R b/R/multistart.R index a788953e..11736670 100644 --- a/R/multistart.R +++ b/R/multistart.R @@ -71,6 +71,7 @@ multistart <- function(object, n = 50, #' @export multistart.saem.mmkin <- function(object, n = 50, cores = 1, cluster = NULL, ...) { + call <- match.call() if (n <= 1) stop("Please specify an n of at least 2") mmkin_parms <- parms(object$mmkin, errparms = FALSE, @@ -90,6 +91,7 @@ multistart.saem.mmkin <- function(object, n = 50, cores = 1, } attr(res, "orig") <- object attr(res, "start_parms") <- start_parms + attr(res, "call") <- call class(res) <- c("multistart.saem.mmkin", "multistart") return(res) } @@ -178,3 +180,30 @@ which.best.default <- function(object, ...) ll <- sapply(object, llfunc) return(which.max(ll)) } + +#' @export +update.multistart <- function(object, ..., evaluate = TRUE) { + call <- attr(object, "call") + # For some reason we get multistart.saem.mmkin in call[[1]] when using multistart + # from the loaded package so we need to fix this so we do not have to export + # multistart.saem.mmkin + call[[1]] <- multistart + + update_arguments <- match.call(expand.dots = FALSE)$... + + if (length(update_arguments) > 0) { + update_arguments_in_call <- !is.na(match(names(update_arguments), names(call))) + } + + for (a in names(update_arguments)[update_arguments_in_call]) { + call[[a]] <- update_arguments[[a]] + } + + update_arguments_not_in_call <- !update_arguments_in_call + if(any(update_arguments_not_in_call)) { + call <- c(as.list(call), update_arguments[update_arguments_not_in_call]) + call <- as.call(call) + } + if(evaluate) eval(call, parent.frame()) + else call +} @@ -313,7 +313,7 @@ saemix_model <- function(object, solution_type = "auto", # Parent only if (length(mkin_model$spec) == 1) { parent_type <- mkin_model$spec[[1]]$type - if (length(odeini_fixed) == 1) { + if (length(odeini_fixed) == 1 && !grepl("_bound$", names(odeini_fixed))) { if (transformations == "saemix") { stop("saemix transformations are not supported for parent fits with fixed initial parent value") } @@ -344,6 +344,9 @@ saemix_model <- function(object, solution_type = "auto", } } } else { + if (length(odeini_fixed) == 2) { + stop("SFORB with fixed initial parent value is not supported") + } if (parent_type == "SFO") { if (transformations == "mkin") { model_function <- function(psi, id, xidep) { @@ -386,6 +389,38 @@ saemix_model <- function(object, solution_type = "auto", transform.par = c(0, 1, 1, 3) } } + if (parent_type == "SFORB") { + if (transformations == "mkin") { + model_function <- function(psi, id, xidep) { + k_12 <- exp(psi[id, 3]) + k_21 <- exp(psi[id, 4]) + k_1output <- exp(psi[id, 2]) + t <- xidep[, "time"] + + sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21) + b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp + b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp + + psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + + ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)) + } + } else { + model_function <- function(psi, id, xidep) { + k_12 <- psi[id, 3] + k_21 <- psi[id, 4] + k_1output <- psi[id, 2] + t <- xidep[, "time"] + + sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21) + b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp + b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp + + psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + + ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)) + } + transform.par = c(0, 1, 1, 1) + } + } if (parent_type == "HS") { if (transformations == "mkin") { model_function <- function(psi, id, xidep) { diff --git a/log/test.log b/log/test.log index 6cd9e6a8..d7d623d7 100644 --- a/log/test.log +++ b/log/test.log @@ -1,11 +1,11 @@ ℹ 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.0s] -✔ | 1 12 | Dimethenamid data from 2018 [31.1s] +✔ | 1 12 | Dimethenamid data from 2018 [31.3s] ──────────────────────────────────────────────────────────────────────────────── 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 @@ -28,24 +28,24 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve ✔ | 8 | mkinmod model generation and printing [0.2s] ✔ | 3 | Model predictions with mkinpredict [0.4s] ✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.8s] -✔ | 9 | Nonlinear mixed-effects models with nlme [8.5s] -✔ | 16 | Plotting [9.8s] +✔ | 9 | Nonlinear mixed-effects models with nlme [8.3s] +✔ | 16 | Plotting [9.7s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 35 | saemix parent models [189.7s] +✔ | 37 | saemix parent models [208.2s] ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s] ✔ | 11 | Processing of residue series ✔ | 7 | Fitting the SFORB model [3.6s] ✔ | 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 [7.8s] +✔ | 9 | Hypothesis tests [7.6s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.1s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 286.0 s +Duration: 304.2 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 253 ] +[ FAIL 0 | WARN 0 | SKIP 2 | PASS 255 ] diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R index ce776bf7..0579f22f 100644 --- a/tests/testthat/test_saemix_parent.R +++ b/tests/testthat/test_saemix_parent.R @@ -101,6 +101,21 @@ test_that("Parent fits using saemix are correctly implemented", { rel_diff_2 <- (s_dfop_s2$confint_back[, "est."] - dfop_pop) / dfop_pop expect_true(all(rel_diff_2 < 0.2)) + # We use constant error for SFORB because tc is overparameterised (b.1 is ill-defined in saem) + mmkin_sforb_2 <- mmkin("SFORB", ds_dfop, quiet = TRUE, error_model = "const", cores = n_cores) + sforb_saemix_1 <- saem(mmkin_sforb_2, quiet = TRUE, + no_random_effect = c("parent_free_0", "k_parent_free_bound"), + transformations = "saemix") + sforb_saemix_2 <- saem(mmkin_sforb_2, quiet = TRUE, + no_random_effect = c("parent_free_0", "log_k_parent_free_bound"), + transformations = "mkin") + expect_equal( + log(endpoints(dfop_saemix_1)$distimes[1:2]), + log(endpoints(sforb_saemix_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) + mmkin_hs_1 <- mmkin("HS", ds_hs, quiet = TRUE, error_model = "const", cores = n_cores) hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE) hs_saem_2 <- saem(mmkin_hs_1, quiet = TRUE, transformations = "mkin") |