diff options
Diffstat (limited to 'CakeSummary.R')
| -rw-r--r-- | CakeSummary.R | 301 |
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) +} |