aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/mkinresplot.R38
-rw-r--r--R/plot.mkinfit.R35
2 files changed, 50 insertions, 23 deletions
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

Contact - Imprint