diff options
Diffstat (limited to 'R/multistart.R')
| -rw-r--r-- | R/multistart.R | 91 |
1 files changed, 36 insertions, 55 deletions
diff --git a/R/multistart.R b/R/multistart.R index 52f279f0..6fd98522 100644 --- a/R/multistart.R +++ b/R/multistart.R @@ -10,24 +10,49 @@ #' Currently, parallel execution of the fits is only supported using #' [parallel::mclapply], i.e. not available on Windows. #' +#' In case the online version of this help page contains error messages +#' in the example code and no plots, this is due to the multistart method +#' not working when called by pkgdown. Please refer to the +#' [online vignette](https://pkgdown.jrwb.de/mkin/dev/articles/web_only/multistart.html) +#' in this case. +#' #' @param object The fit object to work with #' @param n How many different combinations of starting parameters should be #' used? #' @param cores How many fits should be run in parallel? -#' @param \dots Passed to the update function, or to the basic plotting -#' function in the case of the graphical functions. +#' @param \dots Passed to the update function. #' @param x The multistart object to print -#' @param breaks Passed to [hist] -#' @param main title of the plot -#' @param lpos Positioning of the legend. #' @return A list of [saem.mmkin] objects, with class attributes #' 'multistart.saem.mmkin' and 'multistart'. +#' @seealso [parhist], [llhist] #' #' @references Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical #' identifiability in the frame of nonlinear mixed effects models: the example #' of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478. #' doi: 10.1186/s12859-021-04373-4. #' @export +#' @examples +#' \dontrun{ +#' library(mkin) +#' dmta_ds <- lapply(1:7, function(i) { +#' ds_i <- dimethenamid_2018$ds[[i]]$data +#' ds_i[ds_i$name == "DMTAP", "name"] <- "DMTA" +#' ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i] +#' ds_i +#' }) +#' names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title) +#' dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]]) +#' dmta_ds[["Elliot 1"]] <- dmta_ds[["Elliot 2"]] <- NULL +#' +#' f_mmkin <- mmkin("DFOP", dmta_ds, error_model = "tc", cores = 7, quiet = TRUE) +#' f_saem_full <- saem(f_mmkin) +#' f_saem_full_multi <- multistart(f_saem_full, n = 16, cores = 16) +#' parhist(f_saem_full_multi, lpos = "bottomright") +#' +#' f_saem_reduced <- update(f_saem_full, covariance.model = diag(c(1, 1, 0, 1))) +#' f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cores = 16) +#' parhist(f_saem_reduced_multi, lpos = "bottomright") +#' } multistart <- function(object, n = 50, cores = 1, ...) { UseMethod("multistart", object) @@ -36,10 +61,13 @@ multistart <- function(object, n = 50, cores = 1, ...) #' @rdname multistart #' @export multistart.saem.mmkin <- function(object, n = 50, cores = 1, ...) { + + mmkin_parms <- parms(object$mmkin, errparms = FALSE, + transformed = object$transformations == "mkin") start_parms <- apply( - parms(object$mmkin, errparms = FALSE), 1, - function(x) stats::runif(n, min(x), max(x)) - ) + mmkin_parms, 1, + function(x) stats::runif(n, min(x), max(x))) + if (n == 1) dim(start_parms) <- c(1, length(start_parms)) res <- parallel::mclapply(1:n, function(x) { update(object, degparms_start = start_parms[x, ], ...) @@ -62,50 +90,3 @@ print.multistart <- function(x, ...) { parms.multistart <- function(object, ...) { t(sapply(object, parms)) } - -#' @rdname multistart -#' @importFrom stats median -#' @export -parhist <- function(object, lpos = "topleft", main = "", ...) { - orig <- attr(object, "orig") - orig_parms <- parms(orig) - start_parms <- orig$mean_dp_start - all_parms <- parms(object) - median_parms <- apply(all_parms, 2, median) - all_scaled_parms <- t(apply(all_parms, 1, function(x) x / median_parms)) - orig_scaled_parms <- orig_parms / median_parms - start_scaled_parms <- rep(NA_real_, length(orig_parms)) - names(start_scaled_parms) <- names(orig_parms) - start_scaled_parms[names(start_parms)] <- - start_parms / median_parms[names(start_parms)] - - boxplot(all_scaled_parms, log = "y", main = main, ...) - points(orig_scaled_parms, col = 2, cex = 2) - points(start_scaled_parms, col = 3, cex = 3) - legend(lpos, inset = c(0.05, 0.05), bty = "n", - pch = 1, col = 3:1, lty = c(NA, NA, 1), - legend = c( - "Starting parameters", - "Converged parameters", - "Multistart runs")) -} - -#' @rdname multistart -#' @importFrom KernSmooth bkde -#' @export -llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "", ...) { - ll <- sapply(object, logLik) - kde <- KernSmooth::bkde(ll) - h <- hist(ll, freq = TRUE, - xlim = range(kde$x), - xlab = "", main = main, - ylab = "Frequency of log likelihoods", breaks = breaks, ...) - - freq_factor <- h$counts[1] / h$density[1] - lines(kde$x, freq_factor * kde$y) - abline(v = logLik(attr(object, "orig")), col = 2) - legend(lpos, inset = c(0.05, 0.05), bty = "n", - lty = 1, col = c(2, 1), - legend = c("original log likelihood", - "kernel density estimate")) -} |