diff options
Diffstat (limited to 'CakeIrlsFit.R')
| -rw-r--r-- | CakeIrlsFit.R | 147 |
1 files changed, 147 insertions, 0 deletions
diff --git a/CakeIrlsFit.R b/CakeIrlsFit.R new file mode 100644 index 0000000..bf20572 --- /dev/null +++ b/CakeIrlsFit.R @@ -0,0 +1,147 @@ +#$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) +} |