# Some of the CAKE R modules are based on mkin, # Based on mcmckinfit as modified by Bayer # Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta # Tessella Project Reference: 6245, 7247, 8361, 7414 # The CAKE R modules are 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 . # Performs an Markov chain Monte Carlo least squares fit on a given CAKE model. # # cake.model: The model to perform the fit on (as generated by CakeModel.R). # observed: Observation data to fit to. # parms.ini: Initial values for the parameters being fitted. # state.ini: Initial state (i.e. initial values for concentration, the dependent variable being modelled). # lower: Lower bounds to apply to parameters. # upper: Upper bound to apply to parameters. # fixed_parms: A vector of names of parameters that are fixed to their initial values. # fixed_initials: A vector of compartments with fixed initial concentrations. # quiet: Whether the internal cost functions should execute more quietly than normal (less output). # niter: The number of MCMC iterations to apply. # atol: The tolerance to apply to the ODE solver. # sannMaxIter: The maximum number of iterations to apply to SANN processes. # control: ... # useExtraSolver: Whether to use the extra solver for this fit (only used for the initial first fit). CakeMcmcFit <- function (cake.model, observed, parms.ini, state.ini, lower, upper, fixed_parms = NULL, fixed_initials, quiet = FALSE, niter = 1000, verbose = TRUE, seed=NULL, atol=1e-6, sannMaxIter = 10000, control=list(), useExtraSolver = FALSE, ...) { NAind <-which(is.na(observed$value)) mod_vars <- names(cake.model$diffs) observed <- subset(observed, name %in% names(cake.model$map)) ERR <- rep(1,nrow(observed)) observed <- cbind(observed,err=ERR) obs_vars = unique(as.character(observed$name)) if (is.null(names(parms.ini))) { names(parms.ini) <- cake.model$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 = cake.model$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] state.ini.fixed <- state.ini[fixed_initials] 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 = "_") } costFunctions <- CakeInternalCostFunctions(cake.model, state.ini.optim, state.ini.optim.boxnames, state.ini.fixed, parms.fixed, observed, mkindiff, quiet, atol=atol) bestIteration <- -1; costWithStatus <- function(P, ...){ r <- costFunctions$cost(P) if(r$cost == costFunctions$get.best.cost()) { bestIteration<<-costFunctions$get.calls(); cat(' MCMC best so far: c', r$cost, 'it:', bestIteration, '\n') } 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 seed = runif(1,0,2^31-1) } seed = as.integer(seed) set.seed(seed) ## ############# Get Initial Paramtervalues ############# ## Start with no weighting if(useExtraSolver) { a <- try(fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='Port',control=control)) if(class(a) == "try-error") { fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='L-BFGS-B',control=control) } } else { # modFit parameter transformations can explode if you put in parameters that are equal to a bound, so we move them away by a tiny amount. all.optim <- ShiftAwayFromBoundaries(c(state.ini.optim, parms.optim), lower, upper) fit <- modFit(costFunctions$cost, all.optim, lower = lower, upper = upper,...) } ## ############## One reweighted estimation ############################ ## Estimate the error variance(sd) tmpres <- fit$residuals oldERR <- observed$err err <- rep(NA,length(cake.model$map)) for(i in 1:length(cake.model$map)) { box <- names(cake.model$map)[i] ind <- which(names(tmpres)==box) tmp <- tmpres[ind] err[i] <- sd(tmp) } names(err) <- names(cake.model$map) ERR <- err[as.character(observed$name)] observed$err <-ERR costFunctions$set.error(ERR) olderr <- rep(1,length(cake.model$map)) diffsigma <- sum((err-olderr)^2) ## At least do one iteration step to get a weighted LS fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, ...) ## Use this as the Input for MCMC algorithm ## ########################## fs <- summary(fit) cov0 <- if(all(is.na(fs$cov.scaled))) NULL else fs$cov.scaled*2.4^2/length(fit$par) var0 <- fit$var_ms_unweighted costFunctions$set.calls(0); costFunctions$reset.best.cost() res <- modMCMC(costWithStatus,fit$par,...,jump=cov0,lower=lower,upper=upper,prior=NULL,var0=var0,wvar0=0.1,niter=niter,outputlength=niter,burninlength=0,updatecov=niter,ntrydr=1,drscale=NULL,verbose=verbose) # Construct the fit object to return fit$start <- data.frame(initial = c(state.ini.optim, parms.optim)) fit$start$type = c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim))) fit$start$lower <- lower fit$start$upper <- upper fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed)) fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) fit$mkindiff <- mkindiff fit$map <- cake.model$map fit$diffs <- cake.model$diffs # Replace mean from modFit with mean from modMCMC fnm <- function(x) mean(res$pars[,x]) fit$par <- sapply(dimnames(res$pars)[[2]],fnm) fit$bestpar <- res$bestpar fit$costfn <- costWithStatus # Disappearence times parms.all <- c(fit$par, parms.fixed) 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(k1 = rep(NA, length(names(cake.model$map))), k2 = rep(NA, length(names(cake.model$map))), row.names = names(cake.model$map)) for (compartment.name in names(cake.model$map)) { type = names(cake.model$map[[compartment.name]])[1] fit$distimes[compartment.name, ] = CakeDT(type,compartment.name,parms.all,sannMaxIter) fit$extraDT50[compartment.name, ] = CakeExtraDT(type, compartment.name, parms.all) } fit$ioreRepDT = CakeIORERepresentativeDT("Parent", parms.all) fit$fomcRepDT = CakeFOMCBackCalculatedDT(parms.all) # Fits to data odeini <- state.ini odeparms <- parms.all # If this step modelled the start amount, include it. Lookup in vectors returns NA not null if missing if(!is.na(odeparms["Parent_0"])) { odeini['Parent'] = odeparms["Parent_0"] } # Solve the system evalparse <- function(string) { eval(parse(text=string), as.list(c(odeparms, odeini))) } # Ensure initial state is at time 0 obstimes <- unique(c(0,observed$time)) odeini <- AdjustOdeInitialValues(odeini, cake.model, odeparms) out <- ode(y = odeini, times = obstimes, func = mkindiff, parms = odeparms, atol = atol) out_transformed <- PostProcessOdeOutput(out, cake.model, atol) fit$predicted <- out_transformed fit$penalties <- CakePenaltiesLong(odeparms, out_transformed, observed) predicted_long <- wide_to_long(out_transformed,time="time") obs_vars = unique(as.character(observed$name)) fit$errmin <- CakeChi2(cake.model, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini, fit$fixed) data<-observed data$err<-rep(NA,length(data$time)) data<-merge(data, predicted_long, by=c("time","name")) names(data)<-c("time", "variable", "observed","err-var", "predicted") data$residual<-data$observed-data$predicted data$variable <- ordered(data$variable, levels = obs_vars) fit$data <- data[order(data$variable, data$time), ] fit$atol <- atol sq <- data$residual * data$residual fit$ssr <- sum(sq) fit$seed = seed fit$res <- res class(fit) <- c("CakeMcmcFit", "mkinfit", "modFit") return(fit) } # Summarise a fit summary.CakeMcmcFit <- 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 mcmc <- object$res covar <- cov(mcmc$pars) rdf <- object$df.residual message <- "ok" rownames(covar) <- colnames(covar) <-pnames #se <- sqrt(diag(covar) * resvar) fnse <- function(x) sd(mcmc$pars[,x]) #/sqrt(length(mcmc$pars[,x])) se <- sapply(dimnames(mcmc$pars)[[2]],fnse) 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" # 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 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 modVariance <- object$ssr / length(object$data$residual) ans <- list(residuals = object$data$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) ans$diffs <- object$diffs ans$data <- object$data ans$additionalstats = CakeAdditionalStats(object$data) ans$seed <- object$seed ans$start <- object$start ans$fixed <- object$fixed ans$errmin <- object$errmin if(distimes){ ans$distimes <- object$distimes ans$extraDT50 <- object$extraDT50 ans$ioreRepDT <- object$ioreRepDT ans$fomcRepDT <- object$fomcRepDT } if(halflives) ans$halflives <- object$halflives if(ff) ans$ff <- object$ff class(ans) <- c("summary.CakeFit","summary.mkinfit", "summary.modFit") return(ans) }