Diagnostic plots for multistart method

1 files changed, 55 insertions, 4 deletions

diff --git a/R/multistart.R b/R/multistart.R
index a3afa08b..94292e82 100644
--- a/R/multistart.R
+++ b/R/multistart.R
@@ -14,7 +14,12 @@
 #' @param n How many different combinations of starting parameters should be
 #' used?
 #' @param cores How many fits should be run in parallel?
-#' @param \dots Passed to the update function.
+#' @param \dots Passed to the update function, or to the basic plotting
+#' function in the case of the graphical function.
+#' @param x The multistart object to print
+#' @param breaks Passed to [hist]
+#' @param main title of the plot
+#' @param lpos Positioning of the legend.
 #' @return A list of [saem.mmkin] objects, with class attributes
 #' 'multistart.saem.mmkin' and 'multistart'.
 #'
@@ -39,6 +44,7 @@ multistart.saem.mmkin <- function(object, n = 50, cores = 1, ...) {
   res <- parallel::mclapply(1:n, function(x) {
     update(object, degparms_start = start_parms[x, ], ...)
   }, mc.cores = cores)
+  attr(res, "orig") <- object
   attr(res, "start_parms") <- start_parms
   class(res) <- c("multistart.saem.mmkin", "multistart")
   return(res)
@@ -53,8 +59,53 @@ print.multistart <- function(x, ...) {
 
 #' @rdname multistart
 #' @export
-summary.multistart.saem.mmkin <- function(object) {
+parms.multistart <- function(object, ...) {
+  t(sapply(object, parms))
+}
+
+#' @rdname multistart
+#' @importFrom stats median
+#' @export
+parhist <- function(object, lpos = "topleft", ...) {
+  orig <- attr(object, "orig")
+  orig_parms <- parms(orig)
+  start_parms <- orig$mean_dp_start
+  all_parms <- parms(object)
+  median_parms <- apply(all_parms, 2, median)
+  all_scaled_parms <- t(apply(all_parms, 1, function(x) x / median_parms))
+  orig_scaled_parms <- orig_parms / median_parms
+  start_scaled_parms <- rep(NA_real_, length(orig_parms))
+  names(start_scaled_parms) <- names(orig_parms)
+  start_scaled_parms[names(start_parms)] <-
+    start_parms / median_parms[names(start_parms)]
+
+  boxplot(all_scaled_parms, log = "y", ...)
+  points(orig_scaled_parms, col = 2, cex = 2)
+  points(start_scaled_parms, col = 3, cex = 3)
+  legend(lpos, inset = c(0.05, 0.05), bty = "n",
+    pch = 1, col = 3:1, lty = c(NA, NA, 1),
+    legend = c(
+      "Starting parameters",
+      "Converged parameters",
+      "Multistart runs"))
+}
+
+#' @rdname multistart
+#' @importFrom KernSmooth bkde
+#' @export
+llhist <- function(object, breaks = "Sturges", main = "", lpos = "topleft", ...) {
+  ll <- sapply(object, logLik)
+  kde <- KernSmooth::bkde(ll)
+  h <- hist(ll, freq = TRUE,
+    xlim = range(kde$x),
+    xlab = "", main = main,
+    ylab = "Frequency of log likelihoods", breaks = breaks, ...)
 
-  parm_matrix <- sapply(object, parms)
-  parm_matrix
+  freq_factor <- h$counts[1] / h$density[1]
+  lines(kde$x, freq_factor * kde$y)
+  abline(v = logLik(attr(object, "orig")), col = 2)
+  legend(lpos, inset = c(0.05, 0.05), bty = "n",
+    lty = 1, col = c(2, 1),
+    legend = c("original log likelihood",
+      "kernel density estimate"))
 }