# Originally: mkinfit.R 87 2010-12-09 07:31:59Z jranke $
# Based on code in mkinfit
# Portions Johannes Ranke 2010
# Contact: mkin-devel@lists.berlios.de
# The summary function is an adapted and extended version of summary.modFit
# from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn
# inspired by summary.nls.lm
# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta
# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
#
# 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.
#
# 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/>.
# Performs an ordinary least squares fit on a given CAKE model.
#
# cake.model: The model to perform the fit on (as generated by CakeModel.R).
# observed: Observation data to fit to.
# parms.ini: Initial values for the parameters being fitted.
# state.ini: Initial state (i.e. initial values for concentration, the dependent variable being modelled).
# lower: Lower bounds to apply to parameters.
# upper: Upper bound to apply to parameters.
# fixed_parms: A vector of names of parameters that are fixed to their initial values.
# fixed_initials: A vector of compartments with fixed initial concentrations.
# quiet: Whether the internal cost functions should execute more quietly than normal (less output).
# atol: The tolerance to apply to the ODE solver.
# dfopDtMaxIter: The maximum number of iterations to apply to DFOP DT calculation.
CakeOlsFit <- function(cake.model,
observed,
parms.ini = rep(0.1, length(cake.model$parms)),
state.ini = c(100, rep(0, length(cake.model$diffs) - 1)),
lower = 0,
upper = Inf,
fixed_parms = NULL,
fixed_initials = names(cake.model$diffs)[-1],
quiet = FALSE,
atol = 1e-6,
dfopDtMaxIter = 10000,
...)
{
mod_vars <- names(cake.model$diffs)
# Subset dataframe with mapped (modelled) variables
observed <- subset(observed, name %in% names(cake.model$map))
# Get names of observed variables
obs_vars <- unique(as.character(observed$name))
# Name the parameters if they are not named yet
if(is.null(names(parms.ini))) names(parms.ini) <- cake.model$parms
# Name the inital parameter values if they are not named yet
if(is.null(names(state.ini))) names(state.ini) <- mod_vars
# Parameters to be optimised
parms.fixed <- parms.ini[fixed_parms]
optim_parms <- setdiff(names(parms.ini), fixed_parms)
parms.optim <- parms.ini[optim_parms]
state.ini.fixed <- state.ini[fixed_initials]
optim_initials <- setdiff(names(state.ini), fixed_initials)
state.ini.optim <- state.ini[optim_initials]
state.ini.optim.boxnames <- names(state.ini.optim)
if(length(state.ini.optim) > 0) {
names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
}
# Decide if the solution of the model can be based on a simple analytical
# formula, the spectral decomposition of the matrix (fundamental system)
# or a numeric ode solver from the deSolve package
if (length(cake.model$map) == 1) {
solution <- "analytical"
} else {
solution <- "deSolve"
}
# Create a function calculating the differentials specified by the model
# if necessary
if(solution == "deSolve") {
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=cake.model$diffs[[box]])))
}
return(list(c(diffs)))
}
}
costFunctions <- CakeInternalCostFunctions(cake.model, state.ini.optim, state.ini.optim.boxnames, state.ini.fixed,
parms.fixed, observed, mkindiff = mkindiff, quiet, atol=atol, solution=solution, err=NULL)
# modFit parameter transformations can explode if you put in parameters that are equal to a bound, so we move them away by a tiny amount.
all.optim <- ShiftAwayFromBoundaries(c(state.ini.optim, parms.optim), lower, upper)
fit <-modFit(costFunctions$cost, all.optim, lower = lower, upper = upper, ...)
# We need to return some more data for summary and plotting
fit$solution <- solution
if (solution == "deSolve") {
fit$mkindiff <- mkindiff
}
# We also need various other information for summary and plotting
fit$map <- cake.model$map
fit$diffs <- cake.model$diffs
# Collect initial parameter values in two dataframes
fit$start <- data.frame(initial = c(state.ini.optim, parms.optim))
fit$start$type <- c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim)))
fit$start$lower <- lower
fit$start$upper <- upper
fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed))
fit$fixed$type <- c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed)))
# Ensure initial state is at time 0
obstimes <- unique(c(0,observed$time))
out_predicted <- CakeOdeSolve(fit, obstimes, solution = solution, atol)
out_transformed <- PostProcessOdeOutput(out_predicted, cake.model, atol)
predicted_long <- wide_to_long(out_transformed, time = "time")
fit$predicted <- out_transformed
# Calculate chi2 error levels according to FOCUS (2006)
fit$errmin <- CakeChi2(cake.model, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini, fit$fixed)
# Calculate dissipation times DT50 and DT90 and formation fractions
parms.all <- c(fit$par, parms.fixed)
fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)),
row.names = obs_vars)
fit$extraDT50<- data.frame(k1 = rep(NA, length(names(cake.model$map))), k2 = rep(NA, length(names(cake.model$map))), row.names = names(cake.model$map))
for (compartment.name in names(cake.model$map)) {
type <- names(cake.model$map[[compartment.name]])[1]
fit$distimes[compartment.name, ] <- CakeDT(type,compartment.name,parms.all,dfopDtMaxIter)
fit$extraDT50[compartment.name, ] <- CakeExtraDT(type, compartment.name, parms.all)
}
fit$ioreRepDT <- CakeIORERepresentativeDT("Parent", parms.all)
fit$fomcRepDT <- CakeFOMCBackCalculatedDT(parms.all)
fit$penalties <- CakePenaltiesLong(parms.all, out_transformed, observed)
# Collect observed, predicted and residuals
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$variable <- ordered(data$variable, levels = obs_vars)
fit$data <- data[order(data$variable, data$time), ]
fit$atol <- atol
fit$topology <- cake.model$topology
class(fit) <- c("CakeFit", "mkinfit", "modFit")
return(fit)
}
