diff options
Diffstat (limited to 'CakeNullPlot.R')
| -rw-r--r-- | CakeNullPlot.R | 167 |
1 files changed, 84 insertions, 83 deletions
diff --git a/CakeNullPlot.R b/CakeNullPlot.R index f5484b4..e0823f2 100644 --- a/CakeNullPlot.R +++ b/CakeNullPlot.R @@ -1,12 +1,11 @@ #$Id$ # 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 +# Some of 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 +# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414 -# This program is free software: you can redistribute it and/or modify +# The CAKE R modules are 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. @@ -17,38 +16,57 @@ # 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/>.” +# 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), sannMaxIter = 10000, ...) +CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(observed$time), + digits = max(3, getOption("digits") - 3), sannMaxIter = 10000, atol = 1e-6, ...) { # 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) + parms.all = c(state.ini, parms.ini) - # Print disappearence times + # Print disappearance 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) - extraDT50<- data.frame(DT50 = rep(NA, 2), row.names = c("k1", "k2")) + extraDT50<- data.frame(k1 = rep(NA, length(obs_vars)), k2 = 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,sannMaxIter) - } - - cat("\nEstimated disappearance times:\n") - print(distimes, digits=digits,...) + type = names(model$map[[obs_var]])[1] + if(exists("DT50")) distimes[obs_var, ] = c(DT50, DT90) + distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter) + extraDT50[obs_var, ] = CakeExtraDT(type, obs_var, parms.all) + } + + cat("\nEstimated disappearance times:\n") + print(distimes, digits=digits,...) - extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.ini) - cat("\nHalf Lives for Extra k values:\n") - print(extraDT50, digits=digits,...) - - # Plot the model + cat("\nHalf Lives for Extra k values:\n") + print(extraDT50, digits=digits,...) + + ioreRepDT = CakeIORERepresentativeDT("Parent", parms.all) + printIoreRepDT <- !is.null(ioreRepDT) + if (printIoreRepDT){ + cat("\nRepresentative Half Life for IORE:", format(signif(ioreRepDT, digits))) + } + + cat("\n") + + fomcRepDT = CakeFOMCBackCalculatedDT(parms.all) + printFomcRepDT <- !is.null(fomcRepDT) + + if (printIoreRepDT){ + cat("\nBack-calculated Half Life for FOMC:", format(signif(fomcRepDT, digits))) + } + + cat("\n") + + # Plot the model - outtimes <- seq(0, xlim[2], length.out=101) + outtimes <- seq(0, xlim[2], length.out=CAKE.plots.resolution) odeini <- state.ini odenames <- model$parms @@ -60,75 +78,58 @@ CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(obs 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))) + 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)) + + odeini <- AdjustOdeInitialValues(odeini, model, odeparms) - # Ensure initial state is at time 0 - obstimes <- unique(c(0,observed$time)) - out <- ode( - y = odeini, - times = obstimes, - func = mkindiff, - parms = odeparms, - ) + # Solve at observation times + out <- ode(y = odeini, times = obstimes, func = mkindiff, parms = odeparms, atol = atol) - # 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]]]) - } - } + out_transformed <- PostProcessOdeOutput(out, model, atol) + + predicted_long <- wide_to_long(out_transformed,time="time") - predicted_long <- mkin_wide_to_long(out_transformed,time="time") # Calculate chi2 error levels according to FOCUS (2006) - errmin<-CakeChi2(model, observed, predicted_long, obs_vars, c("auto"), c(), state.ini, c("auto")) - cat("\nChi2 error levels in percent:\n") - errmin$err.min <- 100 * errmin$err.min - print(errmin, digits=digits,...) + errmin<-CakeChi2(model, observed, predicted_long, obs_vars, c(), c(), state.ini, parms.ini, fixed) + 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) - # 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, - ) + # Solve at output times + 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]]]) - } - } + out_transformed <- PostProcessOdeOutput(out, model, atol) - cat("\n<PLOT MODEL START>\n") - print(out_transformed) - cat("<PLOT MODEL END>\n") + cat("\n<PLOT MODEL START>\n") + print(out_transformed) + cat("<PLOT MODEL END>\n") } |