diff options
| -rw-r--r-- | CakeCost.R | 38 | ||||
| -rw-r--r-- | CakeInit.R | 5 | ||||
| -rw-r--r-- | CakeIrlsFit.R | 146 | ||||
| -rw-r--r-- | CakeIrlsPlot.R | 4 | ||||
| -rw-r--r-- | CakeMcmcFit.R | 104 | ||||
| -rw-r--r-- | CakeModel.R | 100 | ||||
| -rw-r--r-- | CakeNullPlot.R | 17 | ||||
| -rw-r--r-- | CakeOlsFit.R | 15 | ||||
| -rw-r--r-- | CakeOlsPlot.R | 72 | ||||
| -rw-r--r-- | CakePenalties.R | 2 | ||||
| -rw-r--r-- | CakeSummary.R | 174 |
11 files changed, 548 insertions, 129 deletions
@@ -23,7 +23,7 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see <http://www.gnu.org/licenses/># # -# $Id: CakeCost.R 216 2011-07-05 14:35:03Z nelr $ +# $Id$ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL, weight = "none", scaleVar = FALSE, cost = NULL, ...) { @@ -180,6 +180,24 @@ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL, Res <- (ModVar- obsdat) res <- Res/Err resScaled <- res*Scale + + # print("==========================="); + # print("#Names[i]"); + # print(Names[i]); + # print("xDat"); + # print(xDat); + # print("obsdat"); + # print(obsdat); + # print("ModVar"); + # print(ModVar); + # print("1/Err"); + # print(1/Err); + # print("Res"); + # print(Res); + # print("res"); + # print(res); + # print("==========================="); + Residual <- rbind(Residual, data.frame( name = Names[i], @@ -191,13 +209,17 @@ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL, res = res)) CostVar <- rbind(CostVar, - data.frame( - name = Names[i], - scale = Scale, - N = length(Res), - SSR.unweighted = sum(Res^2), - SSR.unscaled = sum(res^2), - SSR = sum(resScaled^2))) + data.frame( + name = Names[i], + scale = Scale, + N = length(Res), + SSR.unweighted = sum(Res^2), + SSR.unscaled = sum(res^2), + SSR = sum(resScaled^2))) + + # print("************************") + # print(CostVar) + # print("************************") } # end loop over all observed variables ## SSR @@ -1,4 +1,4 @@ -#$Id: CakeInit.R 221 2011-11-01 14:12:45Z nelr $ +#$Id$ # Purpose: Load the modules required by CAKE # Call with chdir=TRUE so it can load our source from the current directory @@ -20,6 +20,7 @@ # along with this program. If not, see <http://www.gnu.org/licenses/>. library(mkin) +library(Rsolnp) source("CakeModel.R") source("CakeOlsFit.R") source("CakeIrlsFit.R") @@ -36,4 +37,4 @@ options(error=recover) options(warn=-1) # When running from C#, this must match the C# CAKE version -CAKE.version<-"1.4" +CAKE.version<-"2.0" diff --git a/CakeIrlsFit.R b/CakeIrlsFit.R index bf20572..05eff0a 100644 --- a/CakeIrlsFit.R +++ b/CakeIrlsFit.R @@ -1,9 +1,9 @@ -#$Id: CakeIrlsFit.R 216 2011-07-05 14:35:03Z nelr $ +#$Id$ # # The CAKE R modules are based on mkin # Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011 Syngenta -# Authors: Rob Nelson, Richard Smith -# Tessella Project Reference: 6245 +# Authors: Rob Nelson, Richard Smith, Tamar Christina +# 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 @@ -22,7 +22,8 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ 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, atol=1e-6, control=list(),...) + scaleVar = FALSE, atol=1e-6, sannMaxIter = 10000, control=list(), + useSolnp = FALSE, method='L-BFGS-B',...) { ### This is a modification based on the "mkinfit" function. @@ -35,21 +36,25 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ observed <- subset(observed, name %in% names(mkinmod$map)) ERR <- rep(1,nrow(observed)) observed <- cbind(observed,err=ERR) + flag <- 0 + obs_vars = unique(as.character(observed$name)) if (is.null(names(parms.ini))) names(parms.ini) <- mkinmod$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 = mkinmod$diffs[[box]]))) } return(list(c(diffs))) } 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] @@ -57,6 +62,7 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ optim_initials <- setdiff(names(state.ini), fixed_initials) state.ini.optim <- state.ini[optim_initials] 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 = "_") @@ -67,8 +73,37 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ ############### Iteratively Reweighted Least Squares############# ## Start with no weighting - fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, - upper = upper,control=control,...) + + ## a prefitting step since this is usually the most effective method + if(useSolnp) + { + pnames=names(c(state.ini.optim, parms.optim)) + fn <- function(P){ + names(P) <- pnames + FF<<-costFunctions$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) + #optimmethod <- method0 + + flag <- 1 + 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') + ## now using submethod already + 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) + } + } + }else{ + fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, + upper = upper,control=control,...) + } + + ## print(fit$hessian) if(length(control)==0) { @@ -81,27 +116,109 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ control <- list(irls.control=irls.control) } } + irls.control <- control$irls.control maxIter <- irls.control$maxIter tol <- irls.control$tol #### + + if(length(mkinmod$map)==1 || useSolnp){ + ## there is only one parent just do one iteration: + maxIter <- 0 + + if(flag==1)## 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 + } + err1 <- sqrt(fit$var_ms_unweighted) + ERR <- err1[as.character(observed$name)] + observed$err <-ERR + } + niter <- 1 ## insure one IRLS iteration diffsigma <- 100 olderr <- rep(1,length(mod_vars)) while(diffsigma>tol & niter<=maxIter) { + 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 + } + + err <- sqrt(fit$var_ms_unweighted) ERR <- err[as.character(observed$name)] costFunctions$set.error(ERR) diffsigma <- sum((err-olderr)^2) cat("IRLS iteration at",niter, "; Diff in error variance ", diffsigma,"\n") olderr <- err - fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...) + + if(useSolnp) + { + flag <- 1 + a <- try(fit <- solnp(fit$par,fun=fn,LB=lower,UB=upper,control=control),silent=TRUE) + + if(class(a) == "try-error") + { + flag <- 0 + print('solnp fails during IRLS iteration, try PORT or other algorithm by users choice.This may takes a while. Do something else!') ## NOTE: because in kingui we switch off the warnings, we need to print out the message instead. + warning('solnp fails during IRLS iteration, switch to PORT or other algorithm by users choice') + + fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, method='Port',control=list()) + } + }else{ + fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...) + } + + ## print(fit$hessian) + niter <- niter+1 ### If not converged, reweight and fit } + + if(flag==1 && useSolnp){ + ## solnp used + optimmethod <- '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 + np <- length(c(state.ini.optim, parms.optim)) + fit$rank <- np + fit$df.residual <- length(fit$residuals) - fit$rank + } ########################################### fit$mkindiff <- mkindiff @@ -125,15 +242,19 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) - fit$errmin <- CakeChi2(observed, predicted_long, obs_vars, parms.optim, state.ini.optim) + fit$errmin <- CakeChi2(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini) parms.all = c(fit$par, parms.fixed) fit$penalties <- CakePenaltiesLong(parms.all, out_predicted, observed) 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) + fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter) } + + fit$extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.all) data <- merge(observed, predicted_long, by = c("time", "name")) @@ -142,6 +263,11 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ data$residual <- data$observed - data$predicted data$variable <- ordered(data$variable, levels = obs_vars) fit$data <- data[order(data$variable, data$time), ] + + fit$costFunctions <- costFunctions + fit$upper <- upper + fit$lower <- lower + class(fit) <- c("CakeFit", "mkinfit", "modFit") return(fit) } diff --git a/CakeIrlsPlot.R b/CakeIrlsPlot.R index 2de8910..5a7992f 100644 --- a/CakeIrlsPlot.R +++ b/CakeIrlsPlot.R @@ -1,4 +1,4 @@ -#$Id: CakeIrlsPlot.R 216 2011-07-05 14:35:03Z nelr $ +#$Id$ # Generates fitted curves so the GUI can plot them # Author: Rob Nelson (Tessella) # Tessella Project Reference: 6245 @@ -20,6 +20,6 @@ CakeIrlsPlot <- function(fit, xlim = range(fit$data$time), ...) { - CakeOlsPlot(fit, xlim, ...) + CakePlotInit(fit, xlim, ...) } diff --git a/CakeMcmcFit.R b/CakeMcmcFit.R index 209f904..f9a340e 100644 --- a/CakeMcmcFit.R +++ b/CakeMcmcFit.R @@ -1,9 +1,9 @@ -# $Id: CakeMcmcFit.R 216 2011-07-05 14:35:03Z nelr $ +# $Id$ # 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 -# Tessella Project Reference: 6245 +# Author: Rob Nelson, Tamar Christina +# 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 @@ -25,7 +25,9 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ 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,...) + ntrydr = 1, drscale = NULL, verbose = TRUE,fitstart=TRUE,seed=NULL,atol=1e-6, + sannMaxIter = 10000, control=list(), + useSolnp = FALSE, method='L-BFGS-B',...) { NAind <-which(is.na(observed$value)) @@ -54,6 +56,7 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ 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), @@ -72,7 +75,7 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ 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 @@ -85,11 +88,44 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ { ## ############# Get Initial Paramtervalues ############# ## Start with no weighting - fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, - upper = upper,...) - + + if(useSolnp) + { + 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' + 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 + } + } + else + { + fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.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) @@ -111,6 +147,21 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ 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 @@ -166,11 +217,13 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ 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) + fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter) } + fit$extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.all) # Fits to data odeini <- state.ini @@ -208,7 +261,7 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ predicted_long <- mkin_wide_to_long(out_transformed,time="time") obs_vars = unique(as.character(observed$name)) - fit$errmin <- CakeChi2(observed, predicted_long, obs_vars, parms.optim, state.ini.optim) + fit$errmin <- CakeChi2(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini) data<-observed data$err<-rep(NA,length(data$time)) @@ -230,7 +283,7 @@ CakeMcmcFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$ # Summarise a fit -summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov = FALSE,...) { +summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, halflives = TRUE, ff = TRUE, cov = FALSE,...) { param <- object$par pnames <- names(param) p <- length(param) @@ -255,12 +308,30 @@ summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, warning(message) } else message <- "ok" - # Parameter fits + # 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) - dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "t value", "Pr(>t)", "Lower CI", "Upper CI")) + 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 (90%)", "Upper CI (90%)", "Lower CI (95%)", "Upper CI (95%)")) # Residuals from mean of MCMC fit resvar <- object$ssr/ rdf @@ -285,6 +356,7 @@ summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, ans$fixed <- object$fixed ans$errmin <- object$errmin if(distimes) ans$distimes <- object$distimes + 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") diff --git a/CakeModel.R b/CakeModel.R index a93d7af..ce6bc01 100644 --- a/CakeModel.R +++ b/CakeModel.R @@ -1,5 +1,5 @@ # Was Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke -# $Id: CakeModel.R 216 2011-07-05 14:35:03Z nelr $ +# $Id$ # The CAKE R modules are based on mkin # Portions Johannes Ranke, 2010 @@ -22,12 +22,13 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see <http://www.gnu.org/licenses/>. -CakeModel <- function(...) +CakeModel <- function(..., use_of_ff = "max") { spec <- list(...) obs_vars <- names(spec) - # The returned model will be a list of character vectors, containing + if (!use_of_ff %in% c("min", "max")) + stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'") # differential equations, parameter names and a mapping from model variables # to observed variables. If possible, a matrix representation of the # differential equations is included @@ -69,21 +70,39 @@ CakeModel <- function(...) # Start a new differential equation for each new box new_diffs <- paste("d_", new_boxes, " =", sep="") - + + # Get the name of the box(es) we are working on for the decline term(s) + box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB + # Turn on sink if not specified otherwise if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE #@@@ADD SFO k HERE !!!!!!!!!!!!! # Construct and add SFO term and add SFO parameters if needed if(spec[[varname]]$type == "SFO") { + if (use_of_ff == "min") { # Minimum use of formation fractions + if(spec[[varname]]$sink) { # From p. 53 of the FOCUS kinetics report + k_term <- paste("k", new_boxes[[1]], "sink", sep="_") + nonlinear_term <- paste(k_term, new_boxes[[1]], sep=" * ") + spec[[varname]]$nlt<-nonlinear_term + new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) + new_parms <- k_term + ff <- vector() + decline_term <- paste(nonlinear_term, "*", new_boxes[[1]]) + } else { # otherwise no decline term needed here + decline_term = "0" + } + } else { k_term <- paste("k", new_boxes[[1]], sep="_") nonlinear_term <- paste(k_term, new_boxes[[1]], sep=" * ") spec[[varname]]$nlt<-nonlinear_term new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) new_parms <- k_term ff <- vector() - } + decline_term <- paste(nonlinear_term, "*", new_boxes[[1]]) + } + } # Construct and add FOMC term and add FOMC parameters if needed if(spec[[varname]]$type == "FOMC") { @@ -225,33 +244,64 @@ CakeModel <- function(...) } } } - model <- list(diffs = diffs, parms = parms, map = map) + model <- list(diffs = diffs, parms = parms, map = map, use_of_ff = use_of_ff) - # Create coefficient matrix if appropriate + # Create coefficient matrix if appropriate#{{{ if (mat) { boxes <- names(diffs) n <- length(boxes) m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes)) - for (from in boxes) { - for (to in boxes) { - if (from == to) { - #@@@@ OMIT NEXT LINE? !!!!! - k.candidate = paste("k", from, c(boxes, "sink"), sep="_") - k.candidate = sub("free.*bound", "free_bound", k.candidate) - k.candidate = sub("bound.*free", "bound_free", k.candidate) - k.effective = intersect(model$parms, k.candidate) - m[from,to] = ifelse(length(k.effective) > 0, - paste("-", k.effective, collapse = " "), "0") - } else { - k.candidate = paste("k", from, to, sep="_") - k.candidate = sub("free.*bound", "free_bound", k.candidate) - k.candidate = sub("bound.*free", "bound_free", k.candidate) - k.effective = intersect(model$parms, k.candidate) - m[to, from] = ifelse(length(k.effective) > 0, - k.effective, "0") + + if (use_of_ff == "min") { # Minimum use of formation fractions + for (from in boxes) { + for (to in boxes) { + if (from == to) { # diagonal elements + k.candidate = paste("k", from, c(boxes, "sink"), sep="_") + k.candidate = sub("free.*bound", "free_bound", k.candidate) + k.candidate = sub("bound.*free", "bound_free", k.candidate) + k.effective = intersect(model$parms, k.candidate) + m[from,to] = ifelse(length(k.effective) > 0, + paste("-", k.effective, collapse = " "), "0") + + } else { # off-diagonal elements + k.candidate = paste("k", from, to, sep="_") + if (sub("_free$", "", from) == sub("_bound$", "", to)) { + k.candidate = paste("k", sub("_free$", "_free_bound", from), sep="_") + } + if (sub("_bound$", "", from) == sub("_free$", "", to)) { + k.candidate = paste("k", sub("_bound$", "_bound_free", from), sep="_") + } + k.effective = intersect(model$parms, k.candidate) + m[to, from] = ifelse(length(k.effective) > 0, + k.effective, "0") + } + } + } + } else { # Use formation fractions where possible + for (from in boxes) { + for (to in boxes) { + if (from == to) { # diagonal elements + k.candidate = paste("k", from, sep="_") + m[from,to] = ifelse(k.candidate %in% model$parms, + paste("-", k.candidate), "0") + if(grepl("_free", from)) { # add transfer to bound compartment for SFORB + m[from,to] = paste(m[from,to], "-", paste("k", from, "bound", sep="_")) + } + if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB + m[from,to] = paste("- k", from, "free", sep="_") + } + m[from,to] = m[from,to] + } else { # off-diagonal elements + f.candidate = paste("f", from, "to", to, sep="_") + k.candidate = paste("k", from, to, sep="_") + # SFORB with maximum use of formation fractions not implemented, see above + m[to, from] = ifelse(f.candidate %in% model$parms, + paste(f.candidate, " * k_", from, sep=""), + ifelse(k.candidate %in% model$parms, k.candidate, "0")) + } } } - } + } model$coefmat <- m } diff --git a/CakeNullPlot.R b/CakeNullPlot.R index 3eacd67..f5484b4 100644 --- a/CakeNullPlot.R +++ b/CakeNullPlot.R @@ -1,4 +1,4 @@ -#$Id: CakeNullPlot.R 216 2011-07-05 14:35:03Z nelr $ +#$Id$ # Generates model curves so the GUI can plot them. Used when all params are fixed so there was no fit # The CAKE R modules are based on mkin # Based on code in IRLSkinfit @@ -19,7 +19,7 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see <http://www.gnu.org/licenses/>. -CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(observed$time), digits = max(3, getOption("digits") - 3), ...) +CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(observed$time), digits = max(3, getOption("digits") - 3), sannMaxIter = 10000, ...) { # Print the fixed values fixed <- data.frame(value = c(state.ini, parms.ini)) @@ -31,13 +31,20 @@ CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(obs obs_vars = unique(as.character(observed$name)) distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), row.names = obs_vars) + extraDT50<- data.frame(DT50 = rep(NA, 2), row.names = c("k1", "k2")) + for (obs_var in obs_vars) { type = names(model$map[[obs_var]])[1] if(exists("DT50")) distimes[obs_var, ] = c(DT50, DT90) - distimes[obs_var, ] = CakeDT(type,obs_var,parms.ini) + distimes[obs_var, ] = CakeDT(type,obs_var,parms.ini,sannMaxIter) } + cat("\nEstimated disappearance times:\n") print(distimes, digits=digits,...) + + extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.ini) + cat("\nHalf Lives for Extra k values:\n") + print(extraDT50, digits=digits,...) # Plot the model @@ -45,7 +52,7 @@ CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(obs odeini <- state.ini odenames <- model$parms - odeparms <- params.ini + odeparms <- parms.ini # Solve the system evalparse <- function(string) @@ -86,7 +93,7 @@ CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(obs predicted_long <- mkin_wide_to_long(out_transformed,time="time") # Calculate chi2 error levels according to FOCUS (2006) - errmin<-CakeChi2(observed, predicted_long, obs_vars, c(), c()) + errmin<-CakeChi2(model, observed, predicted_long, obs_vars, c("auto"), c(), state.ini, c("auto")) cat("\nChi2 error levels in percent:\n") errmin$err.min <- 100 * errmin$err.min print(errmin, digits=digits,...) diff --git a/CakeOlsFit.R b/CakeOlsFit.R index 75ac471..626a595 100644 --- a/CakeOlsFit.R +++ b/CakeOlsFit.R @@ -7,11 +7,11 @@ # from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn # inspired by summary.nls.lm -#$Id: CakeOlsFit.R 216 2011-07-05 14:35:03Z nelr $ +#$Id$ # This version has been modified to expect SFO parameterised as k and flow fractions # Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011 Syngenta -# Authors: Rob Nelson, Richard Smith -# Tessella Project Reference: 6245 +# Authors: Rob Nelson, Richard Smith, Tamar Christina +# 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 @@ -36,6 +36,8 @@ CakeOlsFit <- function(mkinmod, observed, plot = FALSE, quiet = FALSE, err = NULL, weight = "none", scaleVar = FALSE, atol = 1e-6, + sannMaxIter = 10000, + useSolnp = FALSE, method='L-BFGS-B', ...) { mod_vars <- names(mkinmod$diffs) @@ -294,20 +296,23 @@ CakeOlsFit <- function(mkinmod, observed, fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) # # Calculate chi2 error levels according to FOCUS (2006) - fit$errmin <- CakeChi2(observed, predicted_long, obs_vars, parms.optim, state.ini.optim) + fit$errmin <- CakeChi2(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini) # Calculate dissipation times DT50 and DT90 and formation fractions parms.all = c(fit$par, parms.fixed) 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")) + fit$ff <- vector() fit$SFORB <- vector() for (obs_var in obs_vars) { type = names(mkinmod$map[[obs_var]])[1] - fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all) + fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter) } + fit$extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.all) fit$penalties <- CakePenaltiesLong(parms.all, out_predicted, observed) # Collect observed, predicted and residuals diff --git a/CakeOlsPlot.R b/CakeOlsPlot.R index 199cc28..e92742c 100644 --- a/CakeOlsPlot.R +++ b/CakeOlsPlot.R @@ -1,4 +1,4 @@ -#$Id: CakeOlsPlot.R 216 2011-07-05 14:35:03Z nelr $ +#$Id$ # Generates fitted curves so the GUI can plot them # Based on code in IRLSkinfit # Author: Rob Nelson (Tessella) @@ -18,10 +18,70 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see <http://www.gnu.org/licenses/>. -CakeOlsPlot <- function(fit, xlim = range(fit$data$time), ...) +CakePlotInit <- function(fit, xlim = range(fit$data$time), ...) { - cat("<PLOT MODEL START>\n") + t.map.names <- names(fit$map) + metabolites <- grep("[A-Z]\\d", t.map.names, value = TRUE) + t.map.rest <- setdiff(t.map.names, metabolites) + + # Generate the normal graphs. + final <- CakeOlsPlot(fit, xlim) + final_scaled <- final + + if(length(metabolites) > 0){ + for(i in 1:length(metabolites)) + { + metabolite <- metabolites[i] + decay_var <- paste("k", metabolite, sep="_") + + # calculate the new ffm and generate the two ffm scale charts + regex <- paste("f_(.+)_to", metabolite, sep="_") + decays = grep(regex, names(fit$par), value = TRUE) + ffm_fitted <- sum(fit$par[decays]) + normal <- final + ffm_scale <- NULL + + # Generate the DT50=1000d and ffm as fitted. + k_new <- fit$par[decay_var]*fit$distimes[metabolite,]["DT50"]/1000; + fit$par[decay_var]<- k_new[metabolite,] + dt50_1000_ffm_fitted <- CakeOlsPlot(fit, xlim)[metabolite] + + naming <- c(names(final), paste(metabolite, "DT50_1000_FFM_FITTED", sep="_")) + normal <- c(final, dt50_1000_ffm_fitted) + names(normal) <- naming + final <- normal + + # Generate the scaled FFM + if(ffm_fitted != 0) + { + ffm_scale <- 1 / ffm_fitted + final_scaled <- final[c("time", metabolite, paste(metabolite, "DT50_1000_FFM_FITTED", sep="_"))] + final_scaled[t.map.rest] <- NULL; + final_frame <- as.data.frame(final_scaled) + new_names <- c(paste(metabolite, "DT50_FITTED_FFM_1", sep="_"), paste(metabolite, "DT50_1000_FFM_1", sep="_")) + names(final_frame) <- c("time", new_names) + final_frame[new_names]<-final_frame[new_names]*ffm_scale; + + cat("<PLOT MODEL START>\n") + + write.table(final_frame, quote=FALSE) + + cat("<PLOT MODEL END>\n") + } + } + } + + cat("<PLOT MODEL START>\n") + + write.table(final, quote=FALSE) + + cat("<PLOT MODEL END>\n") + + # View(final) +} +CakeOlsPlot <- function(fit, xlim = range(fit$data$time), scale_x = 1.0, ...) +{ solution = fit$solution if ( is.null(solution) ) { solution <- "deSolve" @@ -40,7 +100,7 @@ CakeOlsPlot <- function(fit, xlim = range(fit$data$time), ...) odeini <- parms.all[ininames] names(odeini) <- names(fit$diffs) - outtimes <- seq(0, xlim[2], length.out=101) + outtimes <- seq(0, xlim[2], length.out=101) * scale_x odenames <- c( rownames(subset(fit$start, type == "deparm")), @@ -111,7 +171,5 @@ CakeOlsPlot <- function(fit, xlim = range(fit$data$time), ...) out_transformed[var] <- rowSums(out[, fit$map[[var]]]) } } - print(out_transformed) - - cat("<PLOT MODEL END>\n") + return(out_transformed) } diff --git a/CakePenalties.R b/CakePenalties.R index 3b2a0df..1feffbe 100644 --- a/CakePenalties.R +++ b/CakePenalties.R @@ -1,4 +1,4 @@ -# $Id: CakePenalties.R 216 2011-07-05 14:35:03Z nelr $ +# $Id$ # The CAKE R modules are based on mkin # CAKE (6245), by Tessella, for Syngenta: Copyright (C) 2011 Syngenta # 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){ |