summaryrefslogtreecommitdiff
path: root/CakeNullPlot.R
diff options
context:
space:
mode:
Diffstat (limited to 'CakeNullPlot.R')
-rw-r--r--CakeNullPlot.R127
1 files changed, 127 insertions, 0 deletions
diff --git a/CakeNullPlot.R b/CakeNullPlot.R
new file mode 100644
index 0000000..3eacd67
--- /dev/null
+++ b/CakeNullPlot.R
@@ -0,0 +1,127 @@
+#$Id: CakeNullPlot.R 216 2011-07-05 14:35:03Z nelr $
+# Generates model curves so the GUI can plot them. Used when all params are fixed so there was no fit
+# The CAKE R modules are based on mkin
+# Based on code in IRLSkinfit
+# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011 Syngenta
+# Author: Rob Nelson (Tessella)
+# Tessella Project Reference: 6245
+
+# This program 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/>.”
+
+CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(observed$time), digits = max(3, getOption("digits") - 3), ...)
+{
+ # Print the fixed values
+ fixed <- data.frame(value = c(state.ini, parms.ini))
+ fixed$type = c(rep("state", length(state.ini)), rep("deparm", length(parms.ini)))
+ cat("\nFixed parameter values:\n")
+ print(fixed)
+
+ # Print disappearence times
+ obs_vars = unique(as.character(observed$name))
+ distimes <- data.frame(DT50 = rep(NA, length(obs_vars)),
+ DT90 = rep(NA, length(obs_vars)), row.names = obs_vars)
+ for (obs_var in obs_vars) {
+ type = names(model$map[[obs_var]])[1]
+ if(exists("DT50")) distimes[obs_var, ] = c(DT50, DT90)
+ distimes[obs_var, ] = CakeDT(type,obs_var,parms.ini)
+ }
+ cat("\nEstimated disappearance times:\n")
+ print(distimes, digits=digits,...)
+
+ # Plot the model
+
+ outtimes <- seq(0, xlim[2], length.out=101)
+
+ odeini <- state.ini
+ odenames <- model$parms
+ odeparms <- params.ini
+
+ # Solve the system
+ evalparse <- function(string)
+ {
+ eval(parse(text=string), as.list(c(odeparms, odeini)))
+ }
+
+ mod_vars <- names(model$diffs)
+ mkindiff <- function(t, state, parms) {
+ time <- t
+ diffs <- vector()
+ for (box in mod_vars) {
+ diffname <- paste("d", box, sep = "_")
+ diffs[diffname] <- with(as.list(c(time, state, parms)),
+ eval(parse(text = model$diffs[[box]])))
+ }
+ return(list(c(diffs)))
+ }
+
+ # Ensure initial state is at time 0
+ obstimes <- unique(c(0,observed$time))
+ out <- ode(
+ y = odeini,
+ times = obstimes,
+ func = mkindiff,
+ parms = odeparms,
+ )
+
+ # Output transformation for models with unobserved compartments like SFORB
+ out_transformed <- data.frame(time = out[,"time"])
+ for (var in names(model$map)) {
+ if(length(model$map[[var]]) == 1) {
+ out_transformed[var] <- out[, var]
+ } else {
+ out_transformed[var] <- rowSums(out[, model$map[[var]]])
+ }
+ }
+
+ predicted_long <- mkin_wide_to_long(out_transformed,time="time")
+ # Calculate chi2 error levels according to FOCUS (2006)
+ errmin<-CakeChi2(observed, predicted_long, obs_vars, c(), c())
+ cat("\nChi2 error levels in percent:\n")
+ errmin$err.min <- 100 * errmin$err.min
+ print(errmin, digits=digits,...)
+
+
+ # Fit to data
+ data<-observed
+ data$err<-rep(NA,length(data$time))
+ data<-merge(data, predicted_long, by=c("time","name"))
+ names(data)<-c("time", "variable", "observed","err-var", "predicted")
+ data$residual<-data$observed-data$predicted
+ data<-data[order(data$variable,data$time),]
+ cat("\nData:\n")
+ print(format(data, digits = digits, scientific = FALSE,...), row.names = FALSE)
+
+
+ out <- ode(
+ y = odeini,
+ times = outtimes,
+ func = mkindiff,
+ parms = odeparms,
+ )
+
+ # Output transformation for models with unobserved compartments like SFORB
+ out_transformed <- data.frame(time = out[,"time"])
+ for (var in names(model$map)) {
+ if(length(model$map[[var]]) == 1) {
+ out_transformed[var] <- out[, var]
+ } else {
+ out_transformed[var] <- rowSums(out[, model$map[[var]]])
+ }
+ }
+
+ cat("\n<PLOT MODEL START>\n")
+ print(out_transformed)
+ cat("<PLOT MODEL END>\n")
+}
+

Contact - Imprint