Version 1.4v1.4

1 files changed, 301 insertions, 0 deletions

diff --git a/CakeSummary.R b/CakeSummary.R
new file mode 100644
index 0000000..ec42e75
--- /dev/null
+++ b/CakeSummary.R
@@ -0,0 +1,301 @@
+# $Id: CakeSummary.R 216 2011-07-05 14:35:03Z nelr $
+# CakeSummary: Functions to calculate statistics used in the summary,
+# and display the summary itself.
+
+# Parts modified from package mkin
+# Parts Tessella (Rob Nelson/Richard Smith) for Syngenta: Copyright (C) 2011  Syngenta
+# 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/>.”
+
+# Calculate the decay times for a compartment
+# Arguments
+#  type - type of compartment
+#  obs_var - compartment name
+#  parms.all - list of parameters
+CakeDT<-function(type, obs_var, parms.all) {
+    if (type == "SFO") {
+      # Changed to new parameterisation
+      #k_names = grep(paste("k", obs_var, sep="_"), names(parms.all), value=TRUE)
+      k_name = paste("k", obs_var, sep="_")
+      k_tot = parms.all[k_name]
+      DT50 = log(2)/k_tot
+      DT90 = log(10)/k_tot
+    }
+    if (type == "FOMC") {
+      alpha = parms.all["alpha"]
+      beta = parms.all["beta"]
+      DT50 = beta * (2^(1/alpha) - 1)
+      DT90 = beta * (10^(1/alpha) - 1)
+    }
+    if (type == "DFOP") {
+      k1 = parms.all["k1"]
+      k2 = parms.all["k2"]
+      g = parms.all["g"]
+      f <- function(t, x) {
+        fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
+        (fraction - (1 - x/100))^2
+      }
+      DTmax <- 1000
+      DT50.o <- optimize(f, c(0.001, DTmax), x=50)$minimum
+      DT50 = ifelse(DTmax - DT50.o < 0.1, NA, DT50.o)
+      DT90.o <- optimize(f, c(0.001, DTmax), x=90)$minimum
+      DT90 = ifelse(DTmax - DT90.o < 0.1, NA, DT90.o)
+    }
+    if (type == "HS") {
+      k1 = parms.all["k1"]
+      k2 = parms.all["k2"]
+      tb = parms.all["tb"]
+      DTx <- function(x) {
+        DTx.a <- (log(100/(100 - x)))/k1
+        DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2
+        if (DTx.a < tb) DTx <- DTx.a
+        else DTx <- DTx.b
+        return(DTx)
+      }
+      DT50 <- DTx(50)
+      DT90 <- DTx(90)
+    }
+    if (type == "SFORB") {
+      # FOCUS kinetics (2006), p. 60 f
+      k_out_names = grep(paste("k", obs_var, "free", sep="_"), names(parms.all), value=TRUE)
+      k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
+      k_1output = sum(parms.all[k_out_names])
+      k_12 = parms.all[paste("k", obs_var, "free", "bound", sep="_")]
+      k_21 = parms.all[paste("k", obs_var, "bound", "free", sep="_")]
+
+      sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
+      b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
+      b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
+
+      SFORB_fraction = function(t) {
+        ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
+        ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)
+      }
+      f_50 <- function(t) (SFORB_fraction(t) - 0.5)^2
+      max_DT <- 1000
+      DT50.o <- optimize(f_50, c(0.01, max_DT))$minimum
+      if (abs(DT50.o - max_DT) < 0.01) DT50 = NA else DT50 = DT50.o
+      f_90 <- function(t) (SFORB_fraction(t) - 0.1)^2
+      DT90.o <- optimize(f_90, c(0.01, 1000))$minimum
+      if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o
+    }
+    ret<-c(DT50, DT90)
+    names(ret)<-c("DT50","DT90")
+    ret
+}
+
+# Compute chi2 etc
+# Arguments
+#  observed - data.frame of observed data
+#  predicted_long - fitted values etc
+#  obs_vars - compartment names
+#  parms.optim - list of fitted parameters
+#  state,ini.optim - list of fitted initial values
+#
+CakeChi2<-function(observed, predicted_long, obs_vars, parms.optim, state.ini.optim) {
+  # Calculate chi2 error levels according to FOCUS (2006)
+  means <- aggregate(value ~ time + name, data = observed, mean, na.rm=TRUE)
+  errdata <- merge(means, predicted_long, by = c("time", "name"), suffixes = c("_mean", "_pred"))
+  errdata <- errdata[order(errdata$time, errdata$name), ]
+  errmin.overall <- mkinerrmin(errdata, length(parms.optim) + length(state.ini.optim))
+
+  errmin <- data.frame(err.min = errmin.overall$err.min, 
+    n.optim = errmin.overall$n.optim, df = errmin.overall$df)
+  rownames(errmin) <- "All data"
+  for (obs_var in obs_vars)
+  {
+    errdata.var <- subset(errdata, name == obs_var)
+    n.k.optim <- length(grep(paste("k", obs_var, sep="_"), names(parms.optim)))
+    n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(state.ini.optim)))
+    n.optim <- n.k.optim + n.initials.optim
+    if ("alpha" %in% names(parms.optim)) n.optim <- n.optim + 1
+    if ("beta" %in% names(parms.optim)) n.optim <- n.optim + 1
+    if ("k1" %in% names(parms.optim)) n.optim <- n.optim + 1
+    if ("k2" %in% names(parms.optim)) n.optim <- n.optim + 1
+    if ("g" %in% names(parms.optim)) n.optim <- n.optim + 1
+    if ("tb" %in% names(parms.optim)) n.optim <- n.optim + 1
+    errmin.tmp <- mkinerrmin(errdata.var, n.optim)
+    errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp
+  }
+  errmin
+}
+
+# Compute additional statistics (r²/efficiency) of obs v pred for each compartment and the whole problem
+CakeAdditionalStats<-function(obsvspred){
+	agg<-aggregate(obsvspred[c('observed', 'err-var', 'predicted', 'residual')], list(name=obsvspred$variable), function(x){x}, simplify=FALSE)
+	frames<-apply(agg, 1, as.data.frame);
+	names(frames)<-agg$name;
+	frames$all<-obsvspred
+	
+	t(sapply(frames, function(frame){
+		# This function takes a data frame for a single variable and makes the stats for it
+		# r²: FOCUS eq 6.5 p. 97
+		# Efficiency: eq 6.6 p. 99
+		do <- frame$observed - mean(frame$observed)
+		dp <- frame$predicted - mean(frame$predicted)
+		r2 <- ( sum(do * dp) / sqrt(sum(do^2) * sum(dp^2)) ) ^ 2
+		eff <- 1 - ( sum((frame$predicted - frame$observed)^2) / sum(do^2) )
+		list(r2=r2, eff=eff)
+	}))
+}
+
+# Summarise a fit - used for CakeIrlsFit and CakeOlsFit
+summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov = FALSE,...) {
+  param  <- object$par
+  pnames <- names(param)
+  p      <- length(param)
+  covar  <- try(solve(0.5*object$hessian), silent = TRUE)   # unscaled covariance
+
+  rdf    <- object$df.residual
+  resvar <- object$ssr / rdf
+  modVariance <- object$ssr / length(object$residuals)
+  
+  if (!is.numeric(covar)) {
+#    message <- "Cannot estimate covariance; system is singular"
+#    warning(message)
+#    covar <- matrix(data = NA, nrow = p, ncol = p)
+   covar=NULL
+  } else {
+    message <- "ok"
+    rownames(covar) <- colnames(covar) <-pnames
+    se     <- sqrt(diag(covar) * resvar)
+    names(se) <- pnames
+    tval      <- param / se
+  }
+
+  if (!all(object$start$lower >=0)) {
+    message <- "Note that the one-sided t-test may not be appropriate if
+      parameter values below zero are possible."
+    warning(message)
+  } else message <- "ok"
+
+  if(cov && !is.null(covar)) {
+    alpha <- 0.05
+    lci <- param + qt(alpha/2,rdf)*se
+    uci <- param + qt(1-alpha/2,rdf)*se
+    param <- cbind(param, se, tval, pt(tval, rdf, lower.tail = FALSE), lci, uci)
+    dimnames(param) <- list(pnames, c("Estimate", "Std. Error",
+                                    "t value", "Pr(>t)", "Lower CI", "Upper CI"))
+    ans <- list(residuals = object$residuals,
+                residualVariance = resvar,
+                sigma = sqrt(resvar),
+                modVariance = modVariance,
+                df = c(p, rdf), cov.unscaled = covar,
+                cov.scaled = covar * resvar,
+                info = object$info, niter = object$iterations,
+                stopmess = message,
+                par = param)
+  } else {
+    param <- cbind(param)  
+    ans <- list(residuals = object$residuals,
+                residualVariance = resvar,
+                sigma = sqrt(resvar),
+                modVariance = modVariance,
+                df = c(p, rdf),
+                info = object$info, niter = object$iterations,
+                stopmess = message,
+                par = param)
+    }
+
+  ans$diffs <- object$diffs
+  if(data) {
+	ans$data <- object$data
+	ans$additionalstats = CakeAdditionalStats(object$data)
+  }
+  ans$start <- object$start
+  ans$fixed <- object$fixed
+  ans$errmin <- object$errmin 
+  ans$penalties <- object$penalties
+  if(distimes) ans$distimes <- object$distimes
+  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)  
+}
+
+# Print a summary. Used for CakeIrlsFit, CakeOlsFit and CakeMcmcFit
+# Expanded from print.summary.modFit
+print.summary.CakeFit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
+  cat("\nEquations:\n")
+  print(noquote(as.character(x[["diffs"]])))
+  df  <- x$df
+  rdf <- df[2]
+
+  cat("\nStarting values for optimised parameters:\n")
+  print(x$start)
+
+  cat("\nFixed parameter values:\n")
+  if(length(x$fixed$value) == 0) cat("None\n")
+  else print(x$fixed)
+  
+  if ( !is.null(x$par) ) {
+    cat("\nOptimised parameters:\n")
+    printCoefmat(x$par, digits = digits, ...)
+  }
+
+  cat("\nResidual standard error:",
+      format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")
+
+  cat("\nChi2 error levels in percent:\n")
+  x$errmin$err.min <- 100 * x$errmin$err.min
+  print(x$errmin, digits=digits,...)
+
+  printdistimes <- !is.null(x$distimes)
+  if(printdistimes){
+    cat("\nEstimated disappearance times:\n")
+    print(x$distimes, digits=digits,...)
+  }    
+
+  printff <- !is.null(x$ff)
+  if(printff){
+    cat("\nEstimated formation fractions:\n")
+    print(data.frame(ff = x$ff), digits=digits,...)
+  }    
+
+  printSFORB <- !is.null(x$SFORB)
+  if(printSFORB){
+    cat("\nEstimated Eigenvalues of SFORB model(s):\n")
+    print(x$SFORB, digits=digits,...)
+  }    
+
+  printcor <- !is.null(x$cov.unscaled)
+  if (printcor){
+    Corr <- cov2cor(x$cov.unscaled)
+    rownames(Corr) <- colnames(Corr) <- rownames(x$par)
+    cat("\nParameter correlation:\n")
+    print(Corr, digits = digits, ...)
+  }
+
+  printdata <- !is.null(x$data)
+  if (printdata){
+    cat("\nData:\n")
+    print(format(x$data, digits = digits, scientific = FALSE,...), row.names = FALSE)
+  }
+  
+  if(!is.null(x$additionalstats)){
+	cat("\nAdditional Stats:\n");
+	print(x$additionalstats, digits=digits, ...)
+  }
+  
+  if(!(is.null(x$penalties) || 0 == dim(x$penalties)[[1]])){
+	cat("\nPenalties:\n");
+	print(x$penalties, digits=digits, row.names = FALSE, ...)
+  }
+  
+  if ( !is.null(x$seed) ) {
+      cat("\nRandom Seed:", x$seed, "\n")
+  }
+  invisible(x)
+}