#$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 # 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 .” CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(observed$time), digits = max(3, getOption("digits") - 3), sannMaxIter = 10000, ...) { # 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) extraDT50<- data.frame(DT50 = rep(NA, 2), row.names = c("k1", "k2")) 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,...) extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.ini) cat("\nHalf Lives for Extra k values:\n") print(extraDT50, digits=digits,...) # Plot the model outtimes <- seq(0, xlim[2], length.out=101) odeini <- state.ini odenames <- model$parms odeparms <- parms.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(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,...) # 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\n") print(out_transformed) cat("\n") }