Fix selecting by log likelihood in parhist

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