aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2022-10-17 10:28:54 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2022-10-17 10:28:54 +0200
commitb848fb360aa865c37298ee7526344b5280c700cc (patch)
treef1403f49672e01baf5f6b6475db6a383b0d60bee
parentc03fa5d4e57033869cb437c1154da31abd96fc50 (diff)
SFORB in saem, update for mhmkin and multistart
-rw-r--r--NAMESPACE2
-rw-r--r--R/mhmkin.R4
-rw-r--r--R/multistart.R29
-rw-r--r--R/saem.R37
-rw-r--r--log/test.log16
-rw-r--r--tests/testthat/test_saemix_parent.R15
6 files changed, 94 insertions, 9 deletions
diff --git a/NAMESPACE b/NAMESPACE
index 3b4d2367..80fc927d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -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)
diff --git a/R/mhmkin.R b/R/mhmkin.R
index 2df01f0c..542bddf7 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -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
+}
diff --git a/R/saem.R b/R/saem.R
index 5256f6b5..dbc19dec 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -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")

Contact - Imprint