diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/mkinerrplot.R | 84 | ||||
| -rw-r--r-- | R/mkinresplot.R | 14 | ||||
| -rw-r--r-- | R/plot.mkinfit.R | 19 | ||||
| -rw-r--r-- | R/plot.mmkin.R | 14 |
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) } |