Version 3.2v3.2

1 files changed, 101 insertions, 74 deletions

diff --git a/CakeIrlsFit.R b/CakeIrlsFit.R
index 05eff0a..92c45da 100644
--- a/CakeIrlsFit.R
+++ b/CakeIrlsFit.R
@@ -1,11 +1,10 @@
 #$Id$
 #
-# The CAKE R modules are based on mkin
-# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011 Syngenta
-# Authors: Rob Nelson, Richard Smith, Tamar Christina
-# Tessella Project Reference: 6245, 7247
+# Some of the CAKE R modules are based on mkin
+# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414
  
-#    This program is free software: you can redistribute it and/or modify
+#    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.
@@ -16,31 +15,52 @@
 #    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/>.”
+#    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, sannMaxIter = 10000, control=list(),
-    useSolnp = FALSE, method='L-BFGS-B',...) 
-{
 
-### 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
+# Performs an iteratively-reweighted least squares fit on a given CAKE model.
+# Remark: this function was originally based on the "mkinfit" function, version 0.1.
+#
+# 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.
+# sannMaxIter: The maximum number of iterations to apply to SANN processes.
+# control: ...
+# useExtraSolver: Whether to use the extra solver for this fit.
+CakeIrlsFit <- 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, 
+                         sannMaxIter = 10000, 
+                         control=list(),
+                         useExtraSolver = FALSE,
+                         ...) 
+{
     NAind <-which(is.na(observed$value))
-    mod_vars <- names(mkinmod$diffs)
-    observed <- subset(observed, name %in% names(mkinmod$map))
+    mod_vars <- names(cake.model$diffs)
+    observed <- subset(observed, name %in% names(cake.model$map))
     ERR <- rep(1,nrow(observed))
     observed <- cbind(observed,err=ERR)
-    flag <- 0
+    fitted_with_extra_solver <- 0
     
     obs_vars = unique(as.character(observed$name))
-    if (is.null(names(parms.ini))) 
-        names(parms.ini) <- mkinmod$parms
+    
+    if (is.null(names(parms.ini))) {
+        names(parms.ini) <- cake.model$parms
+    }
         
     mkindiff <- function(t, state, parms) {
         time <- t
@@ -48,12 +68,14 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
         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]])))
+                                    eval(parse(text = cake.model$diffs[[box]])))
         }
         return(list(c(diffs)))
     }
-    if (is.null(names(state.ini))) 
+    
+    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)
@@ -67,15 +89,15 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
         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)
-	
+    
+    costFunctions <- CakeInternalCostFunctions(cake.model, state.ini.optim, state.ini.optim.boxnames, 
+            state.ini.fixed, parms.fixed, observed, mkindiff, quiet, atol=atol)
+    
     ############### Iteratively Reweighted Least Squares#############
     ## Start with no weighting
     
     ## a prefitting step since this is usually the most effective method
-    if(useSolnp)
+    if(useExtraSolver)
     {
         pnames=names(c(state.ini.optim, parms.optim))
         fn <- function(P){
@@ -83,28 +105,28 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
             FF<<-costFunctions$cost(P)
             return(FF$model)}
         a <- try(fit <- solnp(c(state.ini.optim, parms.optim),fun=fn,LB=lower,UB=upper,control=control),silent=TRUE)
-        #optimmethod <- method0
 
-        flag <- 1
+        fitted_with_extra_solver <- 1
         if(class(a) == "try-error")
         {
             print('solnp fails, try PORT or other algorithm by users choice, might take longer time. Do something else!')
             warning('solnp fails, switch to  PORT or other algorithm by users choice')
             ## now using submethod already
             a <- try(fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='Port',control=control))
-            flag <- 0
+            fitted_with_extra_solver <- 0
             if(class(a) == "try-error")
             {
-                fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method=method,control=control)
+                fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='L-BFGS-B',control=control)
             }
         }
-    }else{        
-        fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, 
+    }else{
+        # 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,control=control,...)
     }
     
-    ## print(fit$hessian)
-	
     if(length(control)==0)
       {
         irls.control <- list(maxIter=10,tol=1e-05)
@@ -122,15 +144,15 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     tol <- irls.control$tol
     ####
     
-    if(length(mkinmod$map)==1 || useSolnp){
+    if(length(cake.model$map)==1){
         ## there is only one parent just do one iteration:
         maxIter <- 0
         
-        if(flag==1)## fit from solnp
+        if(fitted_with_extra_solver==1)## managed to fit with extra solver
         {
             fit$ssr <- fit$values[length(fit$values)]
             fit$residuals <-FF$residual$res
-            ## mean square per varaible
+            ## mean square per variable
             if (class(FF) == "modCost") {
                 names(fit$residuals)  <- FF$residuals$name
                 fit$var_ms            <- FF$var$SSR/FF$var$N
@@ -149,14 +171,14 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     niter <- 1
     ## insure one IRLS iteration
     diffsigma <- 100
-    olderr <- rep(1,length(mod_vars))
+    olderr <- rep(1,length(cake.model$map))
     while(diffsigma>tol & niter<=maxIter)
-      {      
-        if(flag==1 && useSolnp)## fit from solnp
+    {      
+        if(fitted_with_extra_solver==1 && useExtraSolver)## managed to fit with extra solver
         {
             fit$ssr <- fit$values[length(fit$values)]
             fit$residuals <-FF$residual$res
-            ## mean square per varaible
+            ## mean square per variable
             if (class(FF) == "modCost") {
                 names(fit$residuals)  <- FF$residuals$name
                 fit$var_ms            <- FF$var$SSR/FF$var$N
@@ -168,7 +190,6 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
             } else fit$var_ms <- fit$var_ms_unweighted <- fit$var_ms_unscaled <- NA
         }
 
-        
         err <- sqrt(fit$var_ms_unweighted)
         ERR <- err[as.character(observed$name)]
         costFunctions$set.error(ERR)
@@ -176,31 +197,32 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
         cat("IRLS iteration at",niter, "; Diff in error variance ", diffsigma,"\n")
         olderr <- err
         
-        if(useSolnp)
+        if(useExtraSolver)
         {
-            flag <- 1
+            fitted_with_extra_solver <- 1
             a <- try(fit <- solnp(fit$par,fun=fn,LB=lower,UB=upper,control=control),silent=TRUE)
 
             if(class(a) == "try-error")
             {
-                flag <- 0
+                fitted_with_extra_solver <- 0
                 print('solnp fails during IRLS iteration, try PORT or other algorithm by users choice.This may takes a while. Do something else!') ## NOTE: because in kingui we switch off the warnings, we need to print out the message instead.
                 warning('solnp fails during IRLS iteration, switch to  PORT or other algorithm by users choice')
             
                 fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, method='Port',control=list())
             }
         }else{
+            # 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.
+            fit$par <- ShiftAwayFromBoundaries(fit$par, lower, upper)
+            
             fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...)
         }
-        
-        ## print(fit$hessian)
                 
         niter <- niter+1       
        
         ### If not converged, reweight and fit                
-      }
+    }
       
-    if(flag==1 && useSolnp){
+    if(fitted_with_extra_solver==1 && useExtraSolver){
         ## solnp used
         optimmethod <- 'solnp'
         fit$ssr <- fit$values[length(fit$values)]
@@ -218,47 +240,50 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
         np <- length(c(state.ini.optim, parms.optim))
         fit$rank <- np
         fit$df.residual <- length(fit$residuals) - fit$rank
+        
+        # solnp can return an incorrect Hessian, so we use another fitting method at the optimised point to determine the Hessian
+        fitForHessian <- modFit(costFunctions$cost, fit$par, lower=lower, upper=upper, method='L-BFGS-B', control=list())
+        
+        fit$solnpHessian <- fit$hessian
+        fit$hessian <- fitForHessian$hessian
     }
-	  
-	###########################################
+      
+    ###########################################
     fit$mkindiff <- mkindiff
-    fit$map <- mkinmod$map
-    fit$diffs <- mkinmod$diffs
-    
-	out_predicted <- costFunctions$get.predicted()
+    fit$map <- cake.model$map
+    fit$diffs <- cake.model$diffs
     
-    # 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)))
+    out_predicted <- costFunctions$get.predicted()
     
-    predicted_long <- mkin_wide_to_long(out_predicted, time = "time")
+    predicted_long <- 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)))
+    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)))
+    rep("deparm", length(parms.fixed)))
 
-    fit$errmin <- CakeChi2(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini)
-    parms.all = c(fit$par, parms.fixed)
-	fit$penalties <- CakePenaltiesLong(parms.all, out_predicted, observed)
+    fit$errmin <- CakeChi2(cake.model, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini, fit$fixed)
+    parms.all = c(fit$par, parms.fixed, state.ini)
+    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)
-    fit$extraDT50<- data.frame(DT50 = rep(NA, 2), row.names = c("k1", "k2"))   
+    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 (obs_var in obs_vars) {
-        type = names(mkinmod$map[[obs_var]])[1]
-        fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter)
+    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,sannMaxIter)
+      fit$extraDT50[compartment.name, ] = CakeExtraDT(type, compartment.name, parms.all)
     }
 
-    fit$extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.all)
+    fit$ioreRepDT = CakeIORERepresentativeDT("Parent", parms.all)
+    fit$fomcRepDT = CakeFOMCBackCalculatedDT(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)
@@ -268,6 +293,8 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     fit$upper <- upper
     fit$lower <- lower
     
+    fit$atol <- atol
+    
     class(fit) <- c("CakeFit", "mkinfit", "modFit") 
     return(fit)
 }