Version 1.4v1.4

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)
+}