project description

Summary
Increasingly, non-linear and complex models are applied as a tool for improving understanding of ecological systems. These statistical models are often used to test hypotheses and make inferences about ecological theories and management decisions based on available data. This explosion in the application of such models is due to rapid and current development of methodology to carryout statistical inference of complex nonlinear models and improvements in computer power (faster and multiple processors). While there are many tools available for statistical inference that differ in their effectiveness for specific applications, no formal comparisons have been conducted between various software packages. It is therefore important to identify which tools are most appropriate for given applications and to demonstrate how such tools can be used most effectively. We evaluate three open source software packages commonly used to carry out statistical inference of complex nonlinear models: OpenBUGS, AD Model Builder, and R. To test the strengths and weaknesses of each package, we will bring together experts in all three software packages and apply a common set of ecological models. Working directly with NCEAS informatics staff, we will produce a web-based guide regarding the utility of each package for particular applications that includes annotated model code for each package, the data sets used in the applications, and peer-reviewed articles. We will also identify how the different packages can be modified to improve their applicability to an array of complex nonlinear models that are essential for advancing ecological research. As statistical models are becoming increasingly more complex and ecologists are faced with a myriad of software options, the results of this project will provide support for ecologists and analysts across a broad spectrum of specialties.

Problem Statement
In the field of ecology, complex nonlinear models are essential tools for testing hypotheses about ecological theories and in making management decisions regarding populations and ecosystems. The increasing use of these models necessitates the development of effective tools for statistical inference. In recent years, with the development of methodology to carryout statistical inference of complex nonlinear models (e.g., Markov-chain Monte Carlo, McMC) and the improvements in computer power (faster and multiple processors), ecology has seen an explosion in the application of nonlinear models (e.g., hierarchical and integrated models) (Maunder et al. 2009). However, statistical inference using nonlinear models can be computationally demanding and the tools available for such calculations can have a dramatic effect on the analyses attempted. Reduced computation time facilitates the investigation, testing, and comparison of models and what can be done guides the development of theory. Historically, nonlinear models were limited to those with convenient algebraic solutions, often through approximation, and the pool of developers was limited by a formidable wall of mathematics (Efron and Tibshirani, 1993). The development of computer-based methods of statistical inference (e.g., McMC, iterative optimization) allowed quantitative scientist to compare their data with the model that they actually want to use rather than a model which has some mathematically convenient form (Clifford 1993, Bolker et al. 2009, Kéry 2010). For example, prior to the widespread introduction of simulation-based methods (e.g., McMC), Bayesian ideas could only be implemented in circumstances in which solutions could be obtained in closed-form using conjugate priors (Lunn et al. 2009). Now the desire for more appropriate and more intuitive priors outweighs the benefit of any such conjugacy (Lunn et al. 2009).

Ecological theory suggests non-linear models are appropriate to describe many systems. Even when linear approximations are adequate, nonlinear models are used to allow clear interpretation of parameters (Seber and Wild 1989). The term ‘non-linear models’ is often used to refer to models that are nonlinear in their parameters and cannot be made linear by transformations of the variables (Seber and Wild 1989). Here we do not hold ourselves to this definition, but instead generalize the terminology to refer to models that do not have a simple closed-form solution for statistical inference in either a Bayesian or frequentist framework and therefore require a more sophisticated, computationally-intensive approach to parameter estimation. These types of models include, for example, models with a nonlinear description of the ecological system, mixed-effects models (e.g., Bolker et al. 2009), and zero-inflated models (e.g., Minami et al. 2007). Many of these models have hierarchical structures (i.e. random effects and state-space models), are autoregessive (e.g., population dynamics models), and include multiple data types (Maunder 2004).

Tools available for statistical inference differ in their ease of use, speed, accuracy and general effectiveness for different applications. It is therefore important to identify which tools are most appropriate for a given application and to demonstrate how the tool can be used most effectively. A student entering the field of quantitative ecology today faces a daunting multiplicity of software environments and methodologies and some guidance of what to use for a particular problem is vital. A large number of software environments have been applied to the quantitative analysis of ecological systems. These include traditional programming languages (e.g., Fortran), common spreadsheet programs (e.g., Excel), and fourth-generation languages that address specific problems for quantitative analysis of ecological systems (e.g., AD Model Builder, R, OpenBUGS) (Maunder et al. 2009). An appropriate choice of software for investigating ecological theory depends on the purpose. For example, some software may be highly efficient for instructional or educational purposes, but inadequate for providing reasonable scientific advice for managers. Different modeling environments vary in their computational efficiency, ease of use, and quality of available documentation. The modern quantitative ecologist needs to balance time between learning new software systems and keeping abreast of growing scientific issues. Therefore, a guide to which software is appropriate for what type of applications and a set of examples for each software package is a huge advantage in deciding on whether a researcher should spend their time learning a new language or deciding between different languages to learn.

We limit the scope of the research to three specific packages because they are free, open source, and commonly used for computationally demanding ecological analysis. Open source is important because it allows for user input into development and improvement of the software. OpenBUGS (http://www.openbugs.info/w/) is a general software package for implementing Bayesian inference, AD Model Builder (ADMB; admb-project.org) is a general software package designed for statistical inference of nonlinear models in either Bayesian or Frequentist frameworks, and R (http://www.r-project.org/) is a general statistical software package that has several tools for statistical inference of nonlinear models.

The three software packages differ in several aspects and each modeling environment has its own advantages and disadvantages. R is a more general package, while OpenBUGS and ADMB are more specialized. ADMB cannot be used to implement models with discrete random effects, while OpenBUGS often struggles on autocorrelated time-series models (Lunn et al. 2009) and R users often encounter memory issues (Fox 2009). The optim() function in R works fairly well particularly if provided the analytical gradients, but these may not always be available and are automatically provided in ADMB. ADMB tends to have a steeper learning curve than the other packages, but does not require knowledge of Bayesian statistics unlike OpenBUGS. Although OpenBUGS provides a wider range of McMC schemes than ADMB, including reversible jump McMC for multiple model inference. Due to ADMB’s limited user and developer bases, the availability of useful functions is limited compared to R. OpenBUGS and particularly ADMB have limited capacity to display and organize results. Therefore, there are several tools for connecting these packages to R and other software packages. For example, BRugs and R2WinBUGS? are both interfaces that allow BUGS to be used interactively from within R.

There have been few comparisons of performance among the different packages and little advice on which package to use for what models. The work comparing ADMB and other packages (Schnute et al. 1998) is long out of date, very limited in scope and needs to be updated. Other studies focus more on the models used than the software that is applied (e.g. NCEAS research - McCarthy? et al. 2004).

Rationale for holding the working group at NCEAS
There are lots of users of these packages, but few developers. The developers are scattered around the world at different government research organizations and universities and in different fields of research. The only way to make a fair comparison among packages is to bring together experts from each software package in a neutral arena like an NCEAS workshop setting and implement a common testing framework. The experience will be invaluable, not only to determine which package to use for what applications, but also for participants to learn from each other and improve the performance of each package. It is important to do this in the context of applications because the improvement of these packages needs interaction with real problems. As analytical methods converge among fields of research, ecology software that crosses these disciplines will become more successful.

Proposed Activities
The goals of the proposal are to evaluate which software packages are appropriate for a type of application, illustrate how they are applied, and identify what improvements could be made. The research will help inform ecologist of the best available software for making inference with respect to a suite of problem driven applications. The results will be published in a series or primary publications and/or reports and a web based guide. We approach this in a series of steps.

Description of software packages
Each software package is first described in general and then how it can be applied to statistical inference using nonlinear models. This is important since most researchers will not be familiar with all three packages. The description will also allow an initial evaluation of the benefits of each package and what types of applications they have been applied to.

Identifications of applications
Each participant will be requested to provide several applications. The working group will evaluate all the applications and from these the working group will select a set of applications that covers the range of applications in ecological research. Care will be taken so that the applications that are selected do not favor one particular software package. The applications will also be chosen to ensure the success of the project and to maximize the benefits to the ecological community. Therefore, preference will be given to applications that have already been implemented in one of the software packages, for which the data can be made publically available, and address real management questions or ecological issues. These applications may include, but not be limited to, population dynamics, mark-recapture, occupancy, zero-inflated count models, and spatial models. See Table 1 for example applications.

Develop simulated data
Both real data and simulated data will be used for each application. The simulated data is necessary to evaluate the accuracy of the software packages. Simulated data will be created for each application based on the real data (e.g. a model will be applied to the real data and the parameter estimates will be used in a model to create the simulated data).

Development of models
For each application, a model will be developed in each software package. The models will be developed by the working group participants that are experts in using that package. This is to minimize differences in results due to the skill of the model developer.

Evaluation
The performance of the software packages will be evaluated for each application. Evaluation will be based on the following criteria:
• Ease of use (this will be subjective)
• Accuracy of both point estimates and estimates of uncertainty (based on simulated data)
• Computational speed and required computer resources (e.g. physical memory) of the applications


Workshops
The work will be conducted in two 5 Day meetings held at NCEAS. The first meeting will be used to summarize the software packages, select the applications to use, and develop the simulators. The participants will then develop and run the models while at their own institutions. The second meeting will provide a forum to present the results. The results will be summarized and manuscripts prepared.

Participants
Ben Bolker (PI), McMaster? University
Liza Comita, NCEAS
Beth Gardner (PI), Patuxent Wildlife Research Center, USGS/ North Carolina State University
Olivier Gimenez , Centre d’Écologie Fonctionnelle et Évolutive
Marc Kéry, Swiss Ornithological Institute
Eun Jung Kim, University of Hawaii
Cleridy Lennert-Cody, Inter-American? Tropical Tuna Commission
Arni Magnusson, Marine Research Institute
Mark Maunder (PI), Inter-American? Tropical Tuna Commission
Mihoko Minami , Keio University
John Nash, University of Ottawa
Anders Nielsen, Danish Technical University
Jim Regetz, NCEAS
Hans Skaug, University of Bergen
Elise Zipkin, Patuxent Wildlife Research Center

Casper W. Berg, Technical University of Denmark



Timetable of Activities
The first meeting will be held at the end of 2010 (e.g., November)
The second meeting will be held in early to mid 2011 (e.g., March)

Anticipated Results and Benefits
A report that compares and contrasts each software package, describes the results of the simulation tests, and describes the future needs for development (July 2011).
One or more manuscripts for publication based on the report (July 2011).
Example code and data for each application (July 2011)
Simulated data sets for each application (July 2011)
A web-based guide regarding the utility of each package for particular applications that includes annotated model code for each package, the data sets used in the applications, and peer-reviewed articles (July 2011).

Description of Software Packages
OpenBUGS
OpenBUGS is a software package that makes Bayesian sampling techniques available to the average analyst. The user specifies a statistical model by simply stating the relationships, which are often probabilistic, between related variables. The BUGS language contains a predefined vocabulary of distributions making coding easy, as long as the required distributions are supported. Technically, a BUGS model must correspond to a directed acyclic graph (DAG), in which random variables successively influence each other to produce the observed data, and Bayesian inference is performed using Gibbs Sampling or a variety of other McMC techniques. Fairly complex models can be developed with only a small amount of code that follows syntax akin to R (see below).

OpenBUGS, along with its predecessor WinBUGS, has been instrumental in raising awareness of Bayesian modelling among both academic and commercial communities internationally, and has enjoyed considerable success over its 20-year life span (Lunn et al. 2009). Currently, the WinBUGS version has over 30,000 registered users (Lunn et al. 2009). BUGS has greatly facilitated the routine analysis of many complex statistical problems including disease mapping , pharmacometrics, ecology, health-economics, genetics, archaeology, psychometrics, coastal engineering, educational performance, behavioural studies, econometrics, automated music transcription, sports modelling, fisheries stock assessment, and actuarial science (Lunn et al. 2009). The software is also used widely for the teaching of Bayesian modelling ideas to students and researchers the world over, and several texts use BUGS extensively for illustrating the Bayesian approach (e.g. King et al. 2009; Kéry 2010). It seems that the most fundamental reasons behind BUGS’ success is its simplicity, flexibility, and generality (Lunn et al. 2009).

R
R is a high-level language for manipulating, analyzing, and displaying data. R provides access to a huge number of algorithms useful for building a variety of statistical models (Ihaka and Gentleman 1996). R has attracted a large and growing user community and developer base. This is much due to the package system in R that permits individuals to participate in the development of R without the direct intervention of the R core group or without changes to the basic R system (Fox 2009). There are nearly 2000 contributed packages to the R project that implement a wide variety of algorithms from diverse fields (Fox 2009). These packages can be written in R script, or in another language (e.g. C++ to improve efficiency) and integrated into R. R also offers beginners a less intimidating environment than third-generation programming languages, like Fortran and C/C++. Numerous books have been published that use R for teaching statistical analysis (e.g. Bolker 2008). As a consequence, quite a few individuals and organizations now have a substantial investment in R. (Fox 2009). For example, the recently-initiated Fisheries Library in R (FLR) provides a collection of R tools that facilitate the construction of models representing fisheries and ecological systems.

R has quite a number of routines that are useful for statistical inference using nonlinear models. For example, the Optim function implements general-purpose optimization based on Nelder–Mead?, quasi-Newton and conjugate-gradient algorithms. It also includes an option for box-constrained optimization and simulated annealing. The nlme package implements Gaussian nonlinear mixed effects models, including models with nested random effects or spatial correlation, while the lme4 and glmmML packages offer options for generalized linear mixed-effect models. The pscl and ZIGP packages fit zero-inflated models appropriate for count data. The MCMCpack package allows the user to perform Bayesian inference via simulation from the posterior distributions of a variety of pre defined models. The metrop function implements Markov chain Monte Carlo for continuous random vector using a Metropolis algorithm.

AD Model Builder
AD Model Builder (ADMB) was specifically designed for complex highly-parameterized nonlinear models. ADMB uses automatic differentiation (AD) to provide the function optimizer with exact derivatives. The AD method, which is based on the chain rule, used by ADMB performs much more quickly and accurately than a numerical derivative calculation based on finite differences. ADMB plays a major role as a software environment for fisheries stock assessment. For example, the ADMB based general model Stock Synthesis (Methot 2009) is use to conduct the assessments of most major fisheries assessments in the United States and is also used in many other countries. ADMB also includes a Bayesian MCMC algorithm, likelihood profile calculations, simulation capabilities, and extensive matrix algebra functionality. A recent extension of ADMB adds the capability of modeling random effects, using the Laplace approximation or importance sampling (Skaug and Fournier 2006). ADMB is based on C++coding with additional structure such as keywords, functions, and operators. Knowledge of C++ is not a prerequisite for a successful ADMB user, because a simplified coding environment shields a user from advanced C++ concepts.

Potential Applications
The following are example applications that could be used by the working group.

Zero-inflated negative binomial mixed effects GAM model.
The example data is for silky shark bycatch data collected from the purse-seine fishery for tunas in the eastern Pacific Ocean. A study of these data has been published previously (Minami et al. 2007). The data set spans 11 years, and has over 32,000 observations. The purpose of the analysis was to estimate a standardized trend in bycatch-per-set for the silky shark, taking into consideration spatial, environmental and gear effects on shark bycatch rates. The bycatch data are counts of sharks, and contain both a large percentage of zeros as well as some very large values. The model was fitted in R, using the mgcv and MASS libraries, plus additional R code. The first few iterations of the fitting procedure used the EM algorithm to obtain good initial values. However, because the EM algorithm can be slow to converge, after obtaining reasonable initial values, the fitting procedure switched to a quasi-Newton method using the optim function of the MASS library to get faster convergence.

Population dynamics model with covariates.
Two recent NCEAS working group papers (Thomson et al. 2010, MacNally? et al. 2010) suggest using a life cycle model to evaluate the impact of different factors on delta smelt. Such an application includes a nonlinear model that is autoregressive and includes random effects (i.e. a state-space model). The data for this application are available directly from an existing NCEAS working group. Alternatively, data for Pacific herring of Prince William Sound could be used (Deriso et al. 2008).

Meta-population model.
Several data sets are available for the Hawaiian population of black-footed albatross (count, bycatch mortality, and banding) and have been used to simultaneously model the dynamics of its sub-populations to determine the impact of fisheries on the population. The model has been developed in ADMB.

Mark recapture model with random effects for survival and the probability of recapture.
This model has been implanted in both WinBUGS (Gimenez et al. 2008) and ADMB (Maunder et al. 2008). Many data sets are available for this type of application including the data sets provided in Maunder et al. (2008) and Gimenez et al. (2008). Other possible mark recapture applications include those that use covariates or individual random effects.

Spatial Poisson regression.
Saracco et al. (2010) uses a CAR model to examine the spatial variation in vital rates, but here we suggest an example that will model spatial variation in count data. This model has spatial random effects and many applicable datasets are available. Marc Kéry would provide a dataset of counts of particular bird species from the Swiss Atlas, covering a large nation-wide sampling area.

Dynamic occupancy model.
This model can be viewed as a type of zero-inflated binomial mixture model, but also as a discrete random effects model (see Royle and Kéry 2007). The dataset, provided by Marc Kéry, is derived from an opportunistic sampling scheme of observations made by volunteers (i.e., citizen scientists) and represents a growing trend in data collection.

Literature Cited
Much of the text comes directly from Maunder M.N., Schnute, J.T., and Ianelli, J. 2009. Computers in Fisheries Population Dynamics. In Megrey, B.A. and Moksness, E. eds.. Computers in Fisheries Research. Springer, pp: 337-372.

Other references include
Bolker, B.M. 2008. Ecological Models and Data in R. Princeton University Press, Princeton, NJ.

Bolker, B.M., Brooks, M.E., Clark, C.J., Geange, S.W., Poulsen, J.R., Stevens, M.H.H., and White, J.S. 2009. Generalized linear mixed models: a practical guide for ecology and evolution. Trends in Ecology and Evolution, 24(3):127–135.

Clifford, P. 1993. Discussion on the meeting on the Gibbs sampler and other Markov chain Monte Carlo methods. Journal of the Royal Statistical Society (Series B) 55:53–102.

Deriso, R.B., Maunder, M.N., and Pearson, W.H. 2008. Incorporating covariates into fisheries stock assessment models with application to Pacific herring of Prince William Sound, Alaska. Ecological Applications 18(5): 1270-1286.

Efron, B.E., Tibshirani, R.J. 1993. An Introduction to the Bootstrap.Chapman&Hall, NewYork?.
Fox, J. 2009. Aspects of the social organization and trajectory of the R project. The R Journal 1(2): 5-13.

Gimenez, O., Bonner, S., King, R., Parker, R.A., Brooks, S., Jamieson, L.E., Grosbois, V., Morgan, B.J.T., and Thomas, L. 2008. WinBUGS for population ecologists: Bayesian modeling using Markov Chain Monte Carlo methods. In: Modeling Demographic Processes in Marked Populations. Eds. Thomson, D.L., Cooch, E.G., and Conroy, M.J. Environmental and Ecological Statistics 3: 883-915.

Kéry, M. 2010. Introduction to WinBUGS for Ecologists. – A Bayesian approach to regression, ANOVA, mixed models and related analyses. Academic Press, Burlington, MA.

King, R., Morgan, B.J.T., Gimenez, O., and Brooks, S.P. 2009. Bayesian Analysis for Population Ecology. CRC Press.

Lunn, D., Spiegelhalter, D., Thomas, A. and Best, N. 2009. The BUGS project: Evolution, critique and future directions (with discussion), Statistics in Medicine 28: 3049--3082.

Mac Nally, R., Thomson, J.R., Kimmerer, W., Feyrer, F., Newman, K.B., Sih, A., Bennett, W., Brown, L., Fleishman, E., Culberson, S.D., Castillo, G. 2010. Analysis of pelagic species decline in the upper San Francisco Estuary using multivariate autoregressive modeling (MAR). Ecological Applications. Vol: 20(5). Pages 1417-1430.

Maunder M.N. 2004. Population Viability Analysis, Based on Combining Integrated, Bayesian, and Hierarchical Analyses. Acta Oecologica 26: 85-94. Special issue for the Extinction Working Group of the National Center for Ecological Synthesis and Analysis.

Maunder, M.N., Skaug, H.J., Fournier, D.A., and Hoyle, S.D. 2008. Comparison of estimators for mark-recapture models: random effects, hierarchical Bayes, and AD Model Builder. In: Modeling Demographic Processes in Marked Populations. Eds. Thomson, D.L., Cooch, E.G., and Conroy, M.J. Environmental and Ecological Statistics 3: 917-948.

McCarthy?, M.A., Keith, D., Tietjen, J., Burgman, M.A., Maunder M.N., Master, L., Brook, B.W., Mace, G., Possingham, H.P., Medellin, R., Andelman, S., Regan, H., Regan, T., and Ruckelshaus, M. 2004. Comparing predictions of extinction risk using models and subjective judgment. Acta Oecologica 26: 67-74. Special issue for the Extinction Working Group of the National Center for Ecological Synthesis and Analysis.

Minami, M, Lennert-Cody?, C.E., Wei, G. and Roman, M.H. 2007. Modeling shark bycatch: The zero-inflated negative binomial regression model with smoothing. Fisheries Research 84(2):210-221

Royle, J. A., and Kéry, M. 2007. A Bayesian state-space formulation of dynamic occupancy models. Ecology 88: 1813–1823.

Saracco, J.F., Royle, J.A., DeSante?, D.F., and Gardner, B. 2010. Modeling spatial variation in avian survival. Ecology 91:1885-1891.

Schnute, J.T., Richards, L.J., Olsen, N. 1998. Statistics, software, and fish stock assessment. In: Funk F, Quinn II TJ, Heifetz J, Ianelli JN, Powers JE, Schweigert JJ, Sullivan PJ, Zhang CI (eds) Fishery Stock Assessment Models. Alaska Sea Grant College Program Report No. AK-SG-98-01, University of Alaska Fairbanks, pp 171–184.

Seber, G. A. F., and C. J. Wild. 1989. Nonlinear regression. Wiley, New York, New York, USA.
Skaug, H. and Fournier, D. 2006. Automatic approximation of the marginal likelihood in nonlinear hierarchical models. Computational Statistics and Data Analysis 51:699–709

Thomson, J., Kimmerer, W., Brown, L., Newman, K. B., Mac Nally, R., Bennett, W., Feyrer, F., Fleishman, E. 2010. Bayesian change point analysis of abundance trends for pelagic fishes in the upper San Francisco Estuary. Ecological Applications. Vol: 20(5). Pages 1431-1448.