From a8514e92fbb23d60db686ddf153592fb28c48a77 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 23 Mar 2023 16:42:31 +0100 Subject: Support covariates in endpoints() --- DESCRIPTION | 2 +- NEWS.md | 4 +++- R/endpoints.R | 33 +++++++++++++++++++++++++++------ R/saem.R | 4 ++-- log/check.log | 13 +++++++------ log/test.log | 35 +++++++++++++++++------------------ man/endpoints.Rd | 9 ++++++++- man/mkinfit.Rd | 4 ++-- man/parms.Rd | 13 ++++++++++++- man/plot.mixed.mmkin.Rd | 22 +++++++++++++++++++--- man/saem.Rd | 6 ------ man/set_nd_nq.Rd | 4 ++-- 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 ’ - -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") -- cgit v1.2.1