diff options
Diffstat (limited to 'CakeMcmcFit.R')
| -rw-r--r-- | CakeMcmcFit.R | 104 |
1 files changed, 88 insertions, 16 deletions
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") |