#$Id: CakeIrlsFit.R 216 2011-07-05 14:35:03Z nelr $
#
# The CAKE R modules are based on mkin
# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011 Syngenta
# Authors: Rob Nelson, Richard Smith
# 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/>.”
#
CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$parms)),
state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), lower = 0,
upper = Inf, fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1],
plot = FALSE, quiet = FALSE, err = NULL, weight = "none",
scaleVar = FALSE, atol=1e-6, control=list(),...)
{
### This is a modification based on the "mkinfit" function.
### version 0.1 July 20
###
# This version has been modified to expect SFO parameterised as k and flow fractions
# Based on code in IRLSkinfit
NAind <-which(is.na(observed$value))
mod_vars <- names(mkinmod$diffs)
observed <- subset(observed, name %in% names(mkinmod$map))
ERR <- rep(1,nrow(observed))
observed <- cbind(observed,err=ERR)
obs_vars = unique(as.character(observed$name))
if (is.null(names(parms.ini)))
names(parms.ini) <- mkinmod$parms
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 = mkinmod$diffs[[box]])))
}
return(list(c(diffs)))
}
if (is.null(names(state.ini)))
names(state.ini) <- mod_vars
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 = "_")
}
costFunctions <- CakeInternalCostFunctions(mkinmod, state.ini.optim, state.ini.optim.boxnames,
state.ini.fixed, parms.fixed, observed, mkindiff, scaleVar, quiet, atol=atol)
############### Iteratively Reweighted Least Squares#############
## Start with no weighting
fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower,
upper = upper,control=control,...)
if(length(control)==0)
{
irls.control <- list(maxIter=10,tol=1e-05)
control <- list(irls.control=irls.control)
}else{
if(is.null(control$irls.control))
{
irls.control <- list(maxIter=10,tol=1e-05)
control <- list(irls.control=irls.control)
}
}
irls.control <- control$irls.control
maxIter <- irls.control$maxIter
tol <- irls.control$tol
####
niter <- 1
## insure one IRLS iteration
diffsigma <- 100
olderr <- rep(1,length(mod_vars))
while(diffsigma>tol & niter<=maxIter)
{
err <- sqrt(fit$var_ms_unweighted)
ERR <- err[as.character(observed$name)]
costFunctions$set.error(ERR)
diffsigma <- sum((err-olderr)^2)
cat("IRLS iteration at",niter, "; Diff in error variance ", diffsigma,"\n")
olderr <- err
fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...)
niter <- niter+1
### If not converged, reweight and fit
}
###########################################
fit$mkindiff <- mkindiff
fit$map <- mkinmod$map
fit$diffs <- mkinmod$diffs
out_predicted <- costFunctions$get.predicted()
# mkin_long_to_wide does not handle ragged data
fit$observed <- reshape(observed, direction="wide", timevar="name", idvar="time")
names(fit$observed) <- c("time", as.vector(unique(observed$name)))
predicted_long <- mkin_wide_to_long(out_predicted, time = "time")
fit$predicted <- out_predicted
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)))
fit$errmin <- CakeChi2(observed, predicted_long, obs_vars, parms.optim, state.ini.optim)
parms.all = c(fit$par, parms.fixed)
fit$penalties <- CakePenaltiesLong(parms.all, out_predicted, observed)
fit$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(mkinmod$map[[obs_var]])[1]
fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all)
}
data <- merge(observed, 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), ]
class(fit) <- c("CakeFit", "mkinfit", "modFit")
return(fit)
}
