Support covariates in endpoints()

12 files changed, 100 insertions, 49 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 9322408d..4440af12 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,7 +2,7 @@ Package: mkin
 Type: Package
 Title: Kinetic Evaluation of Chemical Degradation Data
 Version: 1.2.3
-Date: 2023-02-17
+Date: 2023-03-23
 Authors@R: c(
   person("Johannes", "Ranke", role = c("aut", "cre", "cph"),
     email = "johannes.ranke@jrwb.de",
diff --git a/NEWS.md b/NEWS.md
index 65bbe185..02550282 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,9 @@
 
 - Depend on upcoming deSolve version with the possibility to avoid resolving symbols in a shared library (compiled models) over and over.
 
-- 'R/mkinerrmin.R': Fix typo in subset (use of = instead of ==), thanks to Sebastian Meyer for spotting this during his work on R 4.3.0
+- 'R/mkinerrmin.R': Fix typo in subset (use of = instead of ==), thanks to Sebastian Meyer for spotting this during his work on R 4.3.0.
+
+- 'R/{endpoints,parms,plot.mixed.mmkin}.R': Calculate parameters and endpoints and plot population curves for specific covariate values, or specific percentiles of covariate values used in saem fits.
 
 # mkin 1.2.2 (unreleased)
 
diff --git a/R/endpoints.R b/R/endpoints.R
index 9ea0e598..0aafd728 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -14,10 +14,11 @@
 #' another object that has list components mkinmod containing an [mkinmod]
 #' degradation model, and two numeric vectors, bparms.optim and bparms.fixed,
 #' that contain parameter values for that model.
-#' @param \dots Possibility to specify values for the covariates in the model.
-#' In case more than one vector is given, they either have to be of the same
-#' length, or of length one, in which case the respective covariate values are
-#' recycled.
+#' @param covariates Numeric vector with covariate values for all variables in
+#' any covariate models in the object. If given, it overrides 'covariate_quantile'.
+#' @param covariate_quantile This argument only has an effect if the fitted
+#' object has covariate models. If so, the default is to show endpoints
+#' for the median of the covariate values (50th percentile).
 #' @importFrom stats optimize
 #' @return A list with a matrix of dissipation times named distimes, and, if
 #' applicable, a vector of formation fractions named ff and, if the SFORB model
@@ -38,14 +39,34 @@
 #'   }
 #'
 #' @export
-endpoints <- function(fit, ..., covariates = mean) {
+endpoints <- function(fit, covariates = NULL, covariate_quantile = 0.5) {
   mkinmod <- fit$mkinmod
   obs_vars <- names(mkinmod$spec)
 
-  degparms <- c(fit$bparms.optim, fit$bparms.fixed)
+  if (!is.null(fit$covariate_models)) {
+    if (is.null(covariates)) {
+      covariates = as.data.frame(
+        apply(fit$covariates, 2, quantile,
+          covariate_quantile, simplify = FALSE))
+    } else {
+      covariates <- data.frame("User" = covariates)
+    }
+    degparms_trans <- parms(fit, covariates = covariates)[, 1]
+    if (inherits(fit, "saem.mmkin") & (fit$transformations == "saemix")) {
+      degparms <- degparms_trans
+    } else {
+      degparms <- backtransform_odeparms(degparms_trans,
+        fit$mkinmod,
+        transform_rates = fit$transform_rates,
+        transform_fractions = fit$transform_fractions)
+    }
+  } else {
+    degparms <- c(fit$bparms.optim, fit$bparms.fixed)
+  }
 
   # Set up object to return
   ep <- list()
+  ep$covariates <- covariates
   ep$ff <- vector()
   ep$SFORB <- vector()
   ep$distimes <- data.frame(
diff --git a/R/saem.R b/R/saem.R
index 865f174e..2fa770bb 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -802,14 +802,14 @@ update.saem.mmkin <- function(object, ..., evaluate = TRUE) {
 }
 
 #' @export
-#' @rdname saem
+#' @rdname parms
 #' @param ci Should a matrix with estimates and confidence interval boundaries
 #' be returned? If FALSE (default), a vector of estimates is returned if no
 #' covariates are given, otherwise a matrix of estimates is returned, with
 #' each column corresponding to a row of the data frame holding the covariates
 #' @param covariates A data frame holding covariate values for which to
 #' return parameter values. Only has an effect if 'ci' is FALSE.
-parms.saem.mmkin <- function(object, ci = FALSE, covariates = NULL, covariate_quantiles = ...) {
+parms.saem.mmkin <- function(object, ci = FALSE, covariates = NULL, ...) {
   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) +
diff --git a/log/check.log b/log/check.log
index b44bc1b4..f7a0d87e 100644
--- a/log/check.log
+++ b/log/check.log
@@ -7,12 +7,12 @@
 * checking extension type ... Package
 * this is package ‘mkin’ version ‘1.2.3’
 * package encoding: UTF-8
-* checking CRAN incoming feasibility ... NOTE
+* checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
 Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’
-
-The Date field is over a month old.
 * checking package namespace information ... OK
-* checking package dependencies ... OK
+* checking package dependencies ...Warning: unable to access index for repository https://cran.rstudio.com/src/contrib:
+  cannot open URL 'https://cran.rstudio.com/src/contrib/PACKAGES'
+ OK
 * checking if this is a source package ... OK
 * checking if there is a namespace ... OK
 * checking for executable files ... OK
@@ -23,7 +23,8 @@ The Date field is over a month old.
 * checking whether package ‘mkin’ can be installed ... OK
 * checking installed package size ... OK
 * checking package directory ... OK
-* checking for future file timestamps ... OK
+* checking for future file timestamps ... NOTE
+unable to verify current time
 * checking ‘build’ directory ... OK
 * checking DESCRIPTION meta-information ... OK
 * checking top-level files ... OK
@@ -65,7 +66,7 @@ The Date field is over a month old.
 * checking tests ... SKIPPED
 * checking for unstated dependencies in vignettes ... OK
 * checking package vignettes in ‘inst/doc’ ... OK
-* checking re-building of vignette outputs ... [12s/12s] OK
+* checking re-building of vignette outputs ... [11s/11s] OK
 * checking PDF version of manual ... OK
 * skipping checking HTML version of manual: no command ‘tidy’ found
 * checking for non-standard things in the check directory ... OK
diff --git a/log/test.log b/log/test.log
index d4fb85a2..d3939c0f 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,14 +1,13 @@
-ℹ Loading mkin
 ℹ Testing mkin
 ✔ | F W S  OK | Context
 ✔ |         5 | AIC calculation
 ✔ |         5 | Analytical solutions for coupled models [2.8s]
 ✔ |         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.3s]
+✔ |        12 | Confidence intervals and p-values [0.9s]
+✔ |     1  12 | Dimethenamid data from 2018 [30.7s]
 ────────────────────────────────────────────────────────────────────────────────
-Skip ('test_dmta.R:99'): Different backends get consistent results for SFO-SFO3+, dimethenamid data
+Skip ('test_dmta.R:99: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 [6.9s]
@@ -18,9 +17,9 @@ Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
 ✔ |         4 | Test fitting the decline of metabolites from their maximum [0.3s]
 ✔ |         1 | Fitting the logistic model [0.2s]
 ✔ |        10 | Batch fitting and diagnosing hierarchical kinetic models [42.7s]
-✔ |     1  11 | Nonlinear mixed-effects models [13.1s]
+✔ |     1  11 | Nonlinear mixed-effects models [13.3s]
 ────────────────────────────────────────────────────────────────────────────────
-Skip ('test_mixed.R:78'): saemix results are reproducible for biphasic fits
+Skip ('test_mixed.R:78: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
@@ -28,27 +27,27 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |         3 | mkinfit features [0.6s]
 ✔ |         8 | mkinmod model generation and printing [0.2s]
 ✔ |         3 | Model predictions with mkinpredict [0.3s]
-✔ |        12 | Multistart method for saem.mmkin models [74.1s]
-✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.9s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [8.5s]
-✔ |        15 | Plotting [11.4s]
+✔ |        12 | Multistart method for saem.mmkin models [73.3s]
+✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.0s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [1107.1s]
+✔ |        15 | Plotting [13.3s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |     1  36 | saemix parent models [73.9s]
+✔ |     1  36 | saemix parent models [74.8s]
 ────────────────────────────────────────────────────────────────────────────────
-Skip ('test_saemix_parent.R:143'): We can also use mkin solution methods for saem
+Skip ('test_saemix_parent.R:143: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.1s]
+✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5s]
 ✔ |        11 | Processing of residue series
-✔ |        10 | Fitting the SFORB model [3.1s]
+✔ |        10 | Fitting the SFORB model [3.3s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         5 | Summary [0.1s]
-✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [1.7s]
-✔ |         9 | Hypothesis tests [6.0s]
-✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [1.9s]
+✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [1.8s]
+✔ |         9 | Hypothesis tests [6.1s]
+✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [1.8s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 285.4 s
+Duration: 1386.7 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
 • Fitting this ODE model with saemix takes about 15 minutes on my system (1)
diff --git a/man/endpoints.Rd b/man/endpoints.Rd
index e3663aec..0df3616c 100644
--- a/man/endpoints.Rd
+++ b/man/endpoints.Rd
@@ -5,13 +5,20 @@
 \title{Function to calculate endpoints for further use from kinetic models fitted
 with mkinfit}
 \usage{
-endpoints(fit)
+endpoints(fit, covariates = NULL, covariate_quantile = 0.5)
 }
 \arguments{
 \item{fit}{An object of class \link{mkinfit}, \link{nlme.mmkin} or \link{saem.mmkin}, or
 another object that has list components mkinmod containing an \link{mkinmod}
 degradation model, and two numeric vectors, bparms.optim and bparms.fixed,
 that contain parameter values for that model.}
+
+\item{covariates}{Numeric vector with covariate values for all variables in
+any covariate models in the object. If given, it overrides 'covariate_quantile'.}
+
+\item{covariate_quantile}{This argument only has an effect if the fitted
+object has covariate models. If so, the default is to show endpoints
+for the median of the covariate values (50th percentile).}
 }
 \value{
 A list with a matrix of dissipation times named distimes, and, if
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index b5b24449..f96b4d22 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -23,8 +23,8 @@ mkinfit(
   atol = 1e-08,
   rtol = 1e-10,
   error_model = c("const", "obs", "tc"),
-  error_model_algorithm = c("auto", "d_3", "direct", "twostep", "threestep",
-    "fourstep", "IRLS", "OLS"),
+  error_model_algorithm = c("auto", "d_3", "direct", "twostep", "threestep", "fourstep",
+    "IRLS", "OLS"),
   reweight.tol = 1e-08,
   reweight.max.iter = 10,
   trace_parms = FALSE,
diff --git a/man/parms.Rd b/man/parms.Rd
index 5c0e8895..b1551374 100644
--- a/man/parms.Rd
+++ b/man/parms.Rd
@@ -1,10 +1,11 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/parms.R
+% Please edit documentation in R/parms.R, R/saem.R
 \name{parms}
 \alias{parms}
 \alias{parms.mkinfit}
 \alias{parms.mmkin}
 \alias{parms.multistart}
+\alias{parms.saem.mmkin}
 \title{Extract model parameters}
 \usage{
 parms(object, ...)
@@ -14,6 +15,8 @@ parms(object, ...)
 \method{parms}{mmkin}(object, transformed = FALSE, errparms = TRUE, ...)
 
 \method{parms}{multistart}(object, exclude_failed = TRUE, ...)
+
+\method{parms}{saem.mmkin}(object, ci = FALSE, covariates = NULL, ...)
 }
 \arguments{
 \item{object}{A fitted model object.}
@@ -28,6 +31,14 @@ 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?}
+
+\item{ci}{Should a matrix with estimates and confidence interval boundaries
+be returned? If FALSE (default), a vector of estimates is returned if no
+covariates are given, otherwise a matrix of estimates is returned, with
+each column corresponding to a row of the data frame holding the covariates}
+
+\item{covariates}{A data frame holding covariate values for which to
+return parameter values. Only has an effect if 'ci' is FALSE.}
 }
 \value{
 Depending on the object, a numeric vector of fitted model parameters,
diff --git a/man/plot.mixed.mmkin.Rd b/man/plot.mixed.mmkin.Rd
index 9c4474ff..1e264db3 100644
--- a/man/plot.mixed.mmkin.Rd
+++ b/man/plot.mixed.mmkin.Rd
@@ -9,10 +9,12 @@
   i = 1:ncol(x$mmkin),
   obs_vars = names(x$mkinmod$map),
   standardized = TRUE,
+  covariates = NULL,
+  covariate_quantiles = c(0.5, 0.05, 0.95),
   xlab = "Time",
   xlim = range(x$data$time),
   resplot = c("predicted", "time"),
-  pop_curve = "auto",
+  pop_curves = "auto",
   pred_over = NULL,
   test_log_parms = FALSE,
   conf.level = 0.6,
@@ -43,6 +45,16 @@ variables in the model.}
 \item{standardized}{Should the residuals be standardized? Only takes effect if
 \code{resplot = "time"}.}
 
+\item{covariates}{Data frame with covariate values for all variables in
+any covariate models in the object. If given, it overrides 'covariate_quantiles'.
+Each line in the data frame will result in a line drawn for the population.
+Rownames are used in the legend to label the lines.}
+
+\item{covariate_quantiles}{This argument only has an effect if the fitted
+object has covariate models. If so, the default is to show three population
+curves, for the 5th percentile, the 50th percentile and the 95th percentile
+of the covariate values used for fitting the model.}
+
 \item{xlab}{Label for the x axis.}
 
 \item{xlim}{Plot range in x direction.}
@@ -50,9 +62,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
+\item{pop_curves}{Per default, one 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.}
+In case there is a covariate model, the behaviour depends on the value
+of 'covariates'}
 
 \item{pred_over}{Named list of alternative predictions as obtained
 from \link{mkinpredict} with a compatible \link{mkinmod}.}
@@ -96,6 +109,9 @@ The function is called for its side effect.
 \description{
 Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object
 }
+\note{
+Covariate models are currently only supported for saem.mmkin objects.
+}
 \examples{
 ds <- lapply(experimental_data_for_UBA_2019[6:10],
  function(x) x$data[c("name", "time", "value")])
diff --git a/man/saem.Rd b/man/saem.Rd
index 89647ff4..fd19fe59 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -6,7 +6,6 @@
 \alias{print.saem.mmkin}
 \alias{saemix_model}
 \alias{saemix_data}
-\alias{parms.saem.mmkin}
 \title{Fit nonlinear mixed models with SAEM}
 \usage{
 saem(object, ...)
@@ -54,8 +53,6 @@ saemix_model(
 )
 
 saemix_data(object, covariates = NULL, verbose = FALSE, ...)
-
-\method{parms}{saem.mmkin}(object, ci = FALSE, ...)
 }
 \arguments{
 \item{object}{An \link{mmkin} row object containing several fits of the same
@@ -124,9 +121,6 @@ and the end of the optimisation process?}
 \item{x}{An saem.mmkin object to print}
 
 \item{digits}{Number of digits to use for printing}
-
-\item{ci}{Should a matrix with estimates and confidence interval boundaries
-be returned? If FALSE (default), a vector of estimates is returned.}
 }
 \value{
 An S3 object of class 'saem.mmkin', containing the fitted
diff --git a/man/set_nd_nq.Rd b/man/set_nd_nq.Rd
index 796c27a2..87a3fae1 100644
--- a/man/set_nd_nq.Rd
+++ b/man/set_nd_nq.Rd
@@ -54,9 +54,9 @@ it automates the proposal of Boesten et al (2015).
 }
 \section{Functions}{
 \itemize{
-\item \code{set_nd_nq_focus}: Set non-detects in residue time series according to FOCUS rules
-}}
+\item \code{set_nd_nq_focus()}: Set non-detects in residue time series according to FOCUS rules
 
+}}
 \examples{
 # FOCUS (2014) p. 75/76 and 131/132
 parent_1 <- c(.12, .09, .05, .03, "nd", "nd", "nd", "nd", "nd", "nd")