diff options
| -rw-r--r-- | R/illparms.R | 3 | ||||
| -rw-r--r-- | R/plot.mixed.mmkin.R | 116 | ||||
| -rw-r--r-- | log/check.log | 2 | ||||
| -rw-r--r-- | log/test.log | 18 | ||||
| -rw-r--r-- | man/illparms.Rd | 6 | ||||
| -rw-r--r-- | man/plot.mixed.mmkin.Rd | 5 | ||||
| -rw-r--r-- | tests/testthat/test_plot.R | 2 |
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) |