Improvements to mean_degparms() and plot.mixed.mmkin()

- New argument 'default_log_parms' for mean_degparms()
- 'plot.mixed.mmkin': Pass the frame argument also to residual plots,
  take the 'default_log_parms' argument for 'mean_degparms' used for
  constructing approximate population curves, plot population curve last
  to avoid that it is covered by data

3 files changed, 26 insertions, 10 deletions

diff --git a/R/dimethenamid_2018.R b/R/dimethenamid_2018.R
index 2fdd1981..00ed9073 100644
--- a/R/dimethenamid_2018.R
+++ b/R/dimethenamid_2018.R
@@ -49,12 +49,16 @@
 #' # look more plausible, but the truth is likely to be in
 #' # between these variants
 #' plot(mixed(dmta_sfo_sfo3p_tc), test_log_parms = TRUE)
-#' # Therefore we use nonlinear mixed-effects models
+#' # We can also specify a default value for the failing
+#' # log parameters, to mimic FOCUS guidance
+#' plot(mixed(dmta_sfo_sfo3p_tc), test_log_parms = TRUE,
+#'   default_log_parms = log(2)/1000)
+#' # As these attempts are not satisfying, we use nonlinear mixed-effects models
 #' # f_dmta_nlme_tc <- nlme(dmta_sfo_sfo3p_tc)
 #' # nlme reaches maxIter = 50 without convergence
 #' f_dmta_saem_tc <- saem(dmta_sfo_sfo3p_tc)
 #' # I am commenting out the convergence plot as rendering them
-#' # with pkgdown fails (at least without further tweaks to the 
+#' # with pkgdown fails (at least without further tweaks to the
 #' # graphics device used)
 #' #saemix::plot(f_dmta_saem_tc$so, plot.type = "convergence")
 #' summary(f_dmta_saem_tc)
@@ -65,6 +69,6 @@
 #' #   covariance.model = diag(c(0, rep(1, 7))))
 #' # saemix::plot(f_dmta_saem_tc_2$so, plot.type = "convergence")
 #' # This does not perform better judged by AIC and BIC
-#' saemix::compare.saemix(f_dmta_saem_tc$so, f_dmta_saem_tc_2$so)
+#' # saemix::compare.saemix(f_dmta_saem_tc$so, f_dmta_saem_tc_2$so)
 #' }
 "dimethenamid_2018"
diff --git a/R/mean_degparms.R b/R/mean_degparms.R
index ec20c068..fdcc5c00 100644
--- a/R/mean_degparms.R
+++ b/R/mean_degparms.R
@@ -11,8 +11,12 @@
 #'   rate constants) pass the t-test for significant difference from zero.
 #' @param conf.level Possibility to adjust the required confidence level
 #'   for parameter that are tested if requested by 'test_log_parms'.
+#' @param default_log_parms If set to a numeric value, this is used
+#'   as a default value for the tested log parameters that failed the
+#'   t-test.
 #' @export
-mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6)
+mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6,
+  default_log_parms = NA)
 {
   if (nrow(object) > 1) stop("Only row objects allowed")
   parm_mat_trans <- sapply(object, parms, transformed = TRUE)
@@ -33,7 +37,7 @@ mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.l
       parm_mat_trans_OK <- parm_mat_trans
       for (trans_parm in log_parm_trans_names) {
         parm_mat_trans_OK[trans_parm, ] <- ifelse(t_test_back_OK[trans_parm, ],
-          parm_mat_trans[trans_parm, ], NA)
+          parm_mat_trans[trans_parm, ], default_log_parms)
       }
     } else {
     parm_mat_trans_OK <- parm_mat_trans
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index 2903a05c..3a444253 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -14,6 +14,8 @@ utils::globalVariables("ds")
 #'   [mixed.mmkin] object
 #' @param conf.level Passed to [mean_degparms] in the case of an
 #'   [mixed.mmkin] object
+#' @param default_log_parms Passed to [mean_degparms] in the case of an
+#'   [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
@@ -69,6 +71,7 @@ plot.mixed.mmkin <- function(x,
   pred_over = NULL,
   test_log_parms = FALSE,
   conf.level = 0.6,
+  default_log_parms = NA,
   ymax = "auto", maxabs = "auto",
   ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)),
   nrow.legend = ceiling((length(i) + 1) / ncol.legend),
@@ -87,7 +90,8 @@ 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)
+    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)), ]))
@@ -247,8 +251,7 @@ plot.mixed.mmkin <- function(x,
       par(mar = c(3.0, 4.1, 1.1, 2.1))
     }
 
-    plot(pred_pop$time, pred_pop[[obs_var]],
-      type = "l", lwd = 2, lty = lty_pop,
+    plot(0, type = "n",
       xlim = xlim, ylim = ylim_row,
       xlab = xlab, ylab = paste("Residues", obs_var), frame = frame)
 
@@ -267,6 +270,9 @@ 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 (identical(maxabs, "auto")) {
       maxabs = max(abs(observed_row$residual), na.rm = TRUE)
     }
@@ -274,7 +280,8 @@ plot.mixed.mmkin <- function(x,
     if (identical(resplot, "time")) {
       plot(0, type = "n", xlim = xlim, xlab = "Time",
         ylim = c(-1.2 * maxabs, 1.2 * maxabs),
-        ylab = if (standardized) "Standardized residual" else "Residual")
+        ylab = if (standardized) "Standardized residual" else "Residual",
+        frame = frame)
 
       abline(h = 0, lty = 2)
 
@@ -289,7 +296,8 @@ plot.mixed.mmkin <- function(x,
         xlim = c(0, max(pred_ds[[obs_var]])),
         xlab = "Predicted",
         ylim = c(-1.2 * maxabs, 1.2 * maxabs),
-        ylab = if (standardized) "Standardized residual" else "Residual")
+        ylab = if (standardized) "Standardized residual" else "Residual",
+        frame = frame)
 
       abline(h = 0, lty = 2)