aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2022-10-17 17:28:40 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2022-10-17 17:28:40 +0200
commit73e12103447295d3d9110249b30edb222aed99fe (patch)
tree6669dbace2a317e04eeb83df51d70718b8421b02
parentb848fb360aa865c37298ee7526344b5280c700cc (diff)
Fix selecting by log likelihood in parhist
-rw-r--r--R/parhist.R29
1 files changed, 21 insertions, 8 deletions
diff --git a/R/parhist.R b/R/parhist.R
index 5d498664..0f9c9964 100644
--- a/R/parhist.R
+++ b/R/parhist.R
@@ -5,6 +5,7 @@
#' or by their medians as proposed in the paper by Duchesne et al. (2021).
#'
#' @param object The [multistart] object
+#' @param llmin The minimum likelihood of objects to be shown
#' @param scale By default, scale parameters using the best available fit.
#' If 'median', parameters are scaled using the median parameters from all fits.
#' @param main Title of the plot
@@ -17,7 +18,7 @@
#' @seealso [multistart]
#' @importFrom stats median
#' @export
-parhist <- function(object, scale = c("best", "median"),
+parhist <- function(object, llmin = -Inf, scale = c("best", "median"),
lpos = "bottomleft", main = "", ...)
{
oldpar <- par(no.readonly = TRUE)
@@ -28,6 +29,18 @@ parhist <- function(object, scale = c("best", "median"),
start_parms <- orig$mean_dp_start
all_parms <- parms(object)
+ if (inherits(object, "multistart.saem.mmkin")) {
+ llfunc <- function(object) {
+ if (inherits(object$so, "try-error")) return(NA)
+ else return(logLik(object$so))
+ }
+ } else {
+ stop("parhist is only implemented for multistart.saem.mmkin objects")
+ }
+ ll <- sapply(object, llfunc)
+ selected <- which(ll > llmin)
+ selected_parms <- all_parms[selected, ]
+
par(las = 1)
if (orig$transformations == "mkin") {
degparm_names_transformed <- names(start_parms)
@@ -45,23 +58,23 @@ parhist <- function(object, scale = c("best", "median"),
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,
+ selected_parms[, degparm_names_transformed] <-
+ t(apply(selected_parms[, degparm_names_transformed], 1, backtransform_odeparms,
orig$mmkin[[1]]$mkinmod,
transform_rates = orig$mmkin[[1]]$transform_rates,
transform_fractions = orig$mmkin[[1]]$transform_fractions))
- colnames(all_parms)[1:length(degparm_names)] <- degparm_names
+ colnames(selected_parms)[1:length(degparm_names)] <- degparm_names
}
scale <- match.arg(scale)
parm_scale <- switch(scale,
- best = all_parms[which.best(object), ],
- median = apply(all_parms, 2, median)
+ best = selected_parms[which.best(object[selected]), ],
+ median = apply(selected_parms, 2, median)
)
# Boxplots of all scaled parameters
- all_scaled_parms <- t(apply(all_parms, 1, function(x) x / parm_scale))
- boxplot(all_scaled_parms, log = "y", main = main, ,
+ selected_scaled_parms <- t(apply(selected_parms, 1, function(x) x / parm_scale))
+ boxplot(selected_scaled_parms, log = "y", main = main, ,
ylab = "Normalised parameters", ...)
# Show starting parameters

Contact - Imprint