Scaled residual plots

10 files changed, 88 insertions, 45 deletions

diff --git a/NAMESPACE b/NAMESPACE
index 01abdcbb..8ea4c684 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -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)
diff --git a/NEWS.md b/NEWS.md
index cfb9fa96..395cd623 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -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
diff --git a/build.log b/build.log
index 41df1d5d..ffd56a90 100644
--- a/build.log
+++ b/build.log
@@ -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’
 
diff --git a/check.log b/check.log
index be2746c3..85055c74 100644
--- a/check.log
+++ b/check.log
@@ -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
diff --git a/test.log b/test.log
index cfc8d633..f0f09d41 100644
--- a/test.log
+++ b/test.log
@@ -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]