\documentclass[fleqn]{article} \usepackage[english]{babel} %% for texi2pdf ~ bug \usepackage{natbib} \usepackage[tableposition=top]{caption} \captionsetup {font=footnotesize,labelfont=bf,labelsep=period,singlelinecheck=false,width=9cm} \usepackage{color}\definecolor{gray}{rgb}{0.5,0.5,0.5} \usepackage{courier} \usepackage{float} \usepackage[height=8.5in]{geometry} \usepackage{graphicx} \usepackage{listings} \lstset{basicstyle=\footnotesize\ttfamily,lineskip=0.1pt, literate={~}{$\sim$}{2}} \usepackage{parskip} \usepackage{setspace}\setstretch{1.25} \usepackage[linkbordercolor={1 1 1},citebordercolor={1 1 1}, urlbordercolor={1 1 1},pdfpagemode=None]{hyperref} \newcommand{\current}{_\mathrm{current}} \newcommand{\init}{_\mathrm{init}} \newcommand{\mc}{\multicolumn} \newlength{\epswidth}\setlength{\epswidth}{30em} \newcommand{\code}[1]{{\tt #1}} \hyphenation{under-estimated} \bibliographystyle{chicago} \begin{document} \title{Orange tree growth: A demonstration and evaluation of nonlinear mixed-effects models in R, ADMB, and BUGS} \author{Arni Magnusson, Mark Maunder, and Ben Bolker} \date{\today} \maketitle \tableofcontents <>= require("knitr") opts_chunk$set(echo=FALSE,fig.align="center",fig.width=6,fig.height=4,out.width="0.7\\textwidth", tidy=FALSE) require("plyr") require("gplots") ## for rich.colors() require("ggplot2") theme_set(theme_bw()) require("directlabels") @ \section{Introduction} This document serves two purposes: \begin{itemize} \item Demonstrate how simple nonlinear mixed-effects models can be fitted in R, AD Model Builder, and BUGS \item Evaluate the estimation performance of models implemented in R and AD Model Builder \end{itemize} \section{Data} The data describe the growth of orange trees (Table \ref{tab:data}, Figure \ref{fig:dataplot}). The trunk circumference of 5 trees is measured at 7 different ages, giving a total of 35 datapoints. These data were used as example data by \citet[Ch.~8]{pinheiro_mixed-effects_2000} and have been discussed extensively elsewhere, e.g. \citet{madsen_introduction_2010}. The \code{Orange} data object is among the core datasets that come with R. \vspace{3ex} <>= Orange2 <- transform(Orange, Tree=factor(Tree,levels=1:5,labels=paste("Tree",1:5))) OTab <- rename(dcast(Orange2,age~Tree, value.var="circumference"), c(age="age (days)")) xtable(OTab,digits=0) @ \begin{table}[H] \centering\small\captionsetup{width=57ex} \caption{Growth of orange trees. Trunk circumference at breast height of 5 trees measured at 7 different ages.} \begin{tabular}{rrrrrr} \hline Age & \mc{5}{c}{Circumference (mm)} \\ (days) & Tree 1 & Tree 2 & Tree 3 & Tree 4 & Tree 5\\ \hline 118 & 30 & 33 & 30 & 32 & 30\\ 484 & 58 & 69 & 51 & 62 & 49\\ 664 & 87 & 111 & 75 & 112 & 81\\ 1004 & 115 & 156 & 108 & 167 & 125\\ 1231 & 120 & 172 & 115 & 179 & 142\\ 1372 & 142 & 203 & 139 & 209 & 174\\ 1582 & 145 & 203 & 140 & 214 & 177\\ \hline \end{tabular} \label{tab:data} \end{table} \vspace{2ex} <>= plotData <- function() { xoff <- 75 yoff <- 2 * c(1, 0, -1, 0, 0) matplot(unique(Orange$age), xtabs(circumference~age+Tree, Orange), type="o", xlim=c(0,1670), ylim=c(0,220), xlab="Age (days)", ylab="Circumference (mm)", pch=16, lty=1, lwd=1, col=rich.colors(5)) A <- max(Orange$age) text(rep(A,5)+xoff, Orange$circumference[Orange$age==A]+yoff, Orange$Tree[Orange$age==A]) } plotData() @ <>= p0 <- ggplot(Orange,aes(age,circumference,color=Tree))+ ylim(c(0,220))+labs(x="Age (days)",y="Circumference (mm)")+ geom_point()+ scale_colour_manual(values=rich.colors(5)) p1 <- p0 + geom_line() direct.label(p1,"last.qp") @ % \begin{figure}[H] % \centering\includegraphics[width=\epswidth]{data} % \caption{Growth of orange trees. Trunk circumference at breast height of 5 % trees measured at 7 different ages. Tree numbers are shown on the right.} % \label{fig:data} % \end{figure} \section{Model} \subsection{Description} \citet[pp.~356, 360]{pinheiro_mixed-effects_2000} % Pinheiro and Bates (2000, pp.~356,360) used the following mixed-effects logistic model to analyze the orange tree growth data, \begin{displaymath} y_{ij} \;=\; \frac{\phi_1+b_i}{\,1+\exp[-(x_{ij}\!-\!\phi_2)/\phi_3]\,} \,+\, \epsilon_{ij}, \qquad \epsilon_{ij}\sim N(0,\sigma^2) \end{displaymath} where $y_{ij}$ is the circumference of tree $i$ at time $j$, $x$ is the age in days, $\phi_1$ is the asymptote (maximum circumference), $\phi_2$ is the inflection point (age when trees have reached half maximum circumference), and $\phi_3$ is a scale parameter (days it takes to grow from 50\% to 73\% of maximum circumference\footnote{The 73\% comes from $1/[1\!+\!\exp\:\!(-1)]$.}). The asymptote varies between trees as a random effect: \begin{displaymath} b_i\sim N(0,\sigma_b^2) \end{displaymath} The traces in Figure~\ref{fig:dataplot} can be used to get sensible starting values for the three parameters and set the initial asymptote $\phi_1\!=\!200$, the inflection point $\phi_2\!=\!800$, and the scale parameter $\phi_3\!=\!400$. Symmetric confidence limits around the $\phi$ regression coefficients are constructed by multiplying the estimated standard error with the normal quantile $z$: \begin{displaymath} \mathrm{CI}_{\:\!\hat\phi_i} \;=\; \hat\phi_i \,\pm\, z_{\alpha/2}\,\widehat\mathrm{SE}_{\hat\phi_i} \end{displaymath} Asymmetric confidence limits around the $\sigma$ and $\sigma_b$ standard deviations are based on the standard error of the log-transformed parameters: \begin{displaymath} \mathrm{CI}_{\:\!\hat\sigma} \;=\; \exp\!\left(\log\hat\sigma \,\pm\, z_{\alpha/2}\,\widehat\mathrm{SE}_{\log\hat\sigma}\!\right) \end{displaymath} \subsection{R} Implementing the model in R is easy after loading the \code{nlme} package: <>= library("nlme") fm <- nlme(circumference~phi1/(1+exp(-(age-phi2)/phi3)), fixed=phi1+phi2+phi3~1, random=phi1~1|Tree, data=Orange, start=c(phi1=200,phi2=800,phi3=400)) @ <>= library("rbenchmark") f_ownfun <- function() nlme(circumference~phi1/(1+exp(-(age-phi2)/phi3)), fixed=phi1+phi2+phi3~1, random=phi1~1|Tree, data=Orange, start=c(phi1=200,phi2=800,phi3=400)) f_ssfun <- function() nlme(circumference~SSlogis(age,phi1,phi2,phi3), fixed=phi1+phi2+phi3~1, random=phi1~1|Tree, data=Orange, start=c(phi1=200,phi2=800,phi3=400)) ## ??? ## f_ssfun_s <- function() nlme(circumference~SSlogis(age,phi1,phi2,phi3), ## fixed=phi1+phi2+phi3~1, random=phi1~1|Tree, ## data=Orange) benchmark(f_ownfun(),f_ssfun(), ## f_ssfun_s(), replications=200, columns = c("test", "replications", "elapsed", "relative")) @ There is also a ``self-starting'' \code{SSlogis()} function in R, specifically for fitting logistic models, but the above is a basic general approach for any nonlinear mixed-effects model. (Using the \code{SSlogis()} function speeds up the fit by about 15\%, because in addition to providing initial conditions \code{SSlogis()} also returns an analytically computed gradient of the sum-of-squares function.) %\newpage \subsection{ADMB} In ADMB, the data and model are in two different text files, and the initial parameter values are in a third text file. The data are in a file called \code{ora.dat}, \lstset{comment=[l]{\#},commentstyle=\color{gray}} \lstinputlisting[frame=lines,aboveskip=1em,belowskip=2em]{../ADMB/ora.dat} the model code is in a file called \code{ora.tpl}, \renewcommand{\bfdefault}{b} \lstset{emph={DATA_SECTION,PARAMETER_SECTION,PROCEDURE_SECTION,GLOBALS_SECTION,REPORT_SECTION}} \lstset{emphstyle=\bfseries} \lstset{emph={[2]init_int,init_vector,init_matrix,init_number, random_effects_vector,sdreport_number,number,vector,matrix,number,int}} \lstset{emphstyle={[2]\color{magenta}}} \lstset{emph={[3]objective_function_value,define}} \lstset{emphstyle={[3]\color{red}}} \lstset{emph={[4]exp,for,sum,square,log,endl}} \lstset{emphstyle={[4]\color{blue}}} \lstset{comment=[l]{//},commentstyle=\color{gray}} \lstset{string=[b]{"},stringstyle=\color{green},showstringspaces=false} \lstinputlisting[frame=lines,aboveskip=1em,belowskip=2em]{../ADMB/ora.tpl} and the initial parameter values are in a file called \code{ora.pin}: \lstset{comment=[l]{\#},commentstyle=\color{gray}} \lstinputlisting[frame=lines,aboveskip=1em,belowskip=2em]{../ADMB/ora.pin} This particular ADMB implementation is simple, not taking advantage of efficiency improvements such as separable functions and estimating unscaled random effects \citep{skaug_automatic_2006,Fournieretal12}. For more advanced usage, see the \code{wildflowers}, \code{owl}, and \code{theta} examples at \url{https://groups.nceas.ucsb.edu/non-linear-modeling/projects}. %(Skaug and Fournier 2006, Fournier et al. 2012). The model is compiled with the shell command \begin{lstlisting}[basicstyle=\small\ttfamily] admb -r ora \end{lstlisting} and then run: \begin{lstlisting}[basicstyle=\small\ttfamily] ora \end{lstlisting} It can also be compiled and run from within R using the \code{R2admb} package: <>= library("R2admb") setup_admb() ## establish path etc. for ADMB @ <>= ## hack! shouldn't have to change directories but we do setwd("../ADMB") compile_admb("ora",re=TRUE) run_admb("ora") admb_res <- read_admb("ora") setwd("../WRITEUP") @ %% now the public code: <>= compile_admb("ora",re=TRUE) run_admb("ora") @ \subsection{BUGS} The BUGS model is also fairly simple in this case. \lstinputlisting[frame=lines,aboveskip=1em,belowskip=2em]{../BUGS/OrangeTreeRE_bugs.txt} It loops over trees \code{i} and observations \code{j}, computing the expected size (\code{eta[j,i]}) incorporating the random effect of the \code{i}\textsuperscript{th} tree (\code{phi[i]}), and specifying that \code{phi1[i]} is normally distributed. The rest of the code specifies the priors. The model can be run from within R using the \code{R2jags} package. As written the BUGS code requires the data to be arranged in \emph{wide} format as in Table~1 (and as the ADMB input is structured), rather than in \emph{long} format of the \code{Orange} variable (preferred by \code{nlme}: we use the \code{reshape} function to get a table whose first column is the age, with the other columns corresponding to measurements from specific trees. <>= OrangeTab <- reshape(Orange,timevar="Tree",direction="wide",idvar="age") @ <>= require("R2jags") BUGSData<-list(n = 7, K = 5, x = OrangeTab[,1], Y = OrangeTab[,-1]) set.seed(1001) ## set RNG seed for R inits<-list(phi1=200, phi2=800, phi3=400,tauC = 0.1, tau1=1, .RNG.name="base::Wichmann-Hill", ## set RNG seed/type for JAGS .RNG.seed=round(runif(1)*1000000)) @ <>= tfit_jags <- jags(model="../BUGS/OrangeTreeRE2_bugs.txt", data=BUGSData, parameters.to.save=c("phi1","phi2","phi3", "tau1","logSigmaB","logSigma","b1"), n.chains=1, inits=list(inits), progress.bar="none") @ \subsection{Simulations} $10\;000$ datasets are generated (Table \ref{tab:parameters}) and the R and ADMB model implementations are evaluated in terms of computational speed, convergence, bias, and coverage probability. \begin{table}[H] \centering\small\captionsetup{width=41ex} \caption{Parameter values used to simulate datasets for the second part of this study.} \begin{tabular}{rrrrrr} \hline Parameter & Value \\ \hline $\phi_1$ & 191.05\\ $\phi_2$ & 722.54\\ $\phi_3$ & 344.15\\ $\sigma$ & 7.85\\ $\sigma_b$ & 31.48\\ \hline \end{tabular} \label{tab:parameters} \end{table} \section{Results} \subsection{Model fit to original data} \subsubsection{R} The R command \code{summary(fm)} summarizes the model fit: <>= summary(fm) @ and \code{ranef(fm)} shows the random effects: <>= ranef(fm) @ \subsubsection{ADMB} The ADMB executable produces several output files. The negative log-likelihood and parameter estimates are in a file called \code{ora.par}, \lstset{comment=[l]{\#},commentstyle=\color{gray}} \lstinputlisting[frame=lines,escapechar=\%,aboveskip=1em,belowskip=2em]{../ADMB/ora.par} and standard errors and correlations are in a file called \code{ora.cor}: \lstinputlisting[frame=lines,escapechar=\%,aboveskip=1em,belowskip=2em]{../ADMB/ora.cor} The predicted values, which we reported in the \verb|REPORT_SECTION| of the TPL file, appear in a file called \code{ora.rep}: \lstinputlisting[frame=lines,escapechar=\%,aboveskip=1em,belowskip=2em]{../ADMB/ora.rep} Alternatively, the point estimates and standard errors (without correlations) can be found in a file called \code{ora.std}. The \code{R2admb} package can be useful for reading ADMB output back into R for plotting and analysis. %% fake, was read in earlier <>= admb_res <- read_admb("ora") @ creates an R object that contains all of the information from the ADMB fit. \subsubsection{BUGS} The default \code{print} method for \code{rjags} objects gives basic information: <>= tfit_jags @ We can use the \code{emdbook} and \code{coda} packages to transform and plot the results (trace and density plots): <>= library(emdbook) m <- as.mcmc.bugs(tfit_jags$BUGSoutput) library(coda) xyplot(m[,colnames(m)!="deviance"],layout=c(3,2),as.table=TRUE) densityplot(m[,colnames(m)!="deviance"],layout=c(3,2),as.table=TRUE) @ The \code{scapeMCMC} package gives alternative (prettier) plots. <>= library(scapeMCMC) plotTrace(m[,colnames(m)!="deviance"]) plotDens(m[,colnames(m)!="deviance"]) @ \subsubsection{Comparison} % In summary, the estimates are quite similar between R and ADMB (Table % \ref{tab:estimates}, Figure \ref{fig:fitplot}). % \vspace{1em} % \begin{table}[H] % \centering\small\captionsetup{width=72ex} % \caption{Estimated parameters, random effects, and negative log likelihood of % the model, as implemented in R and ADMB.} % \begin{tabular}{crr@{~}rrr@{~}r} % \hline % ~ & \mc{1}{l}{R} &&& ADMB\\ % ~ & \mc{1}{l}{Value} & \mc{2}{c}{95\% CI} % & \mc{1}{l}{Value} & \mc{2}{c}{95\% CI}\\ % \hline % $\phi_1$ & $191.05$ & $(159.41,$ & $222.69)$ % & $192.05$ & $(161.36,$ & $222.74)$\\ % $\phi_2$ & $722.54$ & $(653.69,$ & $791.38)$ % & $727.91$ & $(658.82,$ & $797.00)$\\ % $\phi_3$ & $344.15$ & $(290.98,$ & $397.32)$ % & $348.07$ & $(294.99,$ & $401.15)$\\ % $\sigma$ & $7.85$ & $ (6.09,$ & $10.11)$ % & $7.84$ & $ (6.09,$ & $10.10)$\\ % $\sigma_b$ & $31.48$ & $ (16.68,$ & $59.43)$ % & $31.65$ & $ (16.76,$ & $59.76)$\\ % \hline % $b_1$ & $-29.41$ &&& $-29.56$ \\ % $b_2$ & $31.56$ &&& $31.73$ \\ % $b_3$ & $-37.00$ &&& $-37.19$ \\ % $b_4$ & $40.02$ &&& $40.23$ \\ % $b_5$ & $-5.18$ &&& $-5.20$ \\ % \hline % $-\log L$ & $131.58$ &&& $131.57$ \\ % \hline % \end{tabular} % \label{tab:estimates} % \end{table} The ADMB and R estimates are extremely similar: the BUGS fits are slightly different, in part because BUGS is using the posterior mean rather than the posterior mode as a point estimate. <>= library("coefplot2") cc1 <- coeftab(admb_res)[1:5,] cc2 <- coeftab(tfit_jags)[rownames(cc1),] ## add Std. Error column: or could drop it from the others cc2 <- data.frame(cc2[,1],NA,cc2[,-1]) names(cc2) <- names(cc1) cc3 <- coeftab(fm,c("fixef","vcov")) cc3[4:5,"Estimate"] <- log(sqrt(cc3[5:4,"Estimate"])) rownames(cc3)[4:5] <- rownames(cc1)[4:5] cc3["logSigma",c(3,6)] <- log(intervals(fm)$sigma)[c("lower","upper")] cc3["logSigmaB",c(3,6)] <- log(intervals(fm)$reStruct$Tree)[c("lower","upper")] library(abind) library(reshape2) ccc <- dcast(melt(abind(ADMB=cc1,R=cc3,BUGS=cc2,along=3)), Var1+Var3~Var2) names(ccc) <- c("var","tool","lwr","lwr2","upr2","upr","est","stder") ccc <- transform(ccc, tool=factor(tool,levels=c("ADMB","R","BUGS")), var=factor(var,levels=c(paste0("phi",1:3), paste0("logSigma",c("","B"))))) ggplot(ccc,aes(x=tool,y=est,ymin=lwr,ymax=upr))+ geom_pointrange()+facet_wrap(~var,scale="free")+xlab("") @ <>= ## FIXME: coeftab needs updating for ADMB; hack for now? cc1R <- coeftab(admb_res)[6:10,] cc2R <- coeftab(tfit_jags)[1:5,] ## add Std. Error column: or could drop it from the others cc2 <- data.frame(cc2[,1],NA,cc2[,-1]) names(cc2) <- names(cc1) cc3R <- coeftab(fm,c("ranef")) cc3[4:5,"Estimate"] <- log(sqrt(cc3[5:4,"Estimate"])) rownames(cc3)[4:5] <- rownames(cc1)[4:5] library(abind) library(reshape2) ccc <- dcast(melt(abind(ADMB=cc1,BUGS=cc2,R=cc3,along=3)), Var1+Var3~Var2) ggplot(ccc,aes(x=Var3,y=Estimate,ymin=`2.5%`,ymax=`97.5%`))+ geom_pointrange()+facet_wrap(~Var1,scale="free")+xlab("") @ \vspace{1em} <>= pframe <- transform(Orange,circumference=predict(fm,level=1)) direct.label(p0 + geom_line(data=pframe),"last.qp") @ % \begin{figure}[H] % \centering\includegraphics[width=\epswidth]{fit} % \caption{Model fit to the data. The fitted values from R and ADMB are % indistinguishable in this figure. Tree numbers are shown on the right.} % \label{fig:fit} % \end{figure} \newpage \subsection{Performance with simulated data} \subsubsection*{Speed} Fitting $10\;000$ models took 7:37 minutes in R (0.05 sec/run) and 56:37 minutes in ADMB (0.34 sec/run) on an old laptop computer, so the model runs seven times faster in R than in ADMB. If a similar model was to be used in a computationally intensive simulation study, it would be worthwhile to take advantage of efficiency improvements such as separable functions and estimating unscaled random effects in ADMB \citep{skaug_automatic_2006,Fournieretal12}. % (Skaug and Fournier 2006, Fournier et al. 2012). \subsubsection*{Convergence} Out of $10\;000$ simulated datasets, 13 model runs did not converge, both in R and ADMB. In R, non-convergence was identified using the \verb|intervals| function, where 13 models returned either an error or an upper $\sigma_b$ confidence limit greater than $10^{10}$. Likewise, non-convergence in ADMB was identified from the standard error of $\log\hat\sigma_b$ being either \verb|NA| or greater than $10^{10}$. This occurred in the same set of 13 simulated datasets. To circumvent the problem of non-convergence, 13 additional simulated datasets were generated. The subsequent analysis of bias and coverage probability is therefore based on $10\;000$ converged simulations. \newpage \subsubsection*{Bias} In both R and ADMB, the $\phi$ parameter estimates are unbiased, but $\sigma$ and $\sigma_b$ are underestimated with a relative bias of around $-0.04$ and $-0.19$ (Table \ref{tab:bias}, Figure \ref{fig:bias}). In terms of bias, the difference between R and ADMB is negligible. \begin{table}[H] \centering\small\captionsetup{width=56ex} \caption{Median of $10\;000$ parameter estimates compared to the true parameter values used from the operating model. The relative bias is calculated as $\mathrm{median}((\hat\theta\!-\!\theta)/\theta)$.} \begin{tabular}{crrrrr} \hline ~ & & \mc{1}{l}{R} & & ADMB \\ ~ & True & Median & Bias & Median & Bias \\ \hline $\phi_1$ & $191.05$ & $190.78$ & $0.00$ & $191.59$ & $0.00$\\ $\phi_2$ & $722.54$ & $718.62$ & $-0.01$ & $723.31$ & $0.00$\\ $\phi_3$ & $344.15$ & $340.72$ & $-0.01$ & $344.04$ & $0.00$\\ $\sigma$ & $7.85$ & $ 7.51$ & $-0.04$ & $7.51$ & $-0.04$\\ $\sigma_b$ & $31.48$ & $ 25.52$ & $-0.19$ & $25.66$ & $-0.18$\\ \hline \end{tabular} \label{tab:bias} \end{table} \begin{figure}[H] \centering\includegraphics[width=0.9\textwidth]{sim-bias} \captionsetup{width=0.86\textwidth} \caption{Distribution of point estimates (Tukey boxplots) compared to the true parameter value (red line) of each parameter. Fourteen $\sigma$ estimates are outside the y-axis limits in the case of ADMB.} \label{fig:bias} \end{figure} \newpage \subsubsection*{Coverage probability} Analysis of $10\;000$ confidence intervals the 90\% confidence level shows that both R and ADMB generate confidence intervals that cover the true parameter less than 90\% of the time (Table \ref{tab:coverage}). The performance is poor for $\sigma_b$ and $\phi_1$, with coverage probability around 78\% and 82\%, but considerably better for the other parameters. \vspace{1em} \begin{table}[H] \centering\small\captionsetup{width=46ex} \caption{Coverage probability of 90\% confidence intervals generated using R and ADMB. Ideally, the coverage probability at this confidence level should be 0.900 for every parameter.} \begin{tabular}{crr} \hline ~ & R & ADMB\\ \hline $\phi_1$ & 0.824 & 0.813\\ $\phi_2$ & 0.884 & 0.879\\ $\phi_3$ & 0.886 & 0.877\\ $\sigma$ & 0.861 & 0.860\\ $\sigma_b$ & 0.783 & 0.787\\ \hline \end{tabular} \label{tab:coverage} \end{table} Analysis of confidence levels ranging from 0 to 99\% shows the same trends (Figure \ref{fig:coverage}). Both R and ADMB generate confidence intervals that are too narrow, especially for $\sigma_b$ and $\phi_1$. Overall, R and ADMB show similar performance in terms of coverage probability: R performs slightly better for the $\phi$ parameters, but ADMB performs slightly better for the problematic $\sigma_b$ parameter. \vspace{1em} \begin{figure}[H] \centering\includegraphics[width=\textwidth]{sim-ci} \captionsetup{width=0.9\textwidth} \caption{Coverage probability for confidence intervals evaluated at several confidence levels (0, 10, 20, 30, 40, 50, 60, 70, 80, 90, and 99\%). Each panel shows the performance of R (circle) and ADMB (cross) for a given parameter.} \label{fig:coverage} \end{figure} \newpage \section{Discussion} We have demonstrated how simple nonlinear mixed-effects models can be fitted in R, AD Model Builder, and BUGS. For this basic example, R and ADMB show similar estimation performance on the whole. To fit a mixed-effects logistic growth model to the orange tree data, it is easier and faster to use the \code{nlme} package in R, yielding similar results as ADMB. One possible reason to use ADMB for analyzing this dataset might be to explore other modelling options (statistical assumptions and methods) that are not provided by the \code{nlme} function in R. %\section{References} \small\sloppy\frenchspacing\setlength{\hyphenpenalty}{1000}\setstretch{1} \bibliography{../../nonlin} \end{document}