aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2022-09-19 12:55:00 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2022-09-19 12:55:00 +0200
commit7a3f2ee22419608a8a634fd4d71d7176303b2f41 (patch)
tree292447a0f0f647253ac269074ba7b6ce4139c576 /R
parent062a1be18677a4e83a64d23622b5214258d9bc7d (diff)
Improve parhist and llhist
In particular, adapt the display of parameter boxplots for saem fits using mkin transformations to the way used for saem fits using saemix transformations, i.e. always show parameters on the natural scale, and normalised them by dividing by the median from the multiple runs.
Diffstat (limited to 'R')
-rw-r--r--R/llhist.R7
-rw-r--r--R/multistart.R4
-rw-r--r--R/parhist.R58
3 files changed, 46 insertions, 23 deletions
diff --git a/R/llhist.R b/R/llhist.R
index 6a6d93cd..b7f6de21 100644
--- a/R/llhist.R
+++ b/R/llhist.R
@@ -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)

Contact - Imprint