diff options
| -rw-r--r-- | NAMESPACE | 1 | ||||
| -rw-r--r-- | NEWS.md | 2 | ||||
| -rw-r--r-- | R/mkinresplot.R | 38 | ||||
| -rw-r--r-- | R/plot.mkinfit.R | 35 | ||||
| -rw-r--r-- | build.log | 2 | ||||
| -rw-r--r-- | check.log | 2 | ||||
| -rw-r--r-- | man/mkinresplot.Rd | 16 | ||||
| -rw-r--r-- | man/plot.mkinfit.Rd | 15 | ||||
| -rw-r--r-- | test.log | 20 | ||||
| -rw-r--r-- | tests/testthat/FOCUS_2006_D.csf | 2 |
10 files changed, 88 insertions, 45 deletions
@@ -84,6 +84,7 @@ importFrom(stats,qchisq) importFrom(stats,qf) importFrom(stats,qnorm) importFrom(stats,qt) +importFrom(stats,residuals) importFrom(stats,rnorm) importFrom(stats,update) importFrom(utils,write.table) @@ -1,5 +1,7 @@ # mkin 0.9.49.8 (unreleased) +- 'plot_res', 'plot_sep' and 'mkinerrplot': Add the possibility to show standardized residuals and make it the default for fits with error models other than 'const' + - 'lrtest.mkinfit': Improve naming of the compared fits in the case of fixed parameters - 'confint.mkinfit': Make the quadratic approximation the default, as the likelihood profiling takes a lot of time, especially if the fit has more than three parameters diff --git a/R/mkinresplot.R b/R/mkinresplot.R index 5377dbf2..0bfdd02f 100644 --- a/R/mkinresplot.R +++ b/R/mkinresplot.R @@ -1,23 +1,25 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual")) #' Function to plot residuals stored in an mkin object -#' +#' #' This function plots the residuals for the specified subset of the observed #' variables from an mkinfit object. A combined plot of the fitted model and #' the residuals can be obtained using \code{\link{plot.mkinfit}} using the #' argument \code{show_residuals = TRUE}. -#' +#' +#' @importFrom stats residuals #' @param object A fit represented in an \code{\link{mkinfit}} object. #' @param obs_vars A character vector of names of the observed variables for #' which residuals should be plotted. Defaults to all observed variables in #' the model #' @param xlim plot range in x direction. -#' @param xlab Label for the x axis. Defaults to "Time [days]". -#' @param ylab Label for the y axis. Defaults to "Residual [\% of applied -#' radioactivity]". +#' @param xlab Label for the x axis. +#' @param standardized Should the residuals be standardized by dividing by the +#' standard deviation given by the error model of the fit? +#' @param ylab Label for the y axis. #' @param maxabs Maximum absolute value of the residuals. This is used for the #' scaling of the y axis and defaults to "auto". -#' @param legend Should a legend be plotted? Defaults to "TRUE". +#' @param legend Should a legend be plotted? #' @param lpos Where should the legend be placed? Default is "topright". Will #' be passed on to \code{\link{legend}}. #' @param col_obs Colors for the observed variables. @@ -28,19 +30,21 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual")) #' effect, namely to produce a plot. #' @author Johannes Ranke #' @seealso \code{\link{mkinplot}}, for a way to plot the data and the fitted -#' lines of the mkinfit object. +#' lines of the mkinfit object, and \code{\link{plot_res}} for a function +#' combining the plot of the fit and the residual plot. #' @examples -#' +#' #' model <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO")) #' fit <- mkinfit(model, FOCUS_2006_D, quiet = TRUE) #' mkinresplot(fit, "m1") -#' +#' #' @export mkinresplot <- function (object, obs_vars = names(object$mkinmod$map), xlim = c(0, 1.1 * max(object$data$time)), - xlab = "Time", ylab = "Residual", - maxabs = "auto", legend= TRUE, lpos = "topright", + standardized = FALSE, + xlab = "Time", ylab = ifelse(standardized, "Standardized residual", "Residual"), + maxabs = "auto", legend = TRUE, lpos = "topright", col_obs = "auto", pch_obs = "auto", frame = TRUE, ...) @@ -51,9 +55,15 @@ mkinresplot <- function (object, obs_vars <- intersect(obs_vars_all, obs_vars) } else obs_vars <- obs_vars_all - residuals <- subset(object$data, variable %in% obs_vars, residual) + if (standardized) { + res_col <- "standardized" + object$data[[res_col]] <- residuals(object, standardized = TRUE) + } else { + res_col <- "residual" + } + res <- subset(object$data, variable %in% obs_vars)[res_col] - if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE) + if (maxabs == "auto") maxabs = max(abs(res), na.rm = TRUE) # Set colors and symbols if (col_obs[1] == "auto") { @@ -71,7 +81,7 @@ mkinresplot <- function (object, ylim = c(-1.2 * maxabs, 1.2 * maxabs), ...) for(obs_var in obs_vars){ - residuals_plot <- subset(object$data, variable == obs_var, c("time", "residual")) + residuals_plot <- subset(object$data, variable == obs_var, c("time", res_col)) points(residuals_plot, pch = pch_obs[obs_var], col = col_obs[obs_var]) } diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R index 16415a3b..e0b3ad13 100644 --- a/R/plot.mkinfit.R +++ b/R/plot.mkinfit.R @@ -30,7 +30,10 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "variable", "obse #' @param show_residuals Should residuals be shown? If only one plot of the #' fits is shown, the residual plot is in the lower third of the plot. #' Otherwise, i.e. if "sep_obs" is given, the residual plots will be located -#' to the right of the plots of the fitted curves. +#' to the right of the plots of the fitted curves. If this is set to +#' 'standardized', a plot of the residuals divided by the standard deviation +#' given by the fitted error model will be shown. +#' @param standardized For #' @param show_errplot Should squared residuals and the error model be shown? #' If only one plot of the fits is shown, this plot is in the lower third of #' the plot. Otherwise, i.e. if "sep_obs" is given, the residual plots will @@ -68,6 +71,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "variable", "obse #' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE, error_model = "tc") #' plot(fit) #' plot_res(fit) +#' plot_res(fit, standardized = FALSE) #' plot_err(fit) #' #' # Show the observed variables separately, with residuals @@ -101,11 +105,19 @@ plot.mkinfit <- function(x, fit = x, show_errmin = FALSE, errmin_digits = 3, frame = TRUE, ...) { + if (identical(show_residuals, "standardized")) { + show_residuals <- TRUE + standardized <- TRUE + } else { + standardized <- FALSE + } + if (add && show_residuals) stop("If adding to an existing plot we can not show residuals") if (add && show_errplot) stop("If adding to an existing plot we can not show the error model plot") if (show_residuals && show_errplot) stop("We can either show residuals over time or the error model plot, not both") if (add && sep_obs) stop("If adding to an existing plot we can not show observed variables separately") + solution_type = fit$solution_type parms.all <- c(fit$bparms.optim, fit$bparms.fixed) @@ -260,14 +272,16 @@ plot.mkinfit <- function(x, fit = x, # Show residuals if requested if (show_residuals) { - mkinresplot(fit, obs_vars = row_obs_vars, pch_obs = pch_obs[row_obs_vars], col_obs = col_obs[row_obs_vars], - legend = FALSE, frame = frame) + mkinresplot(fit, obs_vars = row_obs_vars, standardized = standardized, + pch_obs = pch_obs[row_obs_vars], col_obs = col_obs[row_obs_vars], + legend = FALSE, frame = frame) } # Show error model plot if requested if (show_errplot) { - mkinerrplot(fit, obs_vars = row_obs_vars, pch_obs = pch_obs[row_obs_vars], col_obs = col_obs[row_obs_vars], - legend = FALSE, frame = frame) + mkinerrplot(fit, obs_vars = row_obs_vars, + pch_obs = pch_obs[row_obs_vars], col_obs = col_obs[row_obs_vars], + legend = FALSE, frame = frame) } } if (do_layout) par(oldpar, no.readonly = TRUE) @@ -275,16 +289,19 @@ plot.mkinfit <- function(x, fit = x, #' @rdname plot.mkinfit #' @export -plot_sep <- function(fit, show_errmin = TRUE, ...) { - plot.mkinfit(fit, sep_obs = TRUE, show_residuals = TRUE, +plot_sep <- function(fit, show_errmin = TRUE, + show_residuals = ifelse(identical(fit$err_mod, "const"), TRUE, "standardized"), ...) { + plot.mkinfit(fit, sep_obs = TRUE, show_residuals = show_residuals, show_errmin = show_errmin, ...) } #' @rdname plot.mkinfit #' @export -plot_res <- function(fit, sep_obs = FALSE, show_errmin = sep_obs, ...) { +plot_res <- function(fit, sep_obs = FALSE, show_errmin = sep_obs, + show_residuals = ifelse(identical(fit$err_mod, "const"), TRUE, "standardized"), ...) +{ plot.mkinfit(fit, sep_obs = sep_obs, show_errmin = show_errmin, - show_residuals = TRUE, row_layout = TRUE, ...) + show_residuals = show_residuals, row_layout = TRUE, ...) } #' @rdname plot.mkinfit @@ -6,5 +6,5 @@ * checking for LF line-endings in source and make files and shell scripts * checking for empty or unneeded directories * looking to see if a ‘data/datalist’ file should be added -* building ‘mkin_0.9.49.7.tar.gz’ +* building ‘mkin_0.9.49.8.tar.gz’ @@ -5,7 +5,7 @@ * using options ‘--no-tests --as-cran’ * checking for file ‘mkin/DESCRIPTION’ ... OK * checking extension type ... Package -* this is package ‘mkin’ version ‘0.9.49.7’ +* this is package ‘mkin’ version ‘0.9.49.8’ * package encoding: UTF-8 * checking CRAN incoming feasibility ... Note_to_CRAN_maintainers Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’ diff --git a/man/mkinresplot.Rd b/man/mkinresplot.Rd index 27e1322f..465b3038 100644 --- a/man/mkinresplot.Rd +++ b/man/mkinresplot.Rd @@ -5,7 +5,8 @@ \title{Function to plot residuals stored in an mkin object} \usage{ mkinresplot(object, obs_vars = names(object$mkinmod$map), xlim = c(0, - 1.1 * max(object$data$time)), xlab = "Time", ylab = "Residual", + 1.1 * max(object$data$time)), standardized = FALSE, xlab = "Time", + ylab = ifelse(standardized, "Standardized residual", "Residual"), maxabs = "auto", legend = TRUE, lpos = "topright", col_obs = "auto", pch_obs = "auto", frame = TRUE, ...) } @@ -18,15 +19,17 @@ the model} \item{xlim}{plot range in x direction.} -\item{xlab}{Label for the x axis. Defaults to "Time [days]".} +\item{standardized}{Should the residuals be standardized by dividing by the +standard deviation given by the error model of the fit?} -\item{ylab}{Label for the y axis. Defaults to "Residual [\% of applied -radioactivity]".} +\item{xlab}{Label for the x axis.} + +\item{ylab}{Label for the y axis.} \item{maxabs}{Maximum absolute value of the residuals. This is used for the scaling of the y axis and defaults to "auto".} -\item{legend}{Should a legend be plotted? Defaults to "TRUE".} +\item{legend}{Should a legend be plotted?} \item{lpos}{Where should the legend be placed? Default is "topright". Will be passed on to \code{\link{legend}}.} @@ -58,7 +61,8 @@ mkinresplot(fit, "m1") } \seealso{ \code{\link{mkinplot}}, for a way to plot the data and the fitted - lines of the mkinfit object. + lines of the mkinfit object, and \code{\link{plot_res}} for a function + combining the plot of the fit and the residual plot. } \author{ Johannes Ranke diff --git a/man/plot.mkinfit.Rd b/man/plot.mkinfit.Rd index 3834eaf5..6824f8a5 100644 --- a/man/plot.mkinfit.Rd +++ b/man/plot.mkinfit.Rd @@ -16,9 +16,13 @@ lpos = "topright", inset = c(0.05, 0.05), show_errmin = FALSE, errmin_digits = 3, frame = TRUE, ...) -plot_sep(fit, show_errmin = TRUE, ...) +plot_sep(fit, show_errmin = TRUE, + show_residuals = ifelse(identical(fit$err_mod, "const"), TRUE, + "standardized"), ...) -plot_res(fit, sep_obs = FALSE, show_errmin = sep_obs, ...) +plot_res(fit, sep_obs = FALSE, show_errmin = sep_obs, + show_residuals = ifelse(identical(fit$err_mod, "const"), TRUE, + "standardized"), ...) plot_err(fit, sep_obs = FALSE, show_errmin = sep_obs, ...) } @@ -54,7 +58,9 @@ corresponding model prediction lines.} \item{show_residuals}{Should residuals be shown? If only one plot of the fits is shown, the residual plot is in the lower third of the plot. Otherwise, i.e. if "sep_obs" is given, the residual plots will be located -to the right of the plots of the fitted curves.} +to the right of the plots of the fitted curves. If this is set to +'standardized', a plot of the residuals divided by the standard deviation + given by the fitted error model will be shown.} \item{show_errplot}{Should squared residuals and the error model be shown? If only one plot of the fits is shown, this plot is in the lower third of @@ -89,6 +95,8 @@ chi2 error percentage.} \item{frame}{Should a frame be drawn around the plots?} \item{\dots}{Further arguments passed to \code{\link{plot}}.} + +\item{standardized}{For} } \value{ The function is called for its side effect. @@ -113,6 +121,7 @@ SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1", full = "Parent"), fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE, error_model = "tc") plot(fit) plot_res(fit) +plot_res(fit, standardized = FALSE) plot_err(fit) # Show the observed variables separately, with residuals @@ -2,30 +2,30 @@ Loading mkin Testing mkin ✔ | OK F W S | Context ✔ | 2 | Export dataset for reading into CAKE -✔ | 10 | Confidence intervals and p-values [9.8 s] -✔ | 14 | Error model fitting [38.8 s] -✔ | 4 | Calculation of FOCUS chi2 error levels [2.3 s] -✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.5 s] -✔ | 6 | Test fitting the decline of metabolites from their maximum [0.8 s] +✔ | 10 | Confidence intervals and p-values [9.7 s] +✔ | 14 | Error model fitting [36.6 s] +✔ | 4 | Calculation of FOCUS chi2 error levels [2.2 s] +✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.3 s] +✔ | 6 | Test fitting the decline of metabolites from their maximum [0.7 s] ✔ | 1 | Fitting the logistic model [0.9 s] ✔ | 1 | Test dataset class mkinds used in gmkin ✔ | 12 | Special cases of mkinfit calls [2.4 s] ✔ | 9 | mkinmod model generation and printing [0.2 s] ✔ | 3 | Model predictions with mkinpredict [0.3 s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [4.2 s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [4.0 s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.3 s] ✔ | 3 | Summary ✔ | 11 | Plotting [0.6 s] ✔ | 4 | AIC calculation ✔ | 2 | Residuals extracted from mkinfit models -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [5.5 s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [5.2 s] ✔ | 4 | Fitting the SFORB model [1.7 s] ✔ | 1 | Summaries of old mkinfit objects -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [7.3 s] -✔ | 6 | Hypothesis tests [32.5 s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [7.1 s] +✔ | 6 | Hypothesis tests [31.2 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 113.3 s +Duration: 108.3 s OK: 132 Failed: 0 diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf index a28aafec..d5d807ab 100644 --- a/tests/testthat/FOCUS_2006_D.csf +++ b/tests/testthat/FOCUS_2006_D.csf @@ -5,7 +5,7 @@ Description: MeasurementUnits: % AR TimeUnits: days Comments: Created using mkin::CAKE_export -Date: 2019-11-01 +Date: 2019-11-04 Optimiser: IRLS [Data] |