Diagnostic plots for multistart method

3 files changed, 90 insertions, 13 deletions

diff --git a/R/aw.R b/R/aw.R
index f46b20ec..b3992f94 100644
--- a/R/aw.R
+++ b/R/aw.R
@@ -30,6 +30,14 @@
 #' @export
 aw <- function(object, ...) UseMethod("aw")
 
+.aw <- function(all_objects) {
+  AIC_all <- sapply(all_objects, AIC)
+  delta_i <- AIC_all - min(AIC_all)
+  denom <- sum(exp(-delta_i/2))
+  w_i <- exp(-delta_i/2) / denom
+  return(w_i)
+}
+
 #' @export
 #' @rdname aw
 aw.mkinfit <- function(object, ...) {
@@ -43,11 +51,7 @@ aw.mkinfit <- function(object, ...) {
     }
   }
   all_objects <- list(object, ...)
-  AIC_all <- sapply(all_objects, AIC)
-  delta_i <- AIC_all - min(AIC_all)
-  denom <- sum(exp(-delta_i/2))
-  w_i <- exp(-delta_i/2) / denom
-  return(w_i)
+  .aw(all_objects)
 }
 
 #' @export
@@ -56,3 +60,25 @@ aw.mmkin <- function(object, ...) {
   if (ncol(object) > 1) stop("Please supply an mmkin column object")
   do.call(aw, object)
 }
+
+#' @export
+#' @rdname aw
+aw.mixed.mmkin <- function(object, ...) {
+  oo <- list(...)
+  data_object <- object$data[c("ds", "name", "time", "value")]
+  for (i in seq_along(oo)) {
+    if (!inherits(oo[[i]], "mixed.mmkin")) stop("Please supply objects inheriting from mixed.mmkin")
+    data_other_object <- oo[[i]]$data[c("ds", "name", "time", "value")]
+    if (!identical(data_object, data_other_object)) {
+      stop("It seems that the mixed.mmkin objects have not all been fitted to the same data")
+    }
+  }
+  all_objects <- list(object, ...)
+  .aw(all_objects)
+}
+
+#' @export
+#' @rdname aw
+aw.multistart <- function(object, ...) {
+  do.call(aw, object)
+}
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"))
 }
diff --git a/R/saem.R b/R/saem.R
index 875456c3..0d0d9b8a 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -568,7 +568,7 @@ saemix_data <- function(object, verbose = FALSE, ...) {
 
   ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")])
   names(ds_list) <- ds_names
-  ds_saemix_all <- purrr::map_dfr(ds_list, function(x) x, .id = "ds")
+  ds_saemix_all <- vctrs::vec_rbind(!!!ds_list, .names_to = "ds")
   ds_saemix <- data.frame(ds = ds_saemix_all$ds,
     name = as.character(ds_saemix_all$variable),
     time = ds_saemix_all$time,
@@ -617,9 +617,9 @@ update.saem.mmkin <- function(object, ..., evaluate = TRUE) {
 #' @rdname saem
 #' @param ci Should a matrix with estimates and confidence interval boundaries
 #' be returned? If FALSE (default), a vector of estimates is returned.
-parms.saem.mmkin <- function(x, ci = FALSE, ...) {
-  conf.int <- x$so@results@conf.int[c("estimate", "lower", "upper")]
-  rownames(conf.int) <- x$so@results@conf.int[["name"]]
+parms.saem.mmkin <- function(object, ci = FALSE, ...) {
+  conf.int <- object$so@results@conf.int[c("estimate", "lower", "upper")]
+  rownames(conf.int) <- object$so@results@conf.int[["name"]]
   conf.int.var <- grepl("^Var\\.", rownames(conf.int))
   conf.int <- conf.int[!conf.int.var, ]
   estimate <- conf.int[, "estimate"]