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.

1 files changed, 38 insertions, 20 deletions

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)