Check slopes in saemix covariate models

4 files changed, 43 insertions, 27 deletions

diff --git a/NEWS.md b/NEWS.md
index c7929fbe..1931142c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -6,6 +6,8 @@
 
 - 'R/parplot.R': Possibility to select the top 'llquant' fraction of the fits for the parameter plots, and improved legend text.
 
+- 'R/illparms.R': Also check if confidence intervals for slope parameters in covariate models include zero. Only implemented for fits obtained with the saemix backend.
+
 # mkin 1.2.1 (2022-11-19)
 
 - '{data,R}/ds_mixed.rda': Include the test data in the package instead of generating it in 'tests/testthat/setup_script.R'. Refactor the generating code to make it consistent and update tests.
diff --git a/R/illparms.R b/R/illparms.R
index 01e75cf1..eef4bd33 100644
--- a/R/illparms.R
+++ b/R/illparms.R
@@ -20,6 +20,9 @@
 #' @param random For hierarchical fits, should random effects be tested?
 #' @param errmod For hierarchical fits, should error model parameters be
 #' tested?
+#' @param slopes For hierarchical [saem] fits using saemix as backend,
+#' should slope parameters in the covariate model(starting with 'beta_') be
+#' tested?
 #' @return For [mkinfit] or [saem] objects, a character vector of parameter
 #' names. For [mmkin] or [mhmkin] objects, a matrix like object of class
 #' 'illparms.mmkin' or 'illparms.mhmkin'.
@@ -92,7 +95,7 @@ print.illparms.mmkin <- function(x, ...) {
 
 #' @rdname illparms
 #' @export
-illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) {
+illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, slopes = TRUE, ...) {
   if (inherits(object$so, "try-error")) {
     ill_parms <- NA
   } else {
@@ -106,6 +109,15 @@ illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod
       ill_parms_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "est."] > 0
       ill_parms <- c(ill_parms, names(which(ill_parms_errmod)))
     }
+    if (slopes) {
+      if (is.null(object$so)) stop("Slope testing is only implemented for the saemix backend")
+      slope_names <- grep("^beta_", object$so@model@name.fixed, value = TRUE)
+      ci <- object$so@results@conf.int
+      rownames(ci) <- ci$name
+      slope_ci <- ci[slope_names, ]
+      ill_parms_slopes <- slope_ci[, "lower"] < 0 & slope_ci[, "estimate"] > 0
+      ill_parms <- c(ill_parms, slope_names[ill_parms_slopes])
+    }
   }
   class(ill_parms) <- "illparms.saem.mmkin"
   return(ill_parms)
diff --git a/log/test.log b/log/test.log
index 7614b136..af2ffc41 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,23 +1,23 @@
 ℹ Testing mkin
 ✔ | F W S  OK | Context
 ✔ |         5 | AIC calculation
-✔ |         5 | Analytical solutions for coupled models [3.2s]
+✔ |         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.1s]
-✔ |     1  12 | Dimethenamid data from 2018 [32.0s]
+✔ |     1  12 | Dimethenamid data from 2018 [32.7s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_dmta.R:98'): 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 [4.8s]
+✔ |        14 | Error model fitting [5.0s]
 ✔ |         5 | Time step normalisation
 ✔ |         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.3s]
+✔ |         4 | Test fitting the decline of metabolites from their maximum [0.4s]
 ✔ |         1 | Fitting the logistic model [0.2s]
-✔ |        10 | Batch fitting and diagnosing hierarchical kinetic models [42.2s]
-✔ |     1  11 | Nonlinear mixed-effects models [13.3s]
+✔ |        10 | Batch fitting and diagnosing hierarchical kinetic models [43.4s]
+✔ |     1  11 | Nonlinear mixed-effects models [13.5s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_mixed.R:78'): saemix results are reproducible for biphasic fits
 Reason: Fitting with saemix takes around 10 minutes when using deSolve
@@ -27,40 +27,31 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |         3 | mkinfit features [0.7s]
 ✔ |         8 | mkinmod model generation and printing [0.2s]
 ✔ |         3 | Model predictions with mkinpredict [0.4s]
-✖ | 1      11 | Multistart method for saem.mmkin models [46.7s]
-────────────────────────────────────────────────────────────────────────────────
-Failure ('test_multistart.R:56'): multistart works for saem.mmkin models
-Snapshot of `testcase` to 'multistart/parplot-for-dfop-sfo-fit.svg' has changed
-Run `testthat::snapshot_review('multistart/')` to review changes
-Backtrace:
- 1. vdiffr::expect_doppelganger("parplot for dfop sfo fit", parplot_dfop_sfo)
-      at test_multistart.R:56:2
- 3. testthat::expect_snapshot_file(...)
-────────────────────────────────────────────────────────────────────────────────
+✔ |        12 | Multistart method for saem.mmkin models [50.1s]
 ✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.2s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [9.4s]
-✔ |        15 | Plotting [10.1s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [9.8s]
+✔ |        15 | Plotting [10.5s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |     1  36 | saemix parent models [73.6s]
+✔ |     1  36 | saemix parent models [85.7s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_saemix_parent.R:143'): 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
-✔ |        10 | Fitting the SFORB model [3.7s]
+✔ |        10 | Fitting the SFORB model [3.9s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         5 | Summary [0.2s]
-✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3s]
-✔ |         9 | Hypothesis tests [8.1s]
+✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
+✔ |         9 | Hypothesis tests [8.2s]
 ✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 260.9 s
+Duration: 280.1 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
 • Fitting this ODE model with saemix takes about 15 minutes on my system (1)
 • Fitting with saemix takes around 10 minutes when using deSolve (1)
 • This still takes almost 2.5 minutes although we do not solve ODEs (1)
 
-[ FAIL 1 | WARN 0 | SKIP 3 | PASS 269 ]
+[ FAIL 0 | WARN 0 | SKIP 3 | PASS 270 ]
diff --git a/man/illparms.Rd b/man/illparms.Rd
index 14be9c35..75eb18f0 100644
--- a/man/illparms.Rd
+++ b/man/illparms.Rd
@@ -22,7 +22,14 @@ illparms(object, ...)
 
 \method{print}{illparms.mmkin}(x, ...)
 
-\method{illparms}{saem.mmkin}(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...)
+\method{illparms}{saem.mmkin}(
+  object,
+  conf.level = 0.95,
+  random = TRUE,
+  errmod = TRUE,
+  slopes = TRUE,
+  ...
+)
 
 \method{print}{illparms.saem.mmkin}(x, ...)
 
@@ -43,6 +50,10 @@ illparms(object, ...)
 
 \item{errmod}{For hierarchical fits, should error model parameters be
 tested?}
+
+\item{slopes}{For hierarchical \link{saem} fits using saemix as backend,
+should slope parameters in the covariate model(starting with 'beta_') be
+tested?}
 }
 \value{
 For \link{mkinfit} or \link{saem} objects, a character vector of parameter