diff options
| -rw-r--r-- | R/mkinpredict.R | 4 | ||||
| -rw-r--r-- | R/saem.R | 97 | ||||
| -rw-r--r-- | man/saem.Rd | 30 | ||||
| -rw-r--r-- | test.log | 14 |
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)) @@ -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{ @@ -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 ] |