summaryrefslogtreecommitdiff
path: root/CakeIrlsFit.R
diff options
context:
space:
mode:
Diffstat (limited to 'CakeIrlsFit.R')
-rw-r--r--CakeIrlsFit.R147
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)
+}

Contact - Imprint