summaryrefslogtreecommitdiff
path: root/CakeSummary.R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2017-10-18 11:18:07 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2017-10-18 11:18:07 +0200
commit3d6b4b4b8293a4a4ab6f06805e1380600373796c (patch)
tree4ec4dd041427b0154684122f5ae2101c7385f5e2 /CakeSummary.R
parentbe6d42ef636e8e1c9fdcfa6f8738ee10e885d75b (diff)
Version 2.0v2.0
Diffstat (limited to 'CakeSummary.R')
-rw-r--r--CakeSummary.R174
1 files changed, 126 insertions, 48 deletions
diff --git a/CakeSummary.R b/CakeSummary.R
index ec42e75..6c7c661 100644
--- a/CakeSummary.R
+++ b/CakeSummary.R
@@ -1,10 +1,10 @@
-# $Id: CakeSummary.R 216 2011-07-05 14:35:03Z nelr $
+# $Id$
# 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
+# Parts Tessella (Rob Nelson/Richard Smith/Tamar Christina) for Syngenta: Copyright (C) 2011 Syngenta
+# Tessella Project Reference: 6245, 7247
# 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
@@ -19,12 +19,32 @@
# 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 extra decay times for the other k values in DFOP
+# and HS
+# Arguments
+# type - type of compartment
+# parms.all - list of parameters
+CakeExtraDT<-function(type, parms.all) {
+ DT50k1 = NA
+ DT50k2 = NA
+ if (!is.na(parms.all["k1"]) && !is.na(parms.all["k2"])) {
+ k1 = parms.all["k1"]
+ k2 = parms.all["k2"]
+ DT50k1 = log(2)/k1
+ DT50k2 = log(2)/k2
+ }
+ ret<-c(DT50k1, DT50k2)
+ names(ret)<-c("k1", "k2")
+ ret
+}
+
# 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) {
+# sannMaxIter - the maximum amount of iterations to do for SANN
+CakeDT<-function(type, obs_var, parms.all, sannMaxIter) {
if (type == "SFO") {
# Changed to new parameterisation
#k_names = grep(paste("k", obs_var, sep="_"), names(parms.all), value=TRUE)
@@ -47,11 +67,17 @@ CakeDT<-function(type, obs_var, parms.all) {
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)
+
+ DTminBounds <- 0.001
+ DTminInitial <- 1
+ DT50.temp <- optim(DTminInitial, f, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 50, lower = DTminBounds)
+ DT50.o = DT50.temp$par
+ DT50.converged = DT50.temp$convergence == 0
+ DT50 = ifelse(!DT50.converged, NA, DT50.o)
+ DT90.temp <- optim(DTminInitial, f, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 90, lower = DTminBounds)
+ DT90.o = DT90.temp$par
+ DT90.converged = DT90.temp$convergence == 0
+ DT90 = ifelse(!DT90.converged, NA, DT90.o)
}
if (type == "HS") {
k1 = parms.all["k1"]
@@ -84,12 +110,16 @@ CakeDT<-function(type, obs_var, parms.all) {
((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
+ DTmin <- 0.001
+ DT50.temp <- optim(DTmin, f_50, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 50, lower = DTmin)
+ DT50.o = DT50.temp$par
+ DT50.converged = DT50.temp$convergence == 0
+ if (!DT50.converged) 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
+ DT90.temp <- optim(DTmin, f_90, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 90, lower = DTmin)
+ DT90.o = DT90.temp$par
+ DT90.converged = DT90.temp$convergence == 0
+ if (!DT90.converged) DT90 = NA else DT90 = DT90.o
}
ret<-c(DT50, DT90)
names(ret)<-c("DT50","DT90")
@@ -104,32 +134,48 @@ CakeDT<-function(type, obs_var, parms.all) {
# 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) {
+CakeChi2<-function(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini) {
# 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
+
+ if(getRversion() >= '2.14.0'){
+ # You would usually call mkinfit to fit the data, however it seems that we've already fit the data
+ # bundle <- mkinfit(mkinmod, observed, parms.ini = parms.optim, state.ini = state.ini) #, err='err')
+ # Somewhere else. Calling mkinfit may result in an error during fitting. So Instead we just create
+ # the values that the new version of mkinerrmin expects. There is no class check so it's fine.
+ bundle<-vector()
+ bundle$predicted <- predicted_long
+ bundle$observed <- observed
+ bundle$obs_vars <- obs_vars
+ bundle$par <- c(parms.optim, state.ini.optim)
+ errmin.overall <- mkinerrmin(bundle)
+
+ errmin.overall
+ } else {
+ 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
}
- errmin
}
# Compute additional statistics (r²/efficiency) of obs v pred for each compartment and the whole problem
@@ -152,15 +198,20 @@ CakeAdditionalStats<-function(obsvspred){
}
# Summarise a fit - used for CakeIrlsFit and CakeOlsFit
-summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov = FALSE,...) {
+summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, halflives = 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
+ # covar <- try(ginv(0.5*object$hessian)) # unscaled covariance, calculate using generic inversion
rdf <- object$df.residual
resvar <- object$ssr / rdf
modVariance <- object$ssr / length(object$residuals)
+
+ # if(!is.numeric(covar)){
+ # covar <- try(ginv(0.5*object$hessian))
+ # }
if (!is.numeric(covar)) {
# message <- "Cannot estimate covariance; system is singular"
@@ -181,13 +232,33 @@ summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov
warning(message)
} else message <- "ok"
- if(cov && !is.null(covar)) {
+ if(cov && !is.null(covar)) {
+
+ # Filter the values for t-test, only apply t-test to k-values
+ t.names <- grep("k(\\d+)|k_(.*)", names(tval), value = TRUE)
+ t.rest <- setdiff(names(tval), t.names)
+ t.values <- c(tval)
+ t.values[t.rest] <- NA
+ t.result <- pt(t.values, rdf, lower.tail = FALSE)
+
+ # Now set the values we're not interested in for the lower
+ # and upper bound we're not interested in to NA
+ t.param = c(param)
+ t.param[t.names] <- NA
+ # calculate the 90% confidence interval
+ alpha <- 0.10
+ lci90 <- t.param + qt(alpha/2,rdf)*se
+ uci90 <- t.param + qt(1-alpha/2,rdf)*se
+
+ # calculate the 95% confidence interval
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)
+ lci95 <- t.param + qt(alpha/2,rdf)*se
+ uci95 <- t.param + qt(1-alpha/2,rdf)*se
+
+ param <- cbind(param, se, tval, t.result, lci90, uci90, lci95, uci95)
dimnames(param) <- list(pnames, c("Estimate", "Std. Error",
- "t value", "Pr(>t)", "Lower CI", "Upper CI"))
+ "t value", "Pr(>t)", "Lower CI (90%)", "Upper CI (90%)", "Lower CI (95%)", "Upper CI (95%)"))
+
ans <- list(residuals = object$residuals,
residualVariance = resvar,
sigma = sqrt(resvar),
@@ -219,14 +290,15 @@ summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov
ans$errmin <- object$errmin
ans$penalties <- object$penalties
if(distimes) ans$distimes <- object$distimes
+ if(halflives) ans$extraDT50 <- object$extraDT50
if(ff) ans$ff <- object$ff
if(length(object$SFORB) != 0) ans$SFORB <- object$SFORB
- class(ans) <- c("summary.CakeFit","summary.mkinfit", "summary.modFit")
+ class(ans) <- c("summary.CakeFit","summary.mkinfit", "summary.cakeModFit")
return(ans)
}
# Print a summary. Used for CakeIrlsFit, CakeOlsFit and CakeMcmcFit
-# Expanded from print.summary.modFit
+# Expanded from print.summary.cakeModFit
print.summary.CakeFit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
cat("\nEquations:\n")
print(noquote(as.character(x[["diffs"]])))
@@ -239,7 +311,7 @@ print.summary.CakeFit <- function(x, digits = max(3, getOption("digits") - 3), .
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, ...)
@@ -256,7 +328,13 @@ print.summary.CakeFit <- function(x, digits = max(3, getOption("digits") - 3), .
if(printdistimes){
cat("\nEstimated disappearance times:\n")
print(x$distimes, digits=digits,...)
- }
+ }
+
+ printextraDT50 <- !is.null(x$extraDT50)
+ if(printextraDT50){
+ cat("\nHalf Lives for Extra k values:\n")
+ print(x$extraDT50, digits=digits,...)
+ }
printff <- !is.null(x$ff)
if(printff){

Contact - Imprint