Rename parhist to parplot and make it generic

That parhist name was not the brightest idea, as it does
not show histograms.

16 files changed, 63 insertions, 51 deletions

diff --git a/NAMESPACE b/NAMESPACE
index a5691eee..4755f631 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -36,6 +36,7 @@ S3method(parms,mkinfit)
 S3method(parms,mmkin)
 S3method(parms,multistart)
 S3method(parms,saem.mmkin)
+S3method(parplot,multistart.saem.mmkin)
 S3method(plot,mixed.mmkin)
 S3method(plot,mkinfit)
 S3method(plot,mmkin)
@@ -126,8 +127,8 @@ export(nafta)
 export(nlme)
 export(nlme_data)
 export(nlme_function)
-export(parhist)
 export(parms)
+export(parplot)
 export(plot_err)
 export(plot_res)
 export(plot_sep)
diff --git a/NEWS.md b/NEWS.md
index d288c32f..c7927559 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,6 @@
 # mkin 1.1.2
 
-- 'R/multistart.R': New method for testing multiple start parameters for hierarchical model fits, with diagnostic plotting functions 'llhist' and 'parhist'.
+- 'R/multistart.R': New method for testing multiple start parameters for hierarchical model fits, with diagnostic plotting functions 'llhist' and 'parplot'.
 
 - 'R/mhmkin.R': New method for performing multiple hierarchical mkin fits in one function call, optionally in parallel.
 
@@ -8,7 +8,7 @@
 
 - 'R/saem.R': 'logLik' and 'update' methods for 'saem.mmkin' objects.
 
-- 'R/convergence.R': New generic to show convergence information with methods for 'mmkin' and 'mhmkin' objects.
+- 'R/status.R': New generic to show status information for fit array objects with methods for 'mmkin' and 'mhmkin' objects.
 
 - 'R/illparms.R': New generic to show ill-defined parameters with methods for 'mkinfit', 'mmkin', 'saem.mmkin' and 'mhmkin' objects.
 
diff --git a/R/llhist.R b/R/llhist.R
index 22e3aa08..e158495d 100644
--- a/R/llhist.R
+++ b/R/llhist.R
@@ -25,6 +25,7 @@ llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "",
     stop("llhist is only implemented for multistart.saem.mmkin objects")
   }
 
+  ll_orig <- logLik(attr(object, "orig"))
   ll <- stats::na.omit(sapply(object, llfunc))
 
   par(las = 1)
@@ -34,7 +35,7 @@ llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "",
 
   freq_factor <- h$counts[1] / h$density[1]
 
-  abline(v = logLik(attr(object, "orig")), col = 2)
+  abline(v = ll_orig, col = 2)
 
   legend(lpos, inset = c(0.05, 0.05), bty = "n",
     lty = 1, col = c(2),
diff --git a/R/multistart.R b/R/multistart.R
index 29ccdc44..14683c11 100644
--- a/R/multistart.R
+++ b/R/multistart.R
@@ -17,7 +17,7 @@
 #' @param x The multistart object to print
 #' @return A list of [saem.mmkin] objects, with class attributes
 #' 'multistart.saem.mmkin' and 'multistart'.
-#' @seealso [parhist], [llhist]
+#' @seealso [parplot], [llhist]
 #'
 #' @references Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical
 #' identifiability in the frame of nonlinear mixed effects models: the example
@@ -40,7 +40,7 @@
 #' 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, lpos = "bottomleft")
+#' parplot(f_saem_full_multi, lpos = "bottomleft")
 #' illparms(f_saem_full)
 #'
 #' f_saem_reduced <- update(f_saem_full, no_random_effect = "log_k2")
@@ -52,7 +52,7 @@
 #' cl <- makePSOCKcluster(12)
 #' clusterExport(cl, "f_mmkin")
 #' f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cluster = cl)
-#' parhist(f_saem_reduced_multi, lpos = "topright")
+#' parplot(f_saem_reduced_multi, lpos = "topright")
 #' }
 multistart <- function(object, n = 50,
   cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(),
diff --git a/R/parhist.R b/R/parplot.R
index 0f9c9964..627a4eb8 100644
--- a/R/parhist.R
+++ b/R/parplot.R
@@ -1,4 +1,4 @@
-#' Plot parameter distributions from multistart objects
+#' Plot parameter variability of multistart objects
 #'
 #' Produces a boxplot with all parameters from the multiple runs, scaled
 #' either by the parameters of the run with the highest likelihood,
@@ -18,7 +18,13 @@
 #' @seealso [multistart]
 #' @importFrom stats median
 #' @export
-parhist <- function(object, llmin = -Inf, scale = c("best", "median"),
+parplot <- function(object, ...) {
+  UseMethod("parplot")
+}
+
+#' @rdname parplot
+#' @export
+parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best", "median"),
   lpos = "bottomleft", main = "", ...)
 {
   oldpar <- par(no.readonly = TRUE)
@@ -35,7 +41,7 @@ parhist <- function(object, llmin = -Inf, scale = c("best", "median"),
       else return(logLik(object$so))
     }
   } else {
-    stop("parhist is only implemented for multistart.saem.mmkin objects")
+    stop("parplot is only implemented for multistart.saem.mmkin objects")
   }
   ll <- sapply(object, llfunc)
   selected <- which(ll > llmin)
diff --git a/R/saem.R b/R/saem.R
index cf67b8e1..090ed3bf 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -726,7 +726,7 @@ saemix_data <- function(object, covariates = NULL, verbose = FALSE, ...) {
 #' @param \dots Passed to [saemix::logLik.SaemixObject]
 #' @param method Passed to [saemix::logLik.SaemixObject]
 #' @export
-logLik.saem.mmkin <- function(object, ..., method = c("lin", "is", "gq")) {
+logLik.saem.mmkin <- function(object, ..., method = c("is", "lin", "gq")) {
   method <- match.arg(method)
   return(logLik(object$so, method = method))
 }
diff --git a/log/test.log b/log/test.log
index 429a2e02..0e2ca6b2 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,53 +1,53 @@
 ℹ Testing mkin
 ✔ | F W S  OK | Context
 ✔ |         5 | AIC calculation
-✔ |         5 | Analytical solutions for coupled models [3.3s]
+✔ |         5 | Analytical solutions for coupled models [3.4s]
 ✔ |         5 | Calculation of Akaike weights
 ✔ |         3 | Export dataset for reading into CAKE
 ✔ |        12 | Confidence intervals and p-values [1.0s]
-✔ |     1  12 | Dimethenamid data from 2018 [43.2s]
+✔ |     1  12 | Dimethenamid data from 2018 [63.1s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_dmta.R:98:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data
 Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
 ────────────────────────────────────────────────────────────────────────────────
-✔ |        14 | Error model fitting [11.5s]
+✔ |        14 | Error model fitting [5.5s]
 ✔ |         5 | Time step normalisation
-✔ |         4 | Calculation of FOCUS chi2 error levels [1.3s]
-✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.9s]
-✔ |         4 | Test fitting the decline of metabolites from their maximum [0.9s]
-✔ |         1 | Fitting the logistic model [0.6s]
-✔ |         7 | Batch fitting and diagnosing hierarchical kinetic models [30.9s]
-✔ |     1  12 | Nonlinear mixed-effects models [0.4s]
+✔ |         4 | Calculation of FOCUS chi2 error levels [0.6s]
+✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8s]
+✔ |         4 | Test fitting the decline of metabolites from their maximum [0.4s]
+✔ |         1 | Fitting the logistic model [0.2s]
+✔ |         7 | Batch fitting and diagnosing hierarchical kinetic models [28.6s]
+✔ |     1  12 | Nonlinear mixed-effects models [0.6s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_mixed.R:74:3): saemix results are reproducible for biphasic fits
 Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ────────────────────────────────────────────────────────────────────────────────
 ✔ |         3 | Test dataset classes mkinds and mkindsg
-✔ |        10 | Special cases of mkinfit calls [0.7s]
-✔ |         3 | mkinfit features [1.0s]
-✔ |         8 | mkinmod model generation and printing [0.3s]
-✔ |         3 | Model predictions with mkinpredict [0.3s]
-✔ |         7 | Multistart method for saem.mmkin models [39.8s]
-✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.6s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [8.8s]
-✔ |        16 | Plotting [10.1s]
+✔ |        10 | Special cases of mkinfit calls [1.1s]
+✔ |         3 | mkinfit features [1.7s]
+✔ |         8 | mkinmod model generation and printing [0.4s]
+✔ |         3 | Model predictions with mkinpredict [0.6s]
+✔ |         7 | Multistart method for saem.mmkin models [93.9s]
+✔ |        16 | Evaluations according to 2015 NAFTA guidance [7.0s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [13.0s]
+✔ |        16 | Plotting [11.3s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |     1  36 | saemix parent models [72.3s]
+✔ |     1  36 | saemix parent models [72.6s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_saemix_parent.R:152:3): We can also use mkin solution methods for saem
 Reason: This still takes almost 2.5 minutes although we do not solve ODEs
 ────────────────────────────────────────────────────────────────────────────────
 ✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
 ✔ |        11 | Processing of residue series
-✔ |         7 | Fitting the SFORB model [3.8s]
+✔ |         7 | Fitting the SFORB model [3.7s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         5 | Summary [0.2s]
 ✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
-✔ |         9 | Hypothesis tests [8.0s]
+✔ |         9 | Hypothesis tests [8.1s]
 ✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 249.3 s
+Duration: 324.3 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
 • Fitting this ODE model with saemix takes about 15 minutes on my system (1)
diff --git a/man/logLik.saem.mmkin.Rd b/man/logLik.saem.mmkin.Rd
index 603f4607..bd0bb72e 100644
--- a/man/logLik.saem.mmkin.Rd
+++ b/man/logLik.saem.mmkin.Rd
@@ -4,7 +4,7 @@
 \alias{logLik.saem.mmkin}
 \title{logLik method for saem.mmkin objects}
 \usage{
-\method{logLik}{saem.mmkin}(object, ..., method = c("lin", "is", "gq"))
+\method{logLik}{saem.mmkin}(object, ..., method = c("is", "lin", "gq"))
 }
 \arguments{
 \item{object}{The fitted \link{saem.mmkin} object}
diff --git a/man/multistart.Rd b/man/multistart.Rd
index ebcc7d80..1f6773fe 100644
--- a/man/multistart.Rd
+++ b/man/multistart.Rd
@@ -77,7 +77,7 @@ 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, lpos = "bottomleft")
+parplot(f_saem_full_multi, lpos = "bottomleft")
 illparms(f_saem_full)
 
 f_saem_reduced <- update(f_saem_full, no_random_effect = "log_k2")
@@ -89,7 +89,7 @@ library(parallel)
 cl <- makePSOCKcluster(12)
 clusterExport(cl, "f_mmkin")
 f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cluster = cl)
-parhist(f_saem_reduced_multi, lpos = "topright")
+parplot(f_saem_reduced_multi, lpos = "topright")
 }
 }
 \references{
@@ -99,5 +99,5 @@ of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478.
 doi: 10.1186/s12859-021-04373-4.
 }
 \seealso{
-\link{parhist}, \link{llhist}
+\link{parplot}, \link{llhist}
 }
diff --git a/man/parhist.Rd b/man/parplot.Rd
index f86ff734..37c5841d 100644
--- a/man/parhist.Rd
+++ b/man/parplot.Rd
@@ -1,10 +1,13 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/parhist.R
-\name{parhist}
-\alias{parhist}
-\title{Plot parameter distributions from multistart objects}
+% Please edit documentation in R/parplot.R
+\name{parplot}
+\alias{parplot}
+\alias{parplot.multistart.saem.mmkin}
+\title{Plot parameter variability of multistart objects}
 \usage{
-parhist(
+parplot(object, ...)
+
+\method{parplot}{multistart.saem.mmkin}(
   object,
   llmin = -Inf,
   scale = c("best", "median"),
@@ -16,6 +19,8 @@ parhist(
 \arguments{
 \item{object}{The \link{multistart} object}
 
+\item{\dots}{Passed to \link{boxplot}}
+
 \item{llmin}{The minimum likelihood of objects to be shown}
 
 \item{scale}{By default, scale parameters using the best available fit.
@@ -24,8 +29,6 @@ If 'median', parameters are scaled using the median parameters from all fits.}
 \item{lpos}{Positioning of the legend.}
 
 \item{main}{Title of the plot}
-
-\item{\dots}{Passed to \link{boxplot}}
 }
 \description{
 Produces a boxplot with all parameters from the multiple runs, scaled
diff --git a/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg
index fa92c92a..fe38865d 100644
--- a/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg
+++ b/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg
@@ -50,6 +50,7 @@
 <rect x='228.40' y='75.47' width='146.00' height='410.67' style='stroke-width: 0.75; fill: #D3D3D3;' />
 <rect x='374.40' y='212.36' width='146.00' height='273.78' style='stroke-width: 0.75; fill: #D3D3D3;' />
 <rect x='520.40' y='212.36' width='146.00' height='273.78' style='stroke-width: 0.75; fill: #D3D3D3;' />
+<line x1='352.00' y1='502.56' x2='352.00' y2='59.04' style='stroke-width: 0.75; stroke: #DF536B;' />
 <line x1='101.38' y1='95.62' x2='122.98' y2='95.62' style='stroke-width: 0.75; stroke: #DF536B;' />
 <text x='133.78' y='99.75' style='font-size: 12.00px; font-family: sans;' textLength='51.35px' lengthAdjust='spacingAndGlyphs'>original fit</text>
 </g>
diff --git a/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg
index 4342df9b..98513d06 100644
--- a/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg
+++ b/tests/testthat/_snaps/multistart/llhist-for-sfo-fit.svg
@@ -50,7 +50,7 @@
 <rect x='228.40' y='75.47' width='146.00' height='410.67' style='stroke-width: 0.75; fill: #D3D3D3;' />
 <rect x='374.40' y='75.47' width='146.00' height='410.67' style='stroke-width: 0.75; fill: #D3D3D3;' />
 <rect x='520.40' y='349.24' width='146.00' height='136.89' style='stroke-width: 0.75; fill: #D3D3D3;' />
-<line x1='2159.04' y1='502.56' x2='2159.04' y2='59.04' style='stroke-width: 0.75; stroke: #DF536B;' />
+<line x1='232.06' y1='502.56' x2='232.06' y2='59.04' style='stroke-width: 0.75; stroke: #DF536B;' />
 <line x1='101.38' y1='95.62' x2='122.98' y2='95.62' style='stroke-width: 0.75; stroke: #DF536B;' />
 <text x='133.78' y='99.75' style='font-size: 12.00px; font-family: sans;' textLength='51.35px' lengthAdjust='spacingAndGlyphs'>original fit</text>
 </g>
diff --git a/tests/testthat/_snaps/multistart/parhist-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg
index 75519b07..75519b07 100644
--- a/tests/testthat/_snaps/multistart/parhist-for-biphasic-saemix-fit.svg
+++ b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg
diff --git a/tests/testthat/_snaps/multistart/parhist-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
index f3373901..f3373901 100644
--- a/tests/testthat/_snaps/multistart/parhist-for-sfo-fit.svg
+++ b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
diff --git a/tests/testthat/test_multistart.R b/tests/testthat/test_multistart.R
index e4885191..56eb140c 100644
--- a/tests/testthat/test_multistart.R
+++ b/tests/testthat/test_multistart.R
@@ -25,14 +25,14 @@ test_that("multistart works for saem.mmkin models", {
   skip_on_travis() # Plots are platform dependent
 
   llhist_sfo <- function() llhist(saem_sfo_s_multi)
-  parhist_sfo <- function() parhist(saem_sfo_s_multi, ylim = c(0.5, 2))
+  parplot_sfo <- function() parplot(saem_sfo_s_multi, ylim = c(0.5, 2))
   vdiffr::expect_doppelganger("llhist for sfo fit", llhist_sfo)
-  vdiffr::expect_doppelganger("parhist for sfo fit", parhist_sfo)
+  vdiffr::expect_doppelganger("parplot for sfo fit", parplot_sfo)
 
   llhist_biphasic <- function() llhist(saem_biphasic_m_multi)
-  parhist_biphasic <- function() parhist(saem_biphasic_m_multi,
+  parplot_biphasic <- function() parplot(saem_biphasic_m_multi,
     ylim = c(0.5, 2))
 
   vdiffr::expect_doppelganger("llhist for biphasic saemix fit", llhist_biphasic)
-  vdiffr::expect_doppelganger("parhist for biphasic saemix fit", parhist_biphasic)
+  vdiffr::expect_doppelganger("parplot for biphasic saemix fit", parplot_biphasic)
 })
diff --git a/vignettes/web_only/multistart.rmd b/vignettes/web_only/multistart.rmd
index 94f4ef98..27a8a96a 100644
--- a/vignettes/web_only/multistart.rmd
+++ b/vignettes/web_only/multistart.rmd
@@ -40,7 +40,7 @@ with different starting values.
 
 ```{r}
 f_saem_full_multi <- multistart(f_saem_full, n = 16, cores = 16)
-parhist(f_saem_full_multi)
+parplot(f_saem_full_multi)
 ```
 
 This confirms that the variance of k2 is the most problematic parameter, so we
@@ -51,7 +51,7 @@ for k2.
 f_saem_reduced <- update(f_saem_full, no_random_effect = "log_k2")
 illparms(f_saem_reduced)
 f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cores = 16)
-parhist(f_saem_reduced_multi, lpos = "topright")
+parplot(f_saem_reduced_multi, lpos = "topright")
 ```
 
 The results confirm that all remaining parameters can be determined with sufficient
@@ -67,7 +67,7 @@ The parameter histograms can be further improved by excluding the result with
 the low likelihood.
 
 ```{r}
-parhist(f_saem_reduced_multi, lpos = "topright", llmin = -326, ylim = c(0.5, 2))
+parplot(f_saem_reduced_multi, lpos = "topright", llmin = -326, ylim = c(0.5, 2))
 ```
 
 We can use the `anova` method to compare the models, including a likelihood ratio