summaryrefslogtreecommitdiff
path: root/CakeErrMin.R
diff options
context:
space:
mode:
Diffstat (limited to 'CakeErrMin.R')
-rwxr-xr-xCakeErrMin.R91
1 files changed, 91 insertions, 0 deletions
diff --git a/CakeErrMin.R b/CakeErrMin.R
new file mode 100755
index 0000000..9a074a4
--- /dev/null
+++ b/CakeErrMin.R
@@ -0,0 +1,91 @@
+## -----------------------------------------------------------------------------
+## Chi squared errmin function.
+## -----------------------------------------------------------------------------
+# Some of the CAKE R modules are based on mkin.
+#
+# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2016 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414
+#
+# The CAKE R modules are free software: you can
+# redistribute them and/or modify them 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("name"))
+
+CakeErrMin <- function(fit, alpha = 0.05)
+{
+ parms.optim <- fit$par
+ kinerrmin <- function(errdata, n.parms) {
+ means.mean <- mean(errdata$value_mean, na.rm=TRUE)
+ df = length(errdata$value_mean) - n.parms
+
+ err.min <- sqrt((1 / qchisq(1 - alpha, df)) *
+ sum((errdata$value_mean - errdata$value_pred)^2)/(means.mean^2))
+
+ return(list(err.min = err.min, n.optim = n.parms, df = df))
+ }
+
+ means <- aggregate(value ~ time + name, data = fit$observed, mean, na.rm=TRUE)
+ errdata <- merge(means, fit$predicted, by = c("time", "name"),
+ suffixes = c("_mean", "_pred"))
+ errdata <- errdata[order(errdata$time, errdata$name), ]
+
+ # Remove values at time zero for variables whose value for state.ini is fixed,
+ # as these will not have any effect in the optimization and should therefore not
+ # be counted as degrees of freedom.
+ fixed_initials = gsub("_0$", "", rownames(subset(fit$fixed, type = "state")))
+ errdata <- subset(errdata, !(time == 0 & name %in% fixed_initials))
+
+ n.optim.overall <- length(parms.optim)
+
+ errmin.overall <- kinerrmin(errdata, n.optim.overall)
+ errmin <- data.frame(err.min = errmin.overall$err.min,
+ n.optim = errmin.overall$n.optim, df = errmin.overall$df)
+ rownames(errmin) <- "All data"
+
+ for (obs_var in fit$obs_vars)
+ {
+ errdata.var <- subset(errdata, name == obs_var)
+
+ # Check if initial value is optimised
+ n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(parms.optim)))
+
+ # Rate constants and DFOP parameters are attributed to the source variable
+ n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim)))
+ n.k1.dfop.optim <- length(grep(paste("^k1", obs_var, sep="_"), names(parms.optim)))
+ n.k2.dfop.optim <- length(grep(paste("^k2", obs_var, sep="_"), names(parms.optim)))
+ n.g.dfop.optim <- length(grep(paste("^g", obs_var, sep="_"), names(parms.optim)))
+
+ # Formation fractions are attributed to the target variable
+ n.ff.optim <- length(grep(paste("^f", ".*", obs_var, "$", sep=""), names(parms.optim)))
+
+ n.optim <- n.k.optim + n.k1.dfop.optim + n.k2.dfop.optim + n.g.dfop.optim + n.initials.optim + n.ff.optim
+
+ # FOMC, HS and IORE parameters are only counted if we are looking at the
+ # first variable in the model which is always the source variable
+ if (obs_var == fit$obs_vars[[1]]) {
+ if ("alpha" %in% names(parms.optim)) n.optim <- n.optim + 1
+ if ("beta" %in% names(parms.optim)) n.optim <- n.optim + 1
+ if ("k1" %in% names(parms.optim)) n.optim <- n.optim + 1
+ if ("k2" %in% names(parms.optim)) n.optim <- n.optim + 1
+ if ("tb" %in% names(parms.optim)) n.optim <- n.optim + 1
+ if ("N" %in% names(parms.optim)) n.optim <- n.optim + 1
+ }
+
+ # Calculate and add a line to the results
+ errmin.tmp <- kinerrmin(errdata.var, n.optim)
+ errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp
+ }
+
+ return(errmin)
+}

Contact - Imprint