Add functionality to plot the error model

by plotting squared residuals against predicted values, and
showing the variance function used in the fitted error model.

Rebuild docs

4 files changed, 116 insertions, 15 deletions

diff --git a/R/mkinerrplot.R b/R/mkinerrplot.R
new file mode 100644
index 00000000..a2adcefa
--- /dev/null
+++ b/R/mkinerrplot.R
@@ -0,0 +1,84 @@
+# Copyright (C) 2008-2014,2019 Johannes Ranke
+# Contact: jranke@uni-bremen.de
+
+# This file is part of the R package mkin
+
+# mkin is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual"))
+
+mkinerrplot <- function (object,
+  obs_vars = names(object$mkinmod$map),
+  xlim = c(0, 1.1 * max(object$data$predicted)),
+  xlab = "Predicted", ylab = "Squared residual",
+  maxy = "auto", legend= TRUE, lpos = "topright", 
+  col_obs = "auto", pch_obs = "auto",
+  ...)
+{
+  obs_vars_all <- as.character(unique(object$data$variable))
+
+  if (length(obs_vars) > 0){
+      obs_vars <- intersect(obs_vars_all, obs_vars)
+  } else obs_vars <- obs_vars_all
+
+  residuals <- subset(object$data, variable %in% obs_vars, residual)
+
+  if (maxy == "auto") maxy = max(residuals^2, na.rm = TRUE)
+
+  # Set colors and symbols
+  if (col_obs[1] == "auto") {
+    col_obs <- 1:length(obs_vars)
+  }
+
+  if (pch_obs[1] == "auto") {
+    pch_obs <- 1:length(obs_vars)
+  }
+  names(col_obs) <- names(pch_obs) <- obs_vars
+
+  plot(0, type = "n",
+       xlab = xlab, ylab = ylab,
+       xlim = xlim,
+       ylim = c(0, 1.2 * maxy), ...)
+
+  for(obs_var in obs_vars){
+    residuals_plot <- subset(object$data, variable == obs_var, c("predicted", "residual"))
+    points(residuals_plot[["predicted"]],
+           residuals_plot[["residual"]]^2,
+           pch = pch_obs[obs_var], col = col_obs[obs_var])
+  }
+
+  if (object$err_mod == "const") {
+    abline(h = object$errparms^2, lty = 2, col = col_obs[obs_var])
+  }
+  if (object$err_mod == "obs") {
+    for (obs_var in obs_vars) {
+      sigma_name = paste0("sigma_", obs_var)
+      abline(h = object$errparms[sigma_name]^2, lty = 2,
+             col = col_obs[obs_var])
+    }
+  }
+  if (object$err_mod == "tc") {
+    sigma_plot <- function(predicted) {
+      sigma_twocomp(predicted,
+                    sigma_low = object$errparms[1],
+                    rsd_high = object$errparms[2])^2
+    }
+    plot(sigma_plot, from = 0, to = max(object$data$predicted),
+         add = TRUE, lty = 2, col = col_obs[obs_var])
+  }
+
+  if (legend == TRUE) {
+    legend(lpos, inset = c(0.05, 0.05), legend = obs_vars,
+      col = col_obs[obs_vars], pch = pch_obs[obs_vars])
+  }
+}
diff --git a/R/mkinresplot.R b/R/mkinresplot.R
index 3650ef4b..739a80e9 100644
--- a/R/mkinresplot.R
+++ b/R/mkinresplot.R
@@ -23,7 +23,7 @@ mkinresplot <- function (object,
   xlab = "Time", ylab = "Residual",

   maxabs = "auto", legend= TRUE, lpos = "topright", ...)

 {

-	obs_vars_all <- as.character(unique(object$data$variable))

+  obs_vars_all <- as.character(unique(object$data$variable))

 

   if (length(obs_vars) > 0){

       obs_vars <- intersect(obs_vars_all, obs_vars)

@@ -33,18 +33,18 @@ mkinresplot <- function (object,
 

   if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE)

 

-	col_obs <- pch_obs <- 1:length(obs_vars)

- 	names(col_obs) <- names(pch_obs) <- obs_vars

+  col_obs <- pch_obs <- 1:length(obs_vars)

+  names(col_obs) <- names(pch_obs) <- obs_vars

 

   plot(0, type = "n",

        xlab = xlab, ylab = ylab,

        xlim = xlim,

        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"))

-		points(residuals_plot, pch = pch_obs[obs_var], col = col_obs[obs_var])

-	}

+  for(obs_var in obs_vars){

+    residuals_plot <- subset(object$data, variable == obs_var, c("time", "residual"))

+    points(residuals_plot, pch = pch_obs[obs_var], col = col_obs[obs_var])

+  }

 

   abline(h = 0, lty = 2)

 

diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R
index ee836eb8..df9888e7 100644
--- a/R/plot.mkinfit.R
+++ b/R/plot.mkinfit.R
@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2016 Johannes Ranke
+# Copyright (C) 2010-2016,2019 Johannes Ranke
 # Contact: jranke@uni-bremen.de
 
 # This file is part of the R package mkin
@@ -26,12 +26,16 @@ plot.mkinfit <- function(x, fit = x,
   pch_obs = col_obs,
   lty_obs = rep(1, length(obs_vars)),
   add = FALSE, legend = !add,
-  show_residuals = FALSE, maxabs = "auto",
+  show_residuals = FALSE,
+  show_errplot = FALSE,
+  maxabs = "auto",
   sep_obs = FALSE, rel.height.middle = 0.9,
   lpos = "topright", inset = c(0.05, 0.05),
   show_errmin = FALSE, errmin_digits = 3, ...)
 {
   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
@@ -68,8 +72,7 @@ plot.mkinfit <- function(x, fit = x,
   # Create a plot layout only if not to be added to an existing plot
   # or only a single plot is requested (e.g. by plot.mmkin)
   do_layout = FALSE
-  if (show_residuals) do_layout = TRUE
-  if (sep_obs) do_layout = TRUE
+  if (show_residuals | sep_obs | show_errplot) do_layout = TRUE
   n_plot_rows = if (sep_obs) length(obs_vars) else 1
 
   if (do_layout) {
@@ -78,7 +81,7 @@ plot.mkinfit <- function(x, fit = x,
 
     # If the observed variables are shown separately, do row layout
     if (sep_obs) {
-      n_plot_cols = if (show_residuals) 2 else 1
+      n_plot_cols = if (show_residuals | show_errplot) 2 else 1
       n_plots = n_plot_rows * n_plot_cols
 
       # Set relative plot heights, so the first and the last plot are the norm
@@ -201,6 +204,12 @@ plot.mkinfit <- function(x, fit = x,
       }
       abline(h = 0, lty = 2)
     }
+
+    # 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)
+    }
   }
   if (do_layout) par(oldpar, no.readonly = TRUE)
 }
diff --git a/R/plot.mmkin.R b/R/plot.mmkin.R
index ee7075d3..c9d98718 100644
--- a/R/plot.mmkin.R
+++ b/R/plot.mmkin.R
@@ -1,4 +1,4 @@
-# Copyright (C) 2015-2016 Johannes Ranke
+# Copyright (C) 2015-2016,2019 Johannes Ranke
 # Contact: jranke@uni-bremen.de
 
 # This file is part of the R package mkin
@@ -16,11 +16,15 @@
 # You should have received a copy of the GNU General Public License along with
 # this program. If not, see <http://www.gnu.org/licenses/>
 
-plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", errmin_digits = 3,
+plot.mmkin <- function(x, main = "auto", legends = 1, 
+                       resplot = c("time", "errmod"),
+                       errmin_var = "All data", errmin_digits = 3,
                        cex = 0.7, rel.height.middle = 0.9, ...) {
   n.m <- nrow(x)
   n.d <- ncol(x)
 
+  resplot <- match.arg(resplot)
+
   # We can handle either a row (different models, same dataset)
   # or a column (same model, different datasets)
   if (n.m > 1 & n.d > 1) stop("Please select fits either for one model or for one dataset")
@@ -90,7 +94,11 @@ plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", e
     }
     mtext(chi2_text, cex = cex, line = 0.4)
 
-    mkinresplot(fit, legend = FALSE, ...)
+    if (resplot == "time") {
+      mkinresplot(fit, legend = FALSE, ...)
+    } else {
+      mkinerrplot(fit, legend = FALSE, ...)
+    }
     mtext(paste(fit_name, "residuals"), cex = cex, line = 0.4)
   }