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