1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
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])
}
}
|
