Fix handling of multistart fits with failures

13 files changed, 159 insertions, 104 deletions

diff --git a/NAMESPACE b/NAMESPACE
index bd12d0db..388a16fd 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -13,6 +13,8 @@ S3method(aw,multistart)
 S3method(confint,mkinfit)
 S3method(convergence,mhmkin)
 S3method(convergence,mmkin)
+S3method(convergence,multistart)
+S3method(convergence,multistart.saem.mmkin)
 S3method(f_time_norm_focus,mkindsg)
 S3method(f_time_norm_focus,numeric)
 S3method(illparms,mhmkin)
@@ -42,6 +44,7 @@ S3method(plot,mmkin)
 S3method(plot,nafta)
 S3method(print,convergence.mhmkin)
 S3method(print,convergence.mmkin)
+S3method(print,convergence.multistart)
 S3method(print,illparms.mhmkin)
 S3method(print,illparms.mmkin)
 S3method(print,mhmkin)
diff --git a/R/llhist.R b/R/llhist.R
index b7f6de21..9ddf5b10 100644
--- a/R/llhist.R
+++ b/R/llhist.R
@@ -16,7 +16,16 @@ llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "", ...)
   oldpar <- par(no.readonly = TRUE)
   on.exit(par(oldpar, no.readonly = TRUE))
 
-  ll <- sapply(object, logLik)
+  if (inherits(object, "multistart.saem.mmkin")) {
+    llfunc <- function(object) {
+      if (inherits(object$so, "try-error")) return(NA)
+      else return(logLik(object$so))
+    }
+  } else {
+    stop("llhist is only implemented for multistart.saem.mmkin objects")
+  }
+
+  ll <- stats::na.omit(sapply(object, llfunc))
   kde <- KernSmooth::bkde(ll)
 
   par(las = 1)
diff --git a/R/mhmkin.R b/R/mhmkin.R
index b72ae318..a1475ef9 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -179,7 +179,7 @@ AIC.mhmkin <- function(object, ..., k = 2) {
 BIC.mhmkin <- function(object, ...) {
   res <- sapply(object, function(x) {
     if (inherits(x, "try-error")) return(NA)
-    else return(BIC(x$so, k = k))
+    else return(BIC(x$so))
   })
   dim(res) <- dim(object)
   dimnames(res) <- dimnames(object)
diff --git a/R/multistart.R b/R/multistart.R
index cc55feae..b65c0bee 100644
--- a/R/multistart.R
+++ b/R/multistart.R
@@ -92,15 +92,51 @@ multistart.saem.mmkin <- function(object, n = 50, cores = 1,
   return(res)
 }
 
-#' @rdname multistart
 #' @export
-print.multistart <- function(x, ...) {
-  cat("Multistart object with", length(x), "fits of the following type:\n\n")
-  print(x[[1]])
+convergence.multistart <- function(object, ...) {
+  all_summary_warnings <- character()
+
+  result <- lapply(object,
+    function(fit) {
+      if (inherits(fit, "try-error")) return("E")
+      else {
+        return("OK")
+      }
+  })
+  result <- unlist(result)
+
+  class(result) <- "convergence.multistart"
+  return(result)
+}
+
+#' @export
+convergence.multistart.saem.mmkin <- function(object, ...) {
+  all_summary_warnings <- character()
+
+  result <- lapply(object,
+    function(fit) {
+      if (inherits(fit$so, "try-error")) return("E")
+      else {
+        return("OK")
+      }
+  })
+  result <- unlist(result)
+
+  class(result) <- "convergence.multistart"
+  return(result)
+}
+
+#' @export
+print.convergence.multistart <- function(x, ...) {
+  class(x) <- NULL
+  print(table(x, dnn = NULL))
+  if (any(x == "OK")) cat("OK: Fit terminated successfully\n")
+  if (any(x == "E")) cat("E: Error\n")
 }
 
 #' @rdname multistart
 #' @export
-parms.multistart <- function(object, ...) {
-  t(sapply(object, parms))
+print.multistart <- function(x, ...) {
+  cat("<multistart> object with", length(x), "fits:\n")
+  print(convergence(x))
 }
diff --git a/R/parhist.R b/R/parhist.R
index 69aafe02..10730873 100644
--- a/R/parhist.R
+++ b/R/parhist.R
@@ -29,20 +29,20 @@ parhist <- function(object, lpos = "bottomleft", main = "", ...) {
     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,
+      orig$mmkin[[1]]$mkinmod,
       transform_rates = orig$mmkin[[1]]$transform_rates,
       transform_fractions = orig$mmkin[[1]]$transform_fractions)
     start_parms <- backtransform_odeparms(start_parms,
-      orig$mmkin$mkinmod,
+      orig$mmkin[[1]]$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,
+      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
diff --git a/R/parms.mkinfit.R b/R/parms.R
index 83766355..bd4e479b 100644
--- a/R/parms.mkinfit.R
+++ b/R/parms.R
@@ -67,3 +67,16 @@ parms.mmkin <- function(object, transformed = FALSE, errparms = TRUE, ...)
   }
   return(res)
 }
+
+#' @param exclude_failed For [multistart] objects, should rows for failed fits
+#' be removed from the returned parameter matrix?
+#' @rdname parms
+#' @export
+parms.multistart <- function(object, exclude_failed = TRUE, ...) {
+  res <- t(sapply(object, parms))
+  successful <- which(!is.na(res[, 1]))
+  first_success <- successful[1]
+  colnames(res) <- names(parms(object[[first_success]]))
+  if (exclude_failed) res <- res[successful, ]
+  return(res)
+}
diff --git a/R/saem.R b/R/saem.R
index 0d0d9b8a..8929cf6f 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -137,9 +137,16 @@ saem.mmkin <- function(object,
     transformations = transformations, ...)
   d_saemix <- saemix_data(object, verbose = verbose)
 
+  fit_failed <- FALSE
+  FIM_failed <- NULL
   fit_time <- system.time({
-    utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet)
-    FIM_failed <- NULL
+    utils::capture.output(f_saemix <- try(saemix::saemix(m_saemix, d_saemix, control)), split = !quiet)
+    if (inherits(f_saemix, "try-error")) fit_failed <- TRUE
+  })
+
+  return_data <- nlme_data(object)
+
+  if (!fit_failed) {
     if (any(is.na(f_saemix@results@se.fixed))) FIM_failed <- c(FIM_failed, "fixed effects")
     if (any(is.na(c(f_saemix@results@se.omega, f_saemix@results@se.respar)))) {
       FIM_failed <- c(FIM_failed, "random effects and residual error parameters")
@@ -147,38 +154,37 @@ saem.mmkin <- function(object,
     if (!is.null(FIM_failed) & fail_with_errors) {
       stop("Could not invert FIM for ", paste(FIM_failed, collapse = " and "))
     }
-  })
 
-  transparms_optim <- f_saemix@results@fixed.effects
-  names(transparms_optim) <- f_saemix@results@name.fixed
+    transparms_optim <- f_saemix@results@fixed.effects
+    names(transparms_optim) <- f_saemix@results@name.fixed
 
-  if (transformations == "mkin") {
-    bparms_optim <- backtransform_odeparms(transparms_optim,
-      object[[1]]$mkinmod,
-      object[[1]]$transform_rates,
-      object[[1]]$transform_fractions)
-  } else {
-    bparms_optim <- transparms_optim
-  }
+    if (transformations == "mkin") {
+      bparms_optim <- backtransform_odeparms(transparms_optim,
+        object[[1]]$mkinmod,
+        object[[1]]$transform_rates,
+        object[[1]]$transform_fractions)
+    } else {
+      bparms_optim <- transparms_optim
+    }
 
-  return_data <- nlme_data(object)
-  saemix_data_ds <- f_saemix@data@data$ds
-  mkin_ds_order <- as.character(unique(return_data$ds))
-  saemix_ds_order <- unique(saemix_data_ds)
-
-  psi <- saemix::psi(f_saemix)
-  rownames(psi) <- saemix_ds_order
-  return_data$predicted <- f_saemix@model@model(
-    psi = psi[mkin_ds_order, ],
-    id = as.numeric(return_data$ds),
-    xidep = return_data[c("time", "name")])
-
-  return_data <- transform(return_data,
-    residual = value - predicted,
-    std = sigma_twocomp(predicted,
-      f_saemix@results@respar[1], f_saemix@results@respar[2]))
-  return_data <- transform(return_data,
-    standardized = residual / std)
+    saemix_data_ds <- f_saemix@data@data$ds
+    mkin_ds_order <- as.character(unique(return_data$ds))
+    saemix_ds_order <- unique(saemix_data_ds)
+
+    psi <- saemix::psi(f_saemix)
+    rownames(psi) <- saemix_ds_order
+    return_data$predicted <- f_saemix@model@model(
+      psi = psi[mkin_ds_order, ],
+      id = as.numeric(return_data$ds),
+      xidep = return_data[c("time", "name")])
+
+    return_data <- transform(return_data,
+      residual = value - predicted,
+      std = sigma_twocomp(predicted,
+        f_saemix@results@respar[1], f_saemix@results@respar[2]))
+    return_data <- transform(return_data,
+      standardized = residual / std)
+  }
 
   result <- list(
     mkinmod = object[[1]]$mkinmod,
@@ -187,15 +193,13 @@ saem.mmkin <- function(object,
     transformations = transformations,
     transform_rates = object[[1]]$transform_rates,
     transform_fractions = object[[1]]$transform_fractions,
+    sm = m_saemix,
     so = f_saemix,
     call = call,
     time = fit_time,
     mean_dp_start = attr(m_saemix, "mean_dp_start"),
-    bparms.optim = bparms_optim,
     bparms.fixed = object[[1]]$bparms.fixed,
     data = return_data,
-    mkin_ds_order = mkin_ds_order,
-    saemix_ds_order = saemix_ds_order,
     err_mod = object[[1]]$err_mod,
     date.fit = date(),
     saemixversion = as.character(utils::packageVersion("saemix")),
@@ -203,6 +207,12 @@ saem.mmkin <- function(object,
     Rversion = paste(R.version$major, R.version$minor, sep=".")
   )
 
+  if (!fit_failed) {
+    result$mkin_ds_order <- mkin_ds_order
+    result$saemix_ds_order <- saemix_ds_order
+    result$bparms.optim <- bparms_optim
+  }
+
   class(result) <- c("saem.mmkin", "mixed.mmkin")
   return(result)
 }
@@ -222,16 +232,20 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
     length(unique(x$data$name)), "variable(s) grouped in",
     length(unique(x$data$ds)), "datasets\n")
 
-  cat("\nLikelihood computed by importance sampling\n")
-  print(data.frame(
-      AIC = AIC(x$so, type = "is"),
-      BIC = BIC(x$so, type = "is"),
-      logLik = logLik(x$so, type = "is"),
-      row.names = " "), digits = digits)
-
-  cat("\nFitted parameters:\n")
-  conf.int <- parms(x, ci = TRUE)
-  print(conf.int, digits = digits)
+  if (inherits(x$so, "try-error")) {
+    cat("\nFit did not terminate successfully\n")
+  } else {
+    cat("\nLikelihood computed by importance sampling\n")
+    print(data.frame(
+        AIC = AIC(x$so, type = "is"),
+        BIC = BIC(x$so, type = "is"),
+        logLik = logLik(x$so, type = "is"),
+        row.names = " "), digits = digits)
+
+    cat("\nFitted parameters:\n")
+    conf.int <- parms(x, ci = TRUE)
+    print(conf.int, digits = digits)
+  }
 
   invisible(x)
 }
@@ -618,12 +632,27 @@ update.saem.mmkin <- function(object, ..., evaluate = TRUE) {
 #' @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(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, ]
+  cov.mod <- object$sm@covariance.model
+  n_cov_mod_parms <- sum(cov.mod[upper.tri(cov.mod, diag = TRUE)])
+  n_parms <- length(object$sm@name.modpar) +
+    n_cov_mod_parms +
+    length(object$sm@name.sigma)
+
+  if (inherits(object$so, "try-error")) {
+    conf.int <- matrix(rep(NA, 3 * n_parms), ncol = 3)
+    colnames(conf.int) <- c("estimate", "lower", "upper")
+  } else {
+    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, ]
+    conf.int.cov <- grepl("^Cov\\.", rownames(conf.int))
+    conf.int <- conf.int[!conf.int.cov, ]
+  }
   estimate <- conf.int[, "estimate"]
+
   names(estimate) <- rownames(conf.int)
+
   if (ci) return(conf.int)
   else return(estimate)
 }
diff --git a/log/check.log b/log/check.log
index 60298c09..8fce3fd1 100644
--- a/log/check.log
+++ b/log/check.log
@@ -48,15 +48,7 @@ Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’
 * checking Rd cross-references ... OK
 * checking for missing documentation entries ... OK
 * checking for code/documentation mismatches ... OK
-* checking Rd \usage sections ... WARNING
-Undocumented arguments in documentation object 'multistart'
-  ‘x’
-
-Functions with \usage entries need to have the appropriate \alias
-entries, and all their arguments documented.
-The \usage entries must correspond to syntactically valid R code.
-See chapter ‘Writing R documentation files’ in the ‘Writing R
-Extensions’ manual.
+* checking Rd \usage sections ... OK
 * checking Rd contents ... OK
 * checking for unstated dependencies in examples ... OK
 * checking contents of ‘data’ directory ... OK
@@ -77,9 +69,5 @@ Extensions’ manual.
 * checking for detritus in the temp directory ... OK
 * DONE
 
-Status: 1 WARNING
-See
-  ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’
-for details.
-
+Status: OK
 
diff --git a/man/multistart.Rd b/man/multistart.Rd
index aa352267..78ff4614 100644
--- a/man/multistart.Rd
+++ b/man/multistart.Rd
@@ -4,7 +4,6 @@
 \alias{multistart}
 \alias{multistart.saem.mmkin}
 \alias{print.multistart}
-\alias{parms.multistart}
 \title{Perform a hierarchical model fit with multiple starting values}
 \usage{
 multistart(
@@ -18,8 +17,6 @@ multistart(
 \method{multistart}{saem.mmkin}(object, n = 50, cores = 1, cluster = NULL, ...)
 
 \method{print}{multistart}(x, ...)
-
-\method{parms}{multistart}(object, ...)
 }
 \arguments{
 \item{object}{The fit object to work with}
diff --git a/man/parms.Rd b/man/parms.Rd
index acae2d91..5c0e8895 100644
--- a/man/parms.Rd
+++ b/man/parms.Rd
@@ -1,9 +1,10 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/parms.mkinfit.R
+% Please edit documentation in R/parms.R
 \name{parms}
 \alias{parms}
 \alias{parms.mkinfit}
 \alias{parms.mmkin}
+\alias{parms.multistart}
 \title{Extract model parameters}
 \usage{
 parms(object, ...)
@@ -11,6 +12,8 @@ parms(object, ...)
 \method{parms}{mkinfit}(object, transformed = FALSE, errparms = TRUE, ...)
 
 \method{parms}{mmkin}(object, transformed = FALSE, errparms = TRUE, ...)
+
+\method{parms}{multistart}(object, exclude_failed = TRUE, ...)
 }
 \arguments{
 \item{object}{A fitted model object.}
@@ -22,6 +25,9 @@ during the optimisation?}
 
 \item{errparms}{Should the error model parameters be returned
 in addition to the degradation parameters?}
+
+\item{exclude_failed}{For \link{multistart} objects, should rows for failed fits
+be removed from the returned parameter matrix?}
 }
 \value{
 Depending on the object, a numeric vector of fitted model parameters,
diff --git a/vignettes/web_only/multistart.R b/vignettes/web_only/multistart.R
deleted file mode 100644
index b995c545..00000000
--- a/vignettes/web_only/multistart.R
+++ /dev/null
@@ -1,26 +0,0 @@
-## -----------------------------------------------------------------------------
-library(mkin)
-dmta_ds <- lapply(1:7, function(i) {
-  ds_i <- dimethenamid_2018$ds[[i]]$data
-  ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
-  ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]
-  ds_i
-})
-names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
-dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
-dmta_ds[["Elliot 1"]] <- dmta_ds[["Elliot 2"]] <- NULL
-
-## -----------------------------------------------------------------------------
-f_mmkin <- mmkin("DFOP", dmta_ds, error_model = "tc", cores = 7, quiet = TRUE)
-f_saem_full <- saem(f_mmkin)
-f_saem_full_multi <- multistart(f_saem_full, n = 16, cores = 16)
-parhist(f_saem_full_multi)
-
-## -----------------------------------------------------------------------------
-f_saem_reduced <- update(f_saem_full, covariance.model = diag(c(1, 1, 0, 1)))
-f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cores = 16)
-parhist(f_saem_reduced_multi, lpos = "topright")
-
-## -----------------------------------------------------------------------------
-llhist(f_saem_reduced_multi)
-
diff --git a/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-2-1.png b/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-2-1.png
new file mode 100644
index 00000000..17168f74
--- /dev/null
+++ b/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-2-1.png
diff --git a/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png b/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png
new file mode 100644
index 00000000..17d99e4e
--- /dev/null
+++ b/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png