%\documentclass[10pt]{article} \documentclass[11pt,letterpaper]{article} \usepackage{Sweave} \usepackage{fullpage} \usepackage{graphicx} %If you want to include postscript graphics \DeclareGraphicsExtensions{.pdf} % for PNG conversion by htlatex \usepackage{wrapfig} \usepackage{verbatim} \usepackage[tableposition=top]{caption} \usepackage{ifthen} \usepackage[utf8]{inputenc} \usepackage{amsmath} \usepackage{amsfonts} \setlength{\parindent}{0pt} \setlength{\parskip}{1ex plus 0.5ex minus 0.2ex} \newcommand{\code}[1]{{\tt #1}} \begin{document} \title{Mineralization of terbuthylazine} \author{Anders Nielsen \& John Nash} \date{\today} \maketitle <>= ## run with cacheSweave if possible: ## require("cacheSweave"); Sweave("min.Rnw",driver=cacheSweaveDriver) if (require(cacheSweave)) { setCacheDir("min_cache") } @ \indent \section{The problem} Terbuthylazine is a herbicide used in agriculture. It is a so-called s-triazin like atrazine, which has been banned in Denmark after suspicion of causing cancer. Terbuthylazine can be bound to the soil, but free terbuthylazine can be washed into the drinking water. Some bacteria can mineralize it. This data is part of a larger experiment to determine the ability of certain bacteria to mineralize terbuthylazine, and to estimate the mineralization rate. \subsection{System} The experiment starts out with all the terbuthylazine being free $F_0=100\%$, $B_0=0\%$, and $M_0=0\%$. Over the time of the experiment the terbuthylazine will gradually be bound to the soil ($B$), gradually be mineralized ($M$), and some of what is bound to the soil is released to be free ($B$). The system will be described via a system of ordinary differential equations: \begin{align*} \frac{dB_t}{dt} &= -k_1 B_t + k_2 F_t, &B_0=0 \\ \frac{dF_t}{dt} &= k_1 B_t - (k_2+k_3) F_t, &F_0=100\\ \frac{dM_t}{dt} &= k_3 F_t, &M_0=0\\ \end{align*} \setkeys{Gin}{width=.8\textwidth} \begin{center} <>= par(mar=c(.1,.1,.1,.1)) plot(NA,NA,xlim=0:1, ylim=0:1, axes=FALSE, xlab='', ylab='') #abline(v=0:10/10, col=grey(.5)) #abline(h=0:10/10, col=grey(.5)) rect(0,.25,.2,.75, lwd=3) rect(.4,.25,.6,.75, lwd=3) rect(.8,.25,1,.75, lwd=3) text (x=c(.1,.5,.9), y = c(.5,.5,.5), labels = c('B','F','M'), adj = .5, cex=4) arrows(.2, .6, .4, .6, lwd=3) text (x=.3, y =.7, labels = expression(k[1]), adj = .5, cex=2) arrows(.4, .4, .2, .4, lwd=3) text (x=.3, y =.3, labels = expression(k[2]), adj = .5, cex=2) arrows(.6, .5, .8, .5, lwd=3) text (x=.7, y =.6, labels = expression(k[3]), adj = .5, cex=2) @ \end{center} The system is closed, so the amount mineralized at any given time is \begin{align*} M_t = 100 - B_t - F_t, \end{align*} so the system can be simplified by defining $X_t = (B_t, F_t)'$. The simplified system is: \begin{align*} \frac{dX_t}{dt} &= \underbrace{\begin{pmatrix}-k_1&k_2\\k_1&-(k_2+k_3)\end{pmatrix}}_{A}X_t, &X_0=\left({0 \atop 100}\right) \end{align*} This system is a linear ODE system, so it can be solved. Here we will express the solution via the matrix exponential function. \begin{align*} X_t= e^{At}X_0 \end{align*} The amount mineralized is measured 26 times throughout a year. <>= library(Hmisc) tab<-read.table("../DATA/min.dat", skip=3, header=FALSE) colnames(tab)<-c("Time", "Mineralized") latex(tab,file="", rowname=NULL,first.hline.double=FALSE, caption="Percentages mineralized at different times.", colnamesTexCmd="bf") @ A simple statistical model for these observations assumes independent normal measurement noise. \begin{align*} M_{t_i} \sim \mathcal{N}(100-\Sigma X_{t_i}, \sigma^2), \quad \mbox{independent, and with $X_{t_i}=e^{At_i}X_0$}. \end{align*} \section{Implementation in AD Model Builder} \subsection{Data} The data file prepared for reading into AD Model Builder is \VerbatimInput[frame=single,numbers=left,samepage=true,fontsize=\scriptsize]{../DATA/min.dat} \subsection{AD Model Builder Code} The implementation follows the typical AD Model Builder template; first data is read in, then model parameters are declared, and finally the negative log likelihood is coded. { \VerbatimInput[frame=single,numbers=left,samepage=true,fontsize=\scriptsize]{../ADMB/min.tpl} %\begin{scriptsize} %\linespread{.8} %\verbatiminput{../ADMB/min.tpl} %\end{scriptsize} } \subsection{Running} The model can be run from the command line by compiling and then executing the produced binary, but this can also be accomplished from within the R console like this: <>= system('admb min') system.time(system('./min')) @ <>= load('../ADMB/fit.RData') admb<-results @ as long as we ensure that the R working directory is set to the directory containing the AD Model Builder code for the problem \verb|min.tpl| and the corresponding data file \verb|min.dat|. The time to optimize the likelihood and calculate the Hessian was \Sexpr{admb$time["fit"]} seconds. \subsection{Results} After running the model we can plot the fit to make sure all went well. \setkeys{Gin}{width=0.7\textwidth} \begin{center} <>= dat<-read.table('../DATA/min.dat', skip=3, head=FALSE) plot(dat,xlab='Time (days)', ylab='M (%)', las=1) load('../ADMB/fit.RData') #min<-read.table('min.std', skip=1, head=FALSE) im<-which(results$fit$names=='M') min<-results$fit$est[im] sdmin<-results$fit$std[im] lines(dat[,1],min, col='red') lines(dat[,1],min+2*sdmin, lty='dashed', col='red') lines(dat[,1],min-2*sdmin, lty='dashed', col='red') @ \end{center} Estimates and standard deviations of model parameters and derived quantities are located in the file \verb|min.std|. From the figure it is evident that a more suitable model for this data set would include serial correlation, but for the purpose of comparing the three statistical tools with a simple we will ignore this for now, and treat it in a later section for just \verb|R| and \verb|ADMB|. % \begin{scriptsize} % \verbatiminput{../ADMB/min.std} % \end{scriptsize} \section{Implementation in R} Using \Sexpr{R.version$version.string} and package versions: <>= usedpkgs <- sort(c("optimx","Matrix","R2jags")) i1 <- installed.packages() print(i1[usedpkgs,"Version"],quote=FALSE) @ \subsection{Data} The same data file is used as for the AD Model Builder implementation, which is read in ignoring the first three lines in the file with the command: <>= dat<-read.table('min.dat', skip=3, header=FALSE) @ \subsection{R code} The first attempt used the "optim" function in R with its default settings. The optimization was initialized by the same values as the AD Model Builder program, and parametrized in exactly the same way. <>= dat<-read.table("../DATA/min.dat", skip=3, header=FALSE) library(Matrix) @ <>= nlogL<-function(theta){ k<-exp(theta[1:3]) sigma<-exp(theta[4]) A<-rbind( c(-k[1], k[2]), c( k[1], -(k[2]+k[3])) ) x0<-c(0,100) sol<-function(t)100-sum(expm(A*t)%*%x0) pred <- sapply(dat[,1],sol) -sum(dnorm(dat[,2],mean=pred,sd=sigma, log=TRUE)) } (s1 <- system.time(fit<-optim(c(-2,-2,-2,-2),nlogL,hessian=TRUE))) fit$value fit$convergence @ The default use of \code{optim} used \Sexpr{round(s1["elapsed"])} seconds, and reported successful completion. This is however not the correct solution. The AD Model Builder program reported a minimum negative log likelihood of 0.939, but R reported a minimum of \Sexpr{round(fit$value,3)}. Of the built-in optimizers in optim only "L-BFGS-B" managed to find the correct solution. <>= system.time(fit<-optim(c(-2,-2,-2,-2),nlogL,hessian=TRUE, method='L-BFGS-B')) fit$value fit$convergence @ The recent add-on package to R called "optimx" allows the user to try out many different build-in optimizers for a given objective function, and compare their solutions, and running time. (In this document, a developmental version of from R-forge is used.) We use the faster form of our function, and have also brought in the gradient expressions derived elsewhere (see the writeup \textbb{optderivs} in the TOOLS section of this project collection). <>= library(optimx) @ <>= fit <- optimx(c(-2,-2,-2,-2),nlogL,hessian=TRUE,control=list(all.methods=TRUE)) @ <>= fit @ <>= optimxres <- subset(data.frame(mnames=unlist(fit$method), values=unlist(fit$fvalues), times=unlist(fit$xtimes), KKT1=unlist(fit$KKT1), KKT2=unlist(fit$KKT2)), values-min(values)<1e-3) @ From this table we see that for this particular model the fastest optimizer in R which actually found the correct minimum was \Sexpr{with(optimxres,mnames[which.min(times)])} and the time was ca.~\Sexpr{round(min(optimxres$times))} seconds. Of the \Sexpr{length(fit$method)} built-in optimizers, \Sexpr{nrow(optimxres)} did find the correct solution. Based on the new convergence criteria (the columns \code{KKT1} and \code{KKT2} where both should be TRUE), only the \code{nlminb} model fit is unproblematic. \section{Implementation in BUGS (jags)} \subsection{Data} The same data file is used as for the AD Model Builder implementation. \subsection{BUGS code} The implementation in BUGS comes in two parts; (1) the BUGS model description (in a separate file), and (2) an R script reading in data and calling the model with data and specifying the starting points for the chains and desired number of iterations. First the model description file: \VerbatimInput[frame=single,numbers=left,samepage=true,fontsize=\scriptsize]{../BUGS/model.txt} % \begin{scriptsize} % \linespread{.8} % \verbatiminput{../BUGS/model.txt} % \end{scriptsize} The next section will show the second part. \subsection{Running} The R commands for reading in data and parsing it to the BUGS model. <>= library(R2jags) load.module('msm') # the module containing the matrix exponential @ <>= dat <- read.table('../DATA/min.dat', skip=3, header=FALSE) time <- dat[,1] M <- dat[,2] noObs <- length(time) set.seed(12345) jags.data <- list("noObs","time","M") jags.params <- c("logK1","logK2","logK3","sigma") jags.inits <- function(){ list("logK1"=-2+.2*rnorm(1),"logK2"=-2+.2*rnorm(1),"logK3"=-2+.2*rnorm(1), "sigma"=exp(-2+.2*rnorm(1)), .RNG.name="base::Wichmann-Hill", .RNG.seed=round(runif(1)*1000000)) } @ <>= time<-unname(system.time( jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params, n.iter=10000, n.thin=1,n.burnin=0,model.file="model.txt") )["elapsed"]) @ Running the 3 chains with 10000 iterations in each took approximately \Sexpr{round(time)} seconds. \subsection{Results} \setkeys{Gin}{width=0.9\textwidth} \begin{center} <>= require(R2jags) load("../BUGS/fit.RData") nf <- layout(rbind(1:4,5:8), heights=c(1,3)) par(mar=c(3,3,3,.5)) plot(density(as.vector(mcmcall[,c(2,7,12)])), main=expression(log(k[1]))) plot(density(as.vector(mcmcall[,c(2,7,12)+1])), main=expression(log(k[2]))) plot(density(as.vector(mcmcall[,c(2,7,12)+2])), main=expression(log(k[3]))) plot(density(as.vector(log(mcmcall[,c(2,7,12)+3]))), main=expression(log(sigma))) x<-mcmcall[,c(2,7,12)] matplot(x, 1:nrow(mcmcall), type='l', main=paste(round(mean(x),2),"(",round(sd(as.vector(x)),3),")", sep='')) x<-mcmcall[,c(2,7,12)+1] matplot(x, 1:nrow(mcmcall), type='l', main=paste(round(mean(x),2),"(",round(sd(as.vector(x)),3),")", sep='')) x<-mcmcall[,c(2,7,12)+2] matplot(x, 1:nrow(mcmcall), type='l', main=paste(round(mean(x),2),"(",round(sd(as.vector(x)),3),")", sep='')) x<-log(mcmcall[,c(2,7,12)+3]) matplot(x, 1:nrow(mcmcall), type='l', main=paste(round(mean(x),2),"(",round(sd(as.vector(x)),3),")", sep='')) @ \end{center} \end{document}