diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/llhist.R | 7 | ||||
| -rw-r--r-- | R/multistart.R | 4 | ||||
| -rw-r--r-- | R/parhist.R | 58 |
3 files changed, 46 insertions, 23 deletions
@@ -1,6 +1,6 @@ #' Plot the distribution of log likelihoods from multistart objects #' -#' Produces a histogram of log-likelihoods, and an overlayed kernel density +#' Produces a histogram of log-likelihoods, and an overlayed kernel density #' estimate. In addition, the likelihood of the original fit is shown as #' a red vertical line. #' @@ -13,8 +13,13 @@ #' @importFrom KernSmooth bkde #' @export llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "", ...) { + oldpar <- par(no.readonly = TRUE) + on.exit(par(oldpar, no.readonly = TRUE)) + ll <- sapply(object, logLik) kde <- KernSmooth::bkde(ll) + + par(las = 1) h <- hist(ll, freq = TRUE, xlim = range(kde$x), xlab = "", main = main, diff --git a/R/multistart.R b/R/multistart.R index 6fd98522..beb8194b 100644 --- a/R/multistart.R +++ b/R/multistart.R @@ -12,7 +12,7 @@ #' #' 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 +#' 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. #' @@ -61,13 +61,13 @@ multistart <- function(object, n = 50, cores = 1, ...) #' @rdname multistart #' @export multistart.saem.mmkin <- function(object, n = 50, cores = 1, ...) { + if (n <= 1) stop("Please specify an n of at least 2") mmkin_parms <- parms(object$mmkin, errparms = FALSE, transformed = object$transformations == "mkin") start_parms <- apply( 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, ], ...) diff --git a/R/parhist.R b/R/parhist.R index 62f67249..69aafe02 100644 --- a/R/parhist.R +++ b/R/parhist.R @@ -1,11 +1,7 @@ #' Plot parameter distributions from multistart objects #' -#' Produces a boxplot with all parameters from the multiple runs, scaled using -#' their medians. If parameter transformations were done by mkin (default in -#' [saem]), then the parameters found by saem are on the transformed scale, and -#' scaling is simply done by subtracting the median. If parameter -#' transformations were done by saemix, scaling is done by dividing by the -#' median, as in the paper by Duchesne et al. (2021). +#' Produces a boxplot with all parameters from the multiple runs, divided by +#' using their medians as in the paper by Duchesne et al. (2021). #' #' @param object The [multistart] object #' @param \dots Passed to [boxplot] @@ -18,28 +14,50 @@ #' @seealso [multistart] #' @importFrom stats median #' @export -parhist <- function(object, lpos = "topleft", main = "", ...) { +parhist <- function(object, lpos = "bottomleft", main = "", ...) { + oldpar <- par(no.readonly = TRUE) + on.exit(par(oldpar, no.readonly = TRUE)) + orig <- attr(object, "orig") orig_parms <- parms(orig) start_parms <- orig$mean_dp_start all_parms <- parms(object) + + par(las = 1) + if (orig$transformations == "mkin") { + degparm_names_transformed <- names(start_parms) + degparm_index <- which(names(orig_parms) %in% degparm_names_transformed) + orig_parms[degparm_names_transformed] <- backtransform_odeparms( + orig_parms[degparm_names_transformed], + orig$mmkin$mkinmod, + transform_rates = orig$mmkin[[1]]$transform_rates, + transform_fractions = orig$mmkin[[1]]$transform_fractions) + start_parms <- backtransform_odeparms(start_parms, + orig$mmkin$mkinmod, + transform_rates = orig$mmkin[[1]]$transform_rates, + transform_fractions = orig$mmkin[[1]]$transform_fractions) + degparm_names <- names(start_parms) + + names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index])) + + all_parms[, degparm_names_transformed] <- + t(apply(all_parms[, degparm_names_transformed], 1, backtransform_odeparms, + orig$mmkin$mkinmod, + transform_rates = orig$mmkin[[1]]$transform_rates, + transform_fractions = orig$mmkin[[1]]$transform_fractions)) + colnames(all_parms)[1:length(degparm_names)] <- degparm_names + } + median_parms <- apply(all_parms, 2, median) start_scaled_parms <- rep(NA_real_, length(orig_parms)) names(start_scaled_parms) <- names(orig_parms) - if (orig$transformations == "saemix") { - orig_scaled_parms <- orig_parms / median_parms - all_scaled_parms <- t(apply(all_parms, 1, function(x) x / median_parms)) - start_scaled_parms[names(start_parms)] <- - start_parms / median_parms[names(start_parms)] - boxplot(all_scaled_parms, log = "y", main = main, ...) - } else { - orig_scaled_parms <- orig_parms - median_parms - all_scaled_parms <- t(apply(all_parms, 1, function(x) x - median_parms)) - start_scaled_parms[names(start_parms)] <- - start_parms - median_parms[names(start_parms)] - boxplot(all_scaled_parms, main = main, ...) - } + orig_scaled_parms <- orig_parms / median_parms + all_scaled_parms <- t(apply(all_parms, 1, function(x) x / median_parms)) + start_scaled_parms[names(start_parms)] <- + start_parms / median_parms[names(start_parms)] + boxplot(all_scaled_parms, log = "y", main = main, , + ylab = "Normalised parameters", ...) points(orig_scaled_parms, col = 2, cex = 2) points(start_scaled_parms, col = 3, cex = 3) |