Version 3.2v3.2

1 files changed, 165 insertions, 196 deletions

diff --git a/CakeMcmcFit.R b/CakeMcmcFit.R
index f9a340e..65c9431 100644
--- a/CakeMcmcFit.R
+++ b/CakeMcmcFit.R
@@ -1,11 +1,9 @@
-# $Id$
-# The CAKE R modules are based on mkin, 
+# Some of the CAKE R modules are based on mkin, 
 # Based on mcmckinfit as modified by Bayer
-# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011  Syngenta
-# Author: Rob Nelson, Tamar Christina
-# Tessella Project Reference: 6245, 7247
+# 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,168 +14,155 @@
 #    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/>.
 
-CakeMcmcFit <- 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, commonsigma=FALSE,jump = NULL,prior = NULL,
-    wvar0=0.1,niter = 1000, 
-    outputlength = niter, burninlength = 0, updatecov = niter, 
-    ntrydr = 1, drscale = NULL, verbose = TRUE,fitstart=TRUE,seed=NULL,atol=1e-6,
-    sannMaxIter = 10000, control=list(),
-    useSolnp = FALSE, method='L-BFGS-B',...) 
+# Performs an Markov chain Monte Carlo least squares fit on a given CAKE model.
+#
+# 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).
+# niter: The number of MCMC iterations to apply.
+# 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 (only used for the initial first fit).
+CakeMcmcFit <- function (cake.model, 
+                         observed, 
+                         parms.ini, 
+                         state.ini, 
+                         lower, 
+                         upper, 
+                         fixed_parms = NULL, 
+                         fixed_initials, 
+                         quiet = FALSE, 
+                         niter = 1000, 
+                         verbose = TRUE, 
+                         seed=NULL, 
+                         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)
     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
         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]])))
+                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)
     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]
-    flag <- 0
     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)
-	bestIteration <- -1;
-	costWithStatus <- function(P, ...){
-		r <- costFunctions$cost(P)
-		if(r$cost == costFunctions$get.best.cost()) {
-			bestIteration<<-costFunctions$get.calls();
-			 cat(' MCMC best so far: c', r$cost, 'it:', bestIteration, '\n')
-		}
-		cat("MCMC Call no.", costFunctions$get.calls(), '\n')
-		return(r)
-	}
+    costFunctions <- CakeInternalCostFunctions(cake.model, state.ini.optim, state.ini.optim.boxnames, 
+                        state.ini.fixed, parms.fixed, observed, mkindiff, quiet, atol=atol)
+    bestIteration <- -1;
+    costWithStatus <- function(P, ...){
+        r <- costFunctions$cost(P)
+        if(r$cost == costFunctions$get.best.cost()) {
+            bestIteration<<-costFunctions$get.calls();
+             cat(' MCMC best so far: c', r$cost, 'it:', bestIteration, '\n')
+        }
+        cat("MCMC Call no.", costFunctions$get.calls(), '\n')
+        return(r)
+    }
 
     # Set the seed
     if ( is.null(seed) ) {
       # No seed so create a random one so there is something to report
       seed = runif(1,0,2^31-1)
     }
+    
     seed = as.integer(seed)
     set.seed(seed)
     
-    if(fitstart==TRUE)
-      {
-     ## ############# Get Initial Paramtervalues   #############
+    ## ############# Get Initial Paramtervalues   #############
     ## Start with no weighting
-    
-    if(useSolnp)
+    if(useExtraSolver)
     {
-        pnames=names(c(state.ini.optim, parms.optim))
-        fn <- function(P){
-            names(P) <- pnames
-            FF<<-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)
-        flag <- 1
-        optimmethod <- 'solnp'
+        a <- try(fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='Port',control=control))
+        
         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')
-            a <- try(fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='Port',control=control))
-            flag <- 0
-            
-            if(class(a) == "try-error")
-            {
-                fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method=method,control=control)
-            }
-            ## now using submethod already
+            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, 
+        # 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,...)
-	}
+    }
     
-    if(commonsigma==TRUE){
-      #observed$err <- rep(1,nrow(observed))
-      if(flag==1 && useSolnp)## fit from solnp
-      {
-          ## run a small loop with 'Marq' or some other method
-          fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, method='Port',control=control)
-      }
-      cov0 <- summary(fit)$cov.scaled*2.4^2/length(fit$par)
-	  costFunctions$set.calls(0); costFunctions$reset.best.cost()
-      res <- modMCMC(costWithStatus,fit$par,...,jump=cov0,lower=lower,upper=upper,prior=prior,var0=var0,wvar0=wvar0,niter=niter,outputlength = outputlength, burninlength = burninlength, updatecov = updatecov,ntrydr=ntrydr,drscale=drscale,verbose=verbose)
-      #res <- modMCMC(cost,fit$par,lower=lower,upper=upper,niter=niter)
-    }else{
-      ## ############## One reweighted estimation ############################ 
-      ## Estimate the error variance(sd)     
-      tmpres <- fit$residuals
-      oldERR <- observed$err
-      err <- rep(NA,length(mod_vars))
-      for(i in 1:length(mod_vars))
-        {
-          box <- mod_vars[i]
-          ind <- which(names(tmpres)==box)
-          tmp <- tmpres[ind]
-          err[i] <- sd(tmp)
-        }
-      names(err) <- mod_vars
-      ERR <- err[as.character(observed$name)]
-      observed$err <-ERR
-	  costFunctions$set.error(ERR)
-      if(flag==1 && useSolnp)## fit from solnp
-      {
-          fit$ssr <- fit$values[length(fit$values)]
-          fit$residuals <-FF$residual$res
-          ## mean square per varaible
-          if (class(FF) == "modCost") {
-              names(fit$residuals)  <- FF$residuals$name
-              fit$var_ms            <- FF$var$SSR/FF$var$N
-              fit$var_ms_unscaled   <- FF$var$SSR.unscaled/FF$var$N
-              fit$var_ms_unweighted <- FF$var$SSR.unweighted/FF$var$N
-      
-              names(fit$var_ms_unweighted) <- names(fit$var_ms_unscaled) <-
-                  names(fit$var_ms) <- FF$var$name
-          } else fit$var_ms <- fit$var_ms_unweighted <- fit$var_ms_unscaled <- NA
-      }
-      olderr <- rep(1,length(mod_vars))
-      diffsigma <- sum((err-olderr)^2)
-      ## At least do one iteration step to get a weighted LS
-      fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, ...)
-      ## Use this as the Input for MCMC algorithm
-      ## ##########################
-	  fs <- summary(fit)
-      cov0 <- if(all(is.na(fs$cov.scaled))) NULL else fs$cov.scaled*2.4^2/length(fit$par)
-      var0 <- fit$var_ms_unweighted
-	  costFunctions$set.calls(0); costFunctions$reset.best.cost()
-      res <- modMCMC(costWithStatus,fit$par,...,jump=cov0,lower=lower,upper=upper,prior=prior,var0=var0,wvar0=wvar0,niter=niter,outputlength = outputlength, burninlength = burninlength, updatecov = updatecov,ntrydr=ntrydr,drscale=drscale,verbose=verbose)
-   
+    ## ############## One reweighted estimation ############################ 
+    ## Estimate the error variance(sd)     
+    tmpres <- fit$residuals
+    oldERR <- observed$err
+    err <- rep(NA,length(cake.model$map))
+    
+    for(i in 1:length(cake.model$map)) {
+        box <- names(cake.model$map)[i]
+        ind <- which(names(tmpres)==box)
+        tmp <- tmpres[ind]
+        err[i] <- sd(tmp)
     }
-   }else {
-		stop('fitstart=FALSE code removed until needed')
-	}
+    
+    names(err) <- names(cake.model$map)
+    ERR <- err[as.character(observed$name)]
+    observed$err <-ERR
+    costFunctions$set.error(ERR)
+    
+    olderr <- rep(1,length(cake.model$map))
+    diffsigma <- sum((err-olderr)^2)
+    ## At least do one iteration step to get a weighted LS
+    fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, ...)
+    ## Use this as the Input for MCMC algorithm
+    ## ##########################
+    fs <- summary(fit)
+    cov0 <- if(all(is.na(fs$cov.scaled))) NULL else fs$cov.scaled*2.4^2/length(fit$par)
+    var0 <- fit$var_ms_unweighted
+    costFunctions$set.calls(0); costFunctions$reset.best.cost()
+    res <- modMCMC(costWithStatus,fit$par,...,jump=cov0,lower=lower,upper=upper,prior=NULL,var0=var0,wvar0=0.1,niter=niter,outputlength=niter,burninlength=0,updatecov=niter,ntrydr=1,drscale=NULL,verbose=verbose)
+  
     # Construct the fit object to return
     fit$start <- data.frame(initial = c(state.ini.optim, parms.optim))
     fit$start$type = c(rep("state", length(state.ini.optim)), 
@@ -189,92 +174,72 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
         rep("deparm", length(parms.fixed)))
 
     fit$mkindiff <- mkindiff
-    fit$map <- mkinmod$map
-    fit$diffs <- mkinmod$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)))
-    fns <- function(str) {
-      fields <- strsplit(str, "\\.")[[1]]
-      if (fields[1] == "value") {
-         return(fields[2])
-      }
-      else {
-         return(str)
-      }
-    }
-    names(fit$observed) <- sapply(names(fit$observed), fns)
+    fit$map <- cake.model$map
+    fit$diffs <- cake.model$diffs
    
-   # Replace mean from modFit with mean from modMCMC
-   fnm <- function(x) mean(res$pars[,x])
-   fit$par <- sapply(dimnames(res$pars)[[2]],fnm)
-   fit$bestpar <- res$bestpar
-   fit$costfn <- costWithStatus
+    # Replace mean from modFit with mean from modMCMC
+    fnm <- function(x) mean(res$pars[,x])
+    fit$par <- sapply(dimnames(res$pars)[[2]],fnm)
+    fit$bestpar <- res$bestpar
+    fit$costfn <- costWithStatus
    
-   # Disappearence times
-   parms.all <- c(fit$par, parms.fixed)
-   obs_vars = unique(as.character(observed$name))
-   fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), 
+    # Disappearence times
+    parms.all <- c(fit$par, parms.fixed)
+    obs_vars = unique(as.character(observed$name))
+    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"))     
-    for (obs_var in obs_vars) {
-        type = names(mkinmod$map[[obs_var]])[1]
-        fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter)
+    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 (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)
     
-   # Fits to data
-   odeini <- state.ini
-   odeparms <- parms.all
-   # If this step modelled the start amount, include it. Lookup in vectors returns NA not null if missing
-   if(!is.na(odeparms["Parent_0"])) { odeini['Parent'] = odeparms["Parent_0"] }
+    fit$ioreRepDT = CakeIORERepresentativeDT("Parent", parms.all)
+    fit$fomcRepDT = CakeFOMCBackCalculatedDT(parms.all)
+    
+    # Fits to data
+    odeini <- state.ini
+    odeparms <- parms.all
+    # If this step modelled the start amount, include it. Lookup in vectors returns NA not null if missing
+    if(!is.na(odeparms["Parent_0"])) { odeini['Parent'] = odeparms["Parent_0"] }
    
-   # Solve the system
-   evalparse <- function(string)
-   {
-      eval(parse(text=string), as.list(c(odeparms, odeini)))
-   }
+    # Solve the system
+    evalparse <- function(string)
+    {
+       eval(parse(text=string), as.list(c(odeparms, odeini)))
+    }
   
-   # Ensure initial state is at time 0
-   obstimes <- unique(c(0,observed$time))
+    # Ensure initial state is at time 0
+    obstimes <- unique(c(0,observed$time))
 
-   out <- ode(
-      y = odeini,
-      times = obstimes,
-      func = mkindiff, 
-      parms = odeparms,
-    )
-    
-   # Output transformation for models with unobserved compartments like SFORB
-   out_transformed <- data.frame(time = out[,"time"])
-   for (var in names(mkinmod$map)) {
-      if(length(mkinmod$map[[var]]) == 1) {
-         out_transformed[var] <- out[, var]
-      } else {
-         out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]])
-      }
-   }
-   fit$predicted <- out_transformed
-   fit$penalties <- CakePenaltiesLong(odeparms, out_transformed, observed)
+    odeini <- AdjustOdeInitialValues(odeini, cake.model, odeparms)
+   
+    out <- ode(y = odeini, times = obstimes, func = mkindiff, parms = odeparms, atol = atol)
+   
+    out_transformed <- PostProcessOdeOutput(out, cake.model, atol)
+   
+    fit$predicted <- out_transformed
+    fit$penalties <- CakePenaltiesLong(odeparms, out_transformed, observed)
 
-   predicted_long <- mkin_wide_to_long(out_transformed,time="time")
-   obs_vars = unique(as.character(observed$name))
-   fit$errmin <- CakeChi2(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini)
+    predicted_long <- wide_to_long(out_transformed,time="time")
+    obs_vars = unique(as.character(observed$name))
+    fit$errmin <- CakeChi2(cake.model, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini, fit$fixed)
 
-   data<-observed
-   data$err<-rep(NA,length(data$time))
-   data<-merge(data, 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), ]
+    data<-observed
+    data$err<-rep(NA,length(data$time))
+    data<-merge(data, 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), ]
+    fit$atol <- atol
 
-   sq <- data$residual * data$residual
-   fit$ssr <- sum(sq)
+    sq <- data$residual * data$residual
+    fit$ssr <- sum(sq)
    
-   fit$seed = seed
+    fit$seed = seed
    
     fit$res <- res
     class(fit) <- c("CakeMcmcFit", "mkinfit", "modFit")
@@ -355,10 +320,14 @@ summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, halflives
   ans$start <- object$start
   ans$fixed <- object$fixed
   ans$errmin <- object$errmin 
-  if(distimes) ans$distimes <- object$distimes
+  if(distimes){ 
+    ans$distimes <- object$distimes
+    ans$extraDT50 <- object$extraDT50
+    ans$ioreRepDT <- object$ioreRepDT
+    ans$fomcRepDT <- object$fomcRepDT
+  }
   if(halflives) ans$halflives <- object$halflives
   if(ff) ans$ff <- object$ff
-  if(length(object$SFORB) != 0) ans$SFORB <- object$SFORB
   class(ans) <- c("summary.CakeFit","summary.mkinfit", "summary.modFit") 
   return(ans)  
 }