aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--R/mkinpredict.R4
-rw-r--r--R/saem.R97
-rw-r--r--man/saem.Rd30
-rw-r--r--test.log14
4 files changed, 78 insertions, 67 deletions
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 4b618490..3402a7ba 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -134,6 +134,8 @@ mkinpredict.mkinmod <- function(x,
dimnames = list(as.character(outtimes), c("time", obs_vars)))
out_obs[, "time"] <- outtimes
+ n_out_na <- 0 # to check if we get NA values with deSolve
+
if (solution_type == "analytical") {
# This is clumsy, as we wanted fast analytical predictions for mkinfit,
# which bypasses mkinpredict in the case of analytical solutions
@@ -222,7 +224,7 @@ mkinpredict.mkinmod <- function(x,
if (map_output) {
# Output transformation for models with unobserved compartments like SFORB
# if not already mapped in analytical solution
- if (!na_stop) {
+ if (n_out_na > 0 & !na_stop) {
available <- c(TRUE, rep(FALSE, length(outtimes) - 1))
} else {
available <- rep(TRUE, length(outtimes))
diff --git a/R/saem.R b/R/saem.R
index fda5569e..a21f9990 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -26,10 +26,6 @@ utils::globalVariables(c("predicted", "std"))
#' automatic choice is not desired
#' @param quiet Should we suppress the messages saemix prints at the beginning
#' and the end of the optimisation process?
-#' @param cores The number of cores to be used for multicore processing using
-#' [parallel::mclapply()]. Using more than 1 core is experimental and may
-#' lead to excessive forking, apparently depending on the BLAS version
-#' used.
#' @param suppressPlot Should we suppress any plotting that is done
#' by the saemix function?
#' @param control Passed to [saemix::saemix]
@@ -43,7 +39,7 @@ utils::globalVariables(c("predicted", "std"))
#' ds <- lapply(experimental_data_for_UBA_2019[6:10],
#' function(x) subset(x$data[c("name", "time", "value")]))
#' names(ds) <- paste("Dataset", 6:10)
-#' f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds, cores = 1,
+#' f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds,
#' state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE)
#' f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed)
#'
@@ -76,9 +72,10 @@ utils::globalVariables(c("predicted", "std"))
#' f_mmkin <- mmkin(list(
#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
#' ds, quiet = TRUE)
-#' # These take about five seconds each on this system, as we use
-#' # analytical solutions written for saemix. When using the analytical
-#' # solutions written for mkin this took around four minutes
+#' # saem fits of SFO-SFO and DFOP-SFO to these data take about five seconds
+#' # each on this system, as we use analytical solutions written for saemix.
+#' # When using the analytical solutions written for mkin this took around
+#' # four minutes
#' f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ])
#' f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
#' # We can use print, plot and summary methods to check the results
@@ -86,10 +83,18 @@ utils::globalVariables(c("predicted", "std"))
#' plot(f_saem_dfop_sfo)
#' summary(f_saem_dfop_sfo, data = TRUE)
#'
-#' # Using a single core, the following takes about 6 minutes as we do not have an
-#' # analytical solution. Using 10 cores it is slower instead of faster
-#' f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1)
-#' plot(f_saem_fomc)
+#' # The following takes about 6 minutes
+#' #f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve",
+#' # control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
+#'
+#' #saemix::compare.saemix(list(
+#' # f_saem_dfop_sfo$so,
+#' # f_saem_dfop_sfo_deSolve$so))
+#'
+#' # If the model supports it, we can also use eigenvalue based solutions, which
+#' # take a similar amount of time
+#' #f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen",
+#' # control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
#' }
#' @export
saem <- function(object, ...) UseMethod("saem")
@@ -115,14 +120,13 @@ saem.mmkin <- function(object,
# that we discard afterwards
tmp <- tempfile()
grDevices::png(tmp)
+ on.exit(grDevices::dev.off())
+ on.exit(unlink(tmp), add = TRUE)
}
fit_time <- system.time({
utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet)
})
- if (suppressPlot) {
- grDevices::dev.off()
- unlink(tmp)
- }
+
transparms_optim <- f_saemix@results@fixed.effects
names(transparms_optim) <- f_saemix@results@name.fixed
@@ -417,44 +421,45 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
uid <- unique(id)
- res_list <- parallel::mclapply(uid, function(i) {
- transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict
- names(transparms_optim) <- names(degparms_optim)
+ res_list <- lapply(uid, function(i) {
- odeini_optim <- transparms_optim[odeini_optim_parm_names]
- names(odeini_optim) <- gsub('_0$', '', odeini_optim_parm_names)
+ transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict
+ names(transparms_optim) <- names(degparms_optim)
- odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)]
+ odeini_optim <- transparms_optim[odeini_optim_parm_names]
+ names(odeini_optim) <- gsub('_0$', '', odeini_optim_parm_names)
- ode_transparms_optim_names <- setdiff(names(transparms_optim), odeini_optim_parm_names)
- odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model,
- transform_rates = object[[1]]$transform_rates,
- transform_fractions = object[[1]]$transform_fractions)
- odeparms <- c(odeparms_optim, odeparms_fixed)
+ odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)]
- xidep_i <- subset(xidep, id == i)
+ ode_transparms_optim_names <- setdiff(names(transparms_optim), odeini_optim_parm_names)
+ odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model,
+ transform_rates = object[[1]]$transform_rates,
+ transform_fractions = object[[1]]$transform_fractions)
+ odeparms <- c(odeparms_optim, odeparms_fixed)
- if (solution_type == "analytical") {
- out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms)
- } else {
+ xidep_i <- subset(xidep, id == i)
- i_time <- xidep_i$time
- i_name <- xidep_i$name
+ if (solution_type == "analytical") {
+ out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms)
+ } else {
- out_wide <- mkinpredict(mkin_model,
- odeparms = odeparms, odeini = odeini,
- solution_type = solution_type,
- outtimes = sort(unique(i_time)),
- na_stop = FALSE
- )
+ i_time <- xidep_i$time
+ i_name <- xidep_i$name
- out_index <- cbind(as.character(i_time), as.character(i_name))
- out_values <- out_wide[out_index]
- }
- return(out_values)
- }, mc.cores = cores)
- res <- unlist(res_list)
- return(res)
+ out_wide <- mkinpredict(mkin_model,
+ odeparms = odeparms, odeini = odeini,
+ solution_type = solution_type,
+ outtimes = sort(unique(i_time)),
+ na_stop = FALSE
+ )
+
+ out_index <- cbind(as.character(i_time), as.character(i_name))
+ out_values <- out_wide[out_index]
+ }
+ return(out_values)
+ })
+ res <- unlist(res_list)
+ return(res)
}
}
diff --git a/man/saem.Rd b/man/saem.Rd
index d8d6ea0f..1fdfb2c4 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -53,11 +53,6 @@ automatic choice is not desired}
\item{control}{Passed to \link[saemix:saemix]{saemix::saemix}}
-\item{cores}{The number of cores to be used for multicore processing using
-\code{\link[parallel:mclapply]{parallel::mclapply()}}. Using more than 1 core is experimental and may
-lead to excessive forking, apparently depending on the BLAS version
-used.}
-
\item{verbose}{Should we print information about created objects of
type \link[saemix:SaemixModel-class]{saemix::SaemixModel} and \link[saemix:SaemixData-class]{saemix::SaemixData}?}
@@ -98,7 +93,7 @@ using \link{mmkin}.
ds <- lapply(experimental_data_for_UBA_2019[6:10],
function(x) subset(x$data[c("name", "time", "value")]))
names(ds) <- paste("Dataset", 6:10)
-f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds, cores = 1,
+f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds,
state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE)
f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed)
@@ -131,9 +126,10 @@ dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
f_mmkin <- mmkin(list(
"SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
ds, quiet = TRUE)
-# These take about five seconds each on this system, as we use
-# analytical solutions written for saemix. When using the analytical
-# solutions written for mkin this took around four minutes
+# saem fits of SFO-SFO and DFOP-SFO to these data take about five seconds
+# each on this system, as we use analytical solutions written for saemix.
+# When using the analytical solutions written for mkin this took around
+# four minutes
f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ])
f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
# We can use print, plot and summary methods to check the results
@@ -141,10 +137,18 @@ print(f_saem_dfop_sfo)
plot(f_saem_dfop_sfo)
summary(f_saem_dfop_sfo, data = TRUE)
-# Using a single core, the following takes about 6 minutes as we do not have an
-# analytical solution. Using 10 cores it is slower instead of faster
-f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1)
-plot(f_saem_fomc)
+# The following takes about 6 minutes
+#f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve",
+# control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
+
+#saemix::compare.saemix(list(
+# f_saem_dfop_sfo$so,
+# f_saem_dfop_sfo_deSolve$so))
+
+# If the model supports it, we can also use eigenvalue based solutions, which
+# take a similar amount of time
+#f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen",
+# control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
}
}
\seealso{
diff --git a/test.log b/test.log
index 6105128f..13fa7e70 100644
--- a/test.log
+++ b/test.log
@@ -6,32 +6,32 @@ Testing mkin
✔ | 2 | Export dataset for reading into CAKE
✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.0 s]
✔ | 4 | Calculation of FOCUS chi2 error levels [0.5 s]
-✔ | 7 | Fitting the SFORB model [3.5 s]
-✔ | 5 | Analytical solutions for coupled models [3.2 s]
+✔ | 7 | Fitting the SFORB model [3.7 s]
+✔ | 5 | Analytical solutions for coupled models [3.3 s]
✔ | 5 | Calculation of Akaike weights
✔ | 14 | Confidence intervals and p-values [1.2 s]
✔ | 14 | Error model fitting [4.7 s]
✔ | 3 | Time step normalisation
✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3 s]
✔ | 1 | Fitting the logistic model [0.2 s]
-✔ | 31 | Nonlinear mixed effects models [27.5 s]
+✔ | 31 | Nonlinear mixed effects models [27.4 s]
✔ | 2 | Test dataset classes mkinds and mkindsg
✔ | 1 | mkinfit features [0.3 s]
✔ | 12 | Special cases of mkinfit calls [0.8 s]
✔ | 8 | mkinmod model generation and printing [0.2 s]
✔ | 3 | Model predictions with mkinpredict [0.4 s]
✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.6 s]
-✔ | 9 | Nonlinear mixed-effects models [7.8 s]
+✔ | 9 | Nonlinear mixed-effects models [7.9 s]
✔ | 16 | Plotting [1.9 s]
✔ | 4 | Residuals extracted from mkinfit models
✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s]
✔ | 4 | Summary [0.1 s]
✔ | 1 | Summaries of old mkinfit objects
✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3 s]
-✔ | 9 | Hypothesis tests [7.3 s]
-✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s]
+✔ | 9 | Hypothesis tests [7.5 s]
+✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.3 s]
══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 69.3 s
+Duration: 69.5 s
[ FAIL 0 | WARN 0 | SKIP 0 | PASS 204 ]

Contact - Imprint