Fix plotting saem fits with covariates

7 files changed, 94 insertions, 58 deletions

diff --git a/R/illparms.R b/R/illparms.R
index 4f5f8b00..c9a4f854 100644
--- a/R/illparms.R
+++ b/R/illparms.R
@@ -45,7 +45,7 @@ illparms.mkinfit <- function(object, conf.level = 0.95, ...) {
   return(ret)
 }
 
-#' @rdname
+#' @rdname illparms
 #' @export
 print.illparms.mkinfit <- function(x, ...) {
   class(x) <- NULL
@@ -111,6 +111,7 @@ illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod
   return(failed)
 }
 
+#' @rdname illparms
 #' @export
 print.illparms.saem.mmkin <- print.illparms.mkinfit
 
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index 6815bfd2..6f613b48 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -4,27 +4,30 @@ utils::globalVariables("ds")
 #'
 #' @param x An object of class [mixed.mmkin], [saem.mmkin] or [nlme.mmkin]
 #' @param i A numeric index to select datasets for which to plot the individual predictions,
-#'   in case plots get too large
+#' in case plots get too large
 #' @inheritParams plot.mkinfit
 #' @param standardized Should the residuals be standardized? Only takes effect if
-#'   `resplot = "time"`.
+#' `resplot = "time"`.
+#' @param pop_curve Per default, a population curve is drawn in case
+#' population parameters are fitted by the model, e.g. for saem objects.
+#' In case there is a covariate model, no population curve is currently shown.
 #' @param pred_over Named list of alternative predictions as obtained
-#'   from [mkinpredict] with a compatible [mkinmod].
+#' from [mkinpredict] with a compatible [mkinmod].
 #' @param test_log_parms Passed to [mean_degparms] in the case of an
-#'   [mixed.mmkin] object
+#' [mixed.mmkin] object
 #' @param conf.level Passed to [mean_degparms] in the case of an
-#'   [mixed.mmkin] object
+#' [mixed.mmkin] object
 #' @param default_log_parms Passed to [mean_degparms] in the case of an
-#'   [mixed.mmkin] object
+#' [mixed.mmkin] object
 #' @param rel.height.legend The relative height of the legend shown on top
 #' @param rel.height.bottom The relative height of the bottom plot row
 #' @param ymax Vector of maximum y axis values
 #' @param ncol.legend Number of columns to use in the legend
 #' @param nrow.legend Number of rows to use in the legend
 #' @param resplot Should the residuals plotted against time or against
-#'   predicted values?
+#' predicted values?
 #' @param col_ds Colors used for plotting the observed data and the
-#'   corresponding model prediction lines for the different datasets.
+#' corresponding model prediction lines for the different datasets.
 #' @param pch_ds Symbols to be used for plotting the data.
 #' @param lty_ds Line types to be used for the model predictions.
 #' @importFrom stats coefficients
@@ -68,6 +71,7 @@ plot.mixed.mmkin <- function(x,
   xlab = "Time",
   xlim = range(x$data$time),
   resplot = c("predicted", "time"),
+  pop_curve = "auto",
   pred_over = NULL,
   test_log_parms = FALSE,
   conf.level = 0.6,
@@ -90,8 +94,15 @@ plot.mixed.mmkin <- function(x,
   backtransform = TRUE
 
   if (identical(class(x), "mixed.mmkin")) {
-    degparms_pop <- mean_degparms(x$mmkin, test_log_parms = test_log_parms,
-      conf.level = conf.level, default_log_parms = default_log_parms)
+    if (identical(pop_curve, "auto")) {
+      pop_curve <- FALSE
+    } else {
+      pop_curve <- TRUE
+    }
+    if (pop_curve) {
+      degparms_pop <- mean_degparms(x$mmkin, test_log_parms = test_log_parms,
+        conf.level = conf.level, default_log_parms = default_log_parms)
+    }
 
     degparms_tmp <- parms(x$mmkin, transformed = TRUE)
     degparms_i <- as.data.frame(t(degparms_tmp[setdiff(rownames(degparms_tmp), names(fit_1$errparms)), ]))
@@ -100,6 +111,11 @@ plot.mixed.mmkin <- function(x,
   }
 
   if (inherits(x, "nlme.mmkin")) {
+    if (identical(pop_curve, "auto")) {
+      pop_curve <- TRUE
+    } else {
+      pop_curve <- FALSE
+    }
     degparms_i <- coefficients(x)
     degparms_pop <- nlme::fixef(x)
     residuals <- residuals(x,
@@ -111,24 +127,21 @@ plot.mixed.mmkin <- function(x,
     psi <- saemix::psi(x$so)
     rownames(psi) <- x$saemix_ds_order
     degparms_i <- psi[ds_names, ]
-    degparms_i_names <- setdiff(x$so@results@name.fixed, names(fit_1$errparms))
-    colnames(degparms_i) <- degparms_i_names
+    degparms_i_names <- colnames(degparms_i)
     residual_type = ifelse(standardized, "standardized", "residual")
     residuals <- x$data[[residual_type]]
-    degparms_pop <- x$so@results@fixed.effects
-    names(degparms_pop) <- degparms_i_names
-  }
 
-  if (inherits(x, "nlmixr.mmkin")) {
-    eta_i <- random.effects(x$nm)[-1]
-    names(eta_i) <- gsub("^eta.", "", names(eta_i))
-    degparms_i <- eta_i
-    degparms_pop <- x$nm$theta
-    for (parm_name in names(degparms_i)) {
-      degparms_i[parm_name] <- eta_i[parm_name] + degparms_pop[parm_name]
+    if (identical(pop_curve, "auto")) {
+      if (nrow(x$sm@covariate.model) == 0) {
+        pop_curve <- TRUE
+      } else {
+        pop_curve <- FALSE
+      }
+    }
+    if (pop_curve) {
+      degparms_pop <- x$so@results@fixed.effects
+      names(degparms_pop) <- degparms_i_names
     }
-    residual_type = ifelse(standardized, "standardized", "residual")
-    residuals <- x$data[[residual_type]]
   }
 
   degparms_fixed <- fit_1$fixed$value
@@ -140,8 +153,6 @@ plot.mixed.mmkin <- function(x,
   degparms_all_names <- c(names(degparms_i), names(degparms_fixed))
   colnames(degparms_all) <- degparms_all_names
 
-  degparms_all_pop <- c(degparms_pop, degparms_fixed)
-
   odeini_names <- grep("_0$", degparms_all_names, value = TRUE)
   odeparms_names <- setdiff(degparms_all_names, odeini_names)
 
@@ -175,24 +186,29 @@ plot.mixed.mmkin <- function(x,
   names(pred_list) <- ds_names[i]
   pred_ds <- vctrs::vec_rbind(!!!pred_list, .names_to = "ds")
 
-  odeparms_pop_trans <- degparms_all_pop[odeparms_names]
+  if (pop_curve) {
+    degparms_all_pop <- c(degparms_pop, degparms_fixed)
+    odeparms_pop_trans <- degparms_all_pop[odeparms_names]
 
-  if (backtransform) {
-    odeparms_pop <- backtransform_odeparms(odeparms_pop_trans,
-      x$mkinmod,
-      transform_rates = fit_1$transform_rates,
-      transform_fractions = fit_1$transform_fractions)
-  } else {
-    odeparms_pop <- odeparms_pop_trans
-  }
+    if (backtransform) {
+      odeparms_pop <- backtransform_odeparms(odeparms_pop_trans,
+        x$mkinmod,
+        transform_rates = fit_1$transform_rates,
+        transform_fractions = fit_1$transform_fractions)
+    } else {
+      odeparms_pop <- odeparms_pop_trans
+    }
 
-  odeini_pop <- degparms_all_pop[odeini_names]
-  names(odeini_pop) <- gsub("_0", "", odeini_names)
+    odeini_pop <- degparms_all_pop[odeini_names]
+    names(odeini_pop) <- gsub("_0", "", odeini_names)
 
-  pred_pop <- as.data.frame(
-    mkinpredict(x$mkinmod, odeparms_pop, odeini_pop,
-      outtimes, solution_type = solution_type,
-      atol = fit_1$atol, rtol = fit_1$rtol))
+    pred_pop <- as.data.frame(
+      mkinpredict(x$mkinmod, odeparms_pop, odeini_pop,
+        outtimes, solution_type = solution_type,
+        atol = fit_1$atol, rtol = fit_1$rtol))
+  } else {
+    pred_pop <- NULL
+  }
 
   # Start of graphical section
   oldpar <- par(no.readonly = TRUE)
@@ -217,12 +233,18 @@ plot.mixed.mmkin <- function(x,
   # Empty plot with legend
   if (!is.null(pred_over)) lty_over <- seq(2, length.out = length(pred_over))
   else lty_over <- NULL
-  n_pop <- 1 + length(lty_over)
-  lty_pop <- c(1, lty_over)
+  if (pop_curve) {
+    lty_pop <- 1
+    names(lty_pop) <- "Population"
+  } else {
+    lty_pop <- NULL
+  }
+  n_pop <- length(lty_pop) + length(lty_over)
+  lty_pop <- c(lty_pop, lty_over)
 
   plot(0, type = "n", axes = FALSE, ann = FALSE)
   legend("center", bty = "n", ncol = ncol.legend,
-    legend = c("Population", names(pred_over), ds_names[i]),
+    legend = c(names(lty_pop), names(pred_over), ds_names[i]),
     lty = c(lty_pop, lty_ds),
     lwd = c(rep(2, n_pop), rep(1, length(i))),
     col = c(rep(1, n_pop), col_ds),
@@ -271,8 +293,10 @@ plot.mixed.mmkin <- function(x,
         col = col_ds[ds_i], lty = lty_ds[ds_i])
     }
 
-    lines(pred_pop$time, pred_pop[[obs_var]],
-      type = "l", lwd = 2, lty = lty_pop)
+    if (pop_curve) {
+      lines(pred_pop$time, pred_pop[[obs_var]],
+        type = "l", lwd = 2, lty = lty_pop)
+    }
 
     if (identical(maxabs, "auto")) {
       maxabs = max(abs(observed_row$residual), na.rm = TRUE)
diff --git a/log/check.log b/log/check.log
index 4092108e..44830443 100644
--- a/log/check.log
+++ b/log/check.log
@@ -59,7 +59,7 @@ The Date field is over a month old.
 * checking data for ASCII and uncompressed saves ... OK
 * checking installed files from ‘inst/doc’ ... OK
 * checking files in ‘vignettes’ ... OK
-* checking examples ... [15s/15s] OK
+* checking examples ... [16s/16s] OK
 * checking for unstated dependencies in ‘tests’ ... OK
 * checking tests ... SKIPPED
 * checking for unstated dependencies in vignettes ... OK
diff --git a/log/test.log b/log/test.log
index ab7402b9..611b66b6 100644
--- a/log/test.log
+++ b/log/test.log
@@ -5,7 +5,7 @@
 ✔ |         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 [31.2s]
+✔ |     1  12 | Dimethenamid data from 2018 [31.9s]
 ────────────────────────────────────────────────────────────────────────────────
 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
@@ -16,7 +16,7 @@ Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
 ✔ |        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.3s]
 ✔ |         1 | Fitting the logistic model [0.2s]
-✔ |         8 | Batch fitting and diagnosing hierarchical kinetic models [14.2s]
+✔ |         8 | Batch fitting and diagnosing hierarchical kinetic models [14.5s]
 ✔ |     1  12 | Nonlinear mixed-effects models [0.3s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_mixed.R:74:3): saemix results are reproducible for biphasic fits
@@ -26,28 +26,28 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |        10 | Special cases of mkinfit calls [0.4s]
 ✔ |         3 | mkinfit features [0.7s]
 ✔ |         8 | mkinmod model generation and printing [0.2s]
-✔ |         3 | Model predictions with mkinpredict [0.3s]
-✔ |         7 | Multistart method for saem.mmkin models [35.3s]
+✔ |         3 | Model predictions with mkinpredict [0.4s]
+✔ |         7 | Multistart method for saem.mmkin models [36.2s]
 ✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.4s]
 ✔ |         9 | Nonlinear mixed-effects models with nlme [8.8s]
-✔ |        16 | Plotting [9.9s]
+✔ |        16 | Plotting [10.0s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |     1  37 | saemix parent models [70.3s]
+✔ |     1  37 | saemix parent models [71.4s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_saemix_parent.R:153: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]
+✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5s]
 ✔ |        11 | Processing of residue series
 ✔ |         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: 203.2 s
+Duration: 206.7 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
 • Fitting this ODE model with saemix takes about 15 minutes on my system (1)
diff --git a/man/illparms.Rd b/man/illparms.Rd
index 98ac7aa3..14be9c35 100644
--- a/man/illparms.Rd
+++ b/man/illparms.Rd
@@ -3,9 +3,11 @@
 \name{illparms}
 \alias{illparms}
 \alias{illparms.mkinfit}
+\alias{print.illparms.mkinfit}
 \alias{illparms.mmkin}
 \alias{print.illparms.mmkin}
 \alias{illparms.saem.mmkin}
+\alias{print.illparms.saem.mmkin}
 \alias{illparms.mhmkin}
 \alias{print.illparms.mhmkin}
 \title{Method to get the names of ill-defined parameters}
@@ -14,12 +16,16 @@ illparms(object, ...)
 
 \method{illparms}{mkinfit}(object, conf.level = 0.95, ...)
 
+\method{print}{illparms.mkinfit}(x, ...)
+
 \method{illparms}{mmkin}(object, conf.level = 0.95, ...)
 
 \method{print}{illparms.mmkin}(x, ...)
 
 \method{illparms}{saem.mmkin}(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...)
 
+\method{print}{illparms.saem.mmkin}(x, ...)
+
 \method{illparms}{mhmkin}(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...)
 
 \method{print}{illparms.mhmkin}(x, ...)
diff --git a/man/plot.mixed.mmkin.Rd b/man/plot.mixed.mmkin.Rd
index 33b4a67f..9c4474ff 100644
--- a/man/plot.mixed.mmkin.Rd
+++ b/man/plot.mixed.mmkin.Rd
@@ -12,6 +12,7 @@
   xlab = "Time",
   xlim = range(x$data$time),
   resplot = c("predicted", "time"),
+  pop_curve = "auto",
   pred_over = NULL,
   test_log_parms = FALSE,
   conf.level = 0.6,
@@ -49,6 +50,10 @@ variables in the model.}
 \item{resplot}{Should the residuals plotted against time or against
 predicted values?}
 
+\item{pop_curve}{Per default, a population curve is drawn in case
+population parameters are fitted by the model, e.g. for saem objects.
+In case there is a covariate model, no population curve is currently shown.}
+
 \item{pred_over}{Named list of alternative predictions as obtained
 from \link{mkinpredict} with a compatible \link{mkinmod}.}
 
diff --git a/tests/testthat/test_plot.R b/tests/testthat/test_plot.R
index 58a00662..13058c00 100644
--- a/tests/testthat/test_plot.R
+++ b/tests/testthat/test_plot.R
@@ -44,7 +44,7 @@ test_that("Plotting mkinfit, mmkin and mixed model objects is reproducible", {
 
   f_uba_dfop_sfo_saem <- saem(f_uba_mmkin["DFOP-SFO", ], quiet = TRUE, transformations = "saemix")
 
-  plot_biphasic_mmkin <- function() plot(f_uba_dfop_sfo_mixed)
+  plot_biphasic_mmkin <- function() plot(f_uba_dfop_sfo_mixed, pop_curve = TRUE)
   vdiffr::expect_doppelganger("mixed model fit for mmkin object", plot_biphasic_mmkin)
 
   plot_biphasic_saem_s <- function() plot(f_uba_dfop_sfo_saem)