unmarked: An R Package for Fitting Hierarchical Models of Wildlife Occurrence and Abundance

Ecological research uses data collection techniques that are prone to substantial and unique types of measurement error to address scientific questions about species abundance and distribution. These data collection schemes include a number of survey methods in which unmarked individuals are counted, or determined to be present, at spatiallyreferenced sites. Examples include site occupancy sampling, repeated counts, distance sampling, removal sampling, and double observer sampling. To appropriately analyze these data, hierarchical models have been developed to separately model explanatory variables of both a latent abundance or occurrence process and a conditional detection process. Because these models have a straightforward interpretation paralleling mechanisms under which the data arose, they have recently gained immense popularity. The common hierarchical structure of these models is well-suited for a unified modeling interface. The R package unmarked provides such a unified modeling framework, including tools for data exploration, model fitting, model criticism, post-hoc analysis, and model comparison.


Imperfect detection and data collection in ecological research
A fundamental goal of ecological research is to understand how environmental variables influence spatial or temporal variation in species abundance or occurrence. Addressing these research questions is complicated by imperfect detection of individuals or species. Individuals may go undetected when present for a variety of reasons including their proximity to the observer, cryptic behavior, or camouflage. Thus, imperfect detection can introduce substantial measurement error and obscure underlying ecological relationships if ignored.  Table 1: An example of the data structure required by unmarked's fitting functions. Counts of organisms were made at M = 4 sites during T = 2 seasons with J = 3 visits ("observations") per season.
To accommodate imperfect detection of individuals or species, ecologists have developed specialized methods to survey wildlife populations such as site occupancy sampling, repeated counts, distance sampling, removal sampling, and double observer sampling (see Section 2 and Williams et al. 2002 for definitions). This myriad of sampling methods are all unified by a common repeated-measures type of sampling design in which J "observations" are made at M spatial sample units during each of T seasons. As an example, Table 1 shows a typical repeated count dataset.
In addition to the repeated-measures structure, there are two important features to note about the data shown in Table 1. First, this design is often referred to as a 'metapopulation design' because the population may be regarded as an aggregation of subpopulations and because the spatial sampling is an explicit component of the problem. This is important because it provides a basis for modeling variation among sites as a function of site-specific covariates, and because metapopulation parameters such as local colonization and extinction can be directly estimated. Second, the counts are of individuals that may not be uniquely recognized, so although it may be possible to determine if three individuals were detected on two different visits (as was the case at Site 1 in Table 1), it is not always possible to determine if they were the same three individuals.
Both the metapopulation design and absence of individual recognition distinguish these data from data collected using another large suite of ecological sampling methods known as 'capture-recapture' methods. A vast number of models have been developed for capturerecapture data, and many of them are implemented in the software program MARK (White and Burnham 1999). The origins of unmarked (the package and its name) arose to fill the void of models for data from studies of unmarked individuals involving explicit spatial sampling.

Hierarchical models for metapopulation designs
While many ecological studies produce data consistent with a metapopulation design, the development and implementation of models for inference under individual sampling protocols has proceeded in a piecemeal fashion, without a coherent analysis platform. Recently, however, a broad class of hierarchical models (see Royle and Dorazio 2008 for a general treatment) has been developed that offer a unified framework for analysis by formally recognizing that observations are generated by a combination of (1) a state process determining abundance or species occurrence at each site and (2) a detection process that yields observations conditional on the state process. The model for the state process describes abundance or occurrence at each site, but due to imperfect detection, these quantities cannot be observed directly and are regarded latent variables. This paper introduces unmarked, an R (R Development Core Team 2011) package that provides a unified approach for fitting this broad class of hierarchical models developed for sampling biological populations. It is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=unmarked. Inference is based on the integrated likelihood wherein the latent state variable is marginalized out of the conditional likelihood.
1.3. Scope and features of unmarked unmarked provides tools to assist researchers with every step of the analysis process, including data manipulation and exploration, model fitting, post-hoc analysis, model criticism, and model selection. unmarked provides a growing list of model-fitting functions designed for specific sampling methods. The fitting functions each find the maximum likelihood estimates of parameters from a particular model (Section 3.3) and return an object that can be easily manipulated. Methods exist for performing numerous post-hoc analyses such as requesting linear combinations of parameters, back-transforming parameters to constrained scales, determining confidence intervals, and evaluating goodness of fit. The model specification syntax of the fitting functions was designed to resemble the syntax of R's common fitting functions such as lm for fitting linear models.
Although there is existing software for fitting some of these models (e.g., Hines and MacKenzie 2002), there are a number of advantages to a unified framework within R. Many researchers are already familiar with R and use its powerful data manipulation and plotting capabilities. Sometimes many species are analyzed in tandem, so that a common method of aggregating and post-processing of results is needed, a task easily accomplished in R. Unlike other available software, unmarked makes it possible to map habitat-specific abundance and species distributions when combined with R's GIS capabilities. Another important advantage of unmarked's approach is that researchers can simulate and analyze data within the same computational environment. This work flow permits simulation studies for power analysis calculations or the effectiveness of future sampling designs. All of this is made much simpler by analyzing the data within R, and using a single environment to complete all phases of the analysis is much less error-prone than switching between applications.
In this paper, Section 2 gives a brief summary of many of the models unmarked is capable of fitting. Section 3 describes general unmarked usage aided by a running data example.

Models implemented in unmarked
The list of models implemented in unmarked continues to grow as new models are developed. Table 2 shows the models available as of version 0.9-0. Rather than describe each fitting function in detail, this section provides a summary of several of the most common sampling techniques and how unmarked can be used to model the resulting data.

Site occupancy models
Single season site occupancy model  otherwise. Much interest is focused on estimating functions of species occurrence (e.g., number of occupied sites) or to identify factors that are associated with changes in the probability of a site being occupied, i.e., ψ = Pr(Z i = 1). To estimate these parameters, researchers employ a sampling design, whereby surveyors visit a sample of M sites and record the binary response Y ij of species detection (Y = 1) or non-detection (Y = 0) during j = 1, . . . , J i visits to the ith site during a 'season' . A key feature of species occurrence surveys is that false absences are possible, so that a species might go undetected (Y ij = 0) even if it is present (Z i = 1). The detection probability parameter p accounts for this observation process, and is defined as the probability of detecting a species that is present.
The key assumptions made when modeling these data are that the occupancy state at a site remains constant throughout the season and repeated visits at a site are independent. In order to meet the first assumption, a season will generally be a very short time frame, such as a few months during a breeding season, or even a few minutes if repeated visits are made in quick succession. Replicate samples (i.e., J > 1) provide information about the detection rate separate from the occupancy rate.
The following hierarchical model describes the joint distribution of the observations conditional on the latent occupancy state, and the marginal distribution of the latent occupancy state variable: Removing the latent Z variables by marginalization yields the likelihood: where I(.) is the indicator function taking the value 1 if its argument is true, and 0 otherwise.
The model and likelihood are easily generalizable to accommodate covariates on detection probability or occupancy. Variables that are related to the occupancy state are modeled as where x i is a vector of site-level covariates and β is a vector of their corresponding effect parameters. Similarly, the probability of detection can be modeled with where v ij is a vector of observation-level covariates and α is a vector of their corresponding effect parameters. Examples of site-level covariates include habitat characteristics such as vegetation height. Observation-level covariates could include time of day, date, wind speed, or other factors that might affect detection probability. The function occu implements this model.

Multi-season site occupancy model
Sometimes the study objective is to understand the dynamics of the occupancy state over time.
To obtain such information researchers conduct repeated occupancy studies (see above) at the same sample of sites over consecutive seasons (MacKenzie et al. 2003) and seek to estimate probabilities of colonization (γ) and extinction ( ), where colonization is the change of an unoccupied site to occupied and extinction if the change of an occupied site to unoccupied. If the occupancy status is assumed to evolve according to a Markov process, then a 2-state finite hidden Markov model describes these data. Let Y ijt denote the observed species occurrence status at visit j during season t to site i. Then This model generalizes the single season site occupancy model by relaxing the population closure assumption. To define the likelihood of this model, let φ 0 = (ψ, 1 − ψ) and Furthermore, let the matrix of one-step Markov transition probabilities be The likelihood is then given by which can be maximized by the colext function in unmarked. Again, each of the four model parameters ψ, γ, , and p can be modeled as functions of covariates on the logit scale.

Repeated count data
Occupancy is a crude summary of population structure and dynamics and in practice considerable effort is focused on obtaining information about abundance (Dorazio 2007). A common sampling design for obtaining abundance information, analogous to the basic occupancy design described above, is to repeatedly visit a sample of M sites J times and record the number of unique individuals observed at each site. Similar assumptions are made as with occurrence data: (1) abundance at a site remains constant during a season and (2) counts at a site are independent. Royle (2004b) presented the following hierarchical model for repeated count data. Let N i be the unobserved total number of individuals using a site and define Y ij as the number of individuals observed during the jth visit. Then, where λ i is the abundance rate at site i and p is the detection probability. f is a discrete distribution with support restricted to N i ≥ 0 and θ are extra parameters of f other than the location parameter, λ. unmarked currently supports f as Poisson or negative binomial. In the negative binomial case, θ is a dispersion parameter, which is useful when overdispersion is suspected. In the Poisson case, the overdispersion parameter θ is constrained to equal zero.
As with the occupancy model, covariates may be included at either the state (here, abundance) or detection levels, but abundance is modeled through a log link to enforce its positivity constraint.
where x i is a vector of site-level covariates and β is a vector of their corresponding effect parameters. Similarly, the probability of detection can be modeled with where v ij is a vector of observation-level covariates and α is a vector of their corresponding effect parameters.
The N i variables are latent and so analysis is carried out based on the integrated likelihood obtained by marginalizing each N i from the conditional likelihood (Royle et al. 2004): This likelihood can be maximized using the pcount function.
The function pcountOpen implements a generalized form of this model designed for 'open' populations with temporal dynamics governed by apparent survival and recruitment parameters (Dail and Madsen 2011). Another related model was described by Royle and Nichols (2003) and can be fit with the occuRN function. This model estimates abundance from site occupancy data by exploiting the link between abundance and detection probability.

General multinomial-Poisson mixture model
Here we discuss a more general class of models that can be customized by the user for a variety of sampling methods, the multinomial-Poisson model (Royle 2004a). The general form of this model is where N i is the latent abundance at site i as with the repeated count model, and π = (π 1 , π 2 , . . . , π J ) is the vector of cell probabilities corresponding to the vector of counts Y i . Examples of sampling methods that produce data of this form include removal sampling, double observer sampling, and distance sampling. In the next sections, we discuss these specific cases. In general, π is determined by the specific sampling method and j π j ≤ 1 because detection is imperfect and π J+1 = 1 − j π j is the probability of the species escaping detection.
Note that the replication here is of a different nature than in previously described modelsnot over time necessarily -but still effectively replication as far as the design goes. That is, there are repeated measurements at each site, but with a multinomial protocol, the replicate counts are dependent instead of independent.
To illustrate the likelihood under a multinomial observation model, suppose that a sampling method produces a multinomial observation with 3 observable frequencies at each site, ie J = 3. As before, N i are latent variables, and so inference is based on the integrated likelihood of Y i which is: For this specific case of a multinomial-Poisson mixture, the likelihood simplifies analytically (to a product-Poisson likelihood). The function multinomPois can be used to maximize this likelihood.
For multinomial-Poisson sampling methods, the actual observations are an underlying categorical detection variable with M ≤ J levels so that the J-dimensional Y i is derived from the M -dimensional raw counts in some sampling method-specific manner. Thus, it is necessary to model the detection at the raw observation level, denoted p ik for k = 1, 2, . . . , M at site i. Then we derive the multinomial cell probabilities π i through the sampling technique-specific relationship π ij = g(p ik ) where p ik is the underlying probability of detection and g is some sampling method-specific function.
Thus, the only two requirements to adapt unmarked's general multinomial Poisson modeling to a new sampling method is to specify g and a binary 0/1 matrix that describes the mapping of elements of p i = (p i1 , . . . , p iR ) to elements of π i . This mapping matrix, referred to in unmarked as obsToY, is necessary to consistently clean missing values from the data and relate observation-level covariates with the responses. The (j,k)th element of obsToY is 1 only if p ik is involved in the computation of π ij . The detection function g is called piFun in unmarked.
Covariates may be included in either the state (here, abundance) or detection models, through where x i is a vector of site-level covariates and β is a vector of their corresponding effect parameters. Similarly, the probability of detection can be modeled with where v ij is a vector of observation-level covariates and α is a vector of their corresponding effect parameters.
We now describe two common sampling methods that can be modeled with the multinomial-Poisson model: removal sampling and double observer sampling. These two methods are included in unmarked, but additional methods, such as mark-recapture samples, may easily be specified by defining a user-specified piFun function and obsToY matrix.
Removal sampling Popular in fisheries, removal sampling is commonly implemented by visiting a sample of M sites J times each and trapping or otherwise capturing and then removing individuals at each visit. Thus, Y ij is the number of individuals captured at the jth visit for j = 1, 2, . . . , J. We can specify g for removal sampling as follows. The probability of an individual at site i being captured on the first visit is π i1 = p i1 . The probability of capture on the jth visit is for j = 2, . . . , J and the probability of not being sampled is Thus, the mapping matrix is an J × J matrix with ones in the upper triangle, Double observer sampling Double observer sampling involves collecting data by a team of two surveyors simultaneously visiting a site. Each observer independently records a list of detected organisms, and at the end of the survey the two observers attempt to reconcile their counts. If individuals are not uniquely marked, this may be a difficult task in practice; however, assuming that individuals can be distinguished, the data at each site are a vector of length three (Y i ), corresponding to the numbers of individuals seen by only observer one, only observer two, and both observers. Thus, for double observer sampling, g is defined as follows The obsToY mapping matrix for double observer sampling is the following 3 × 2 matrix Distance sampling One of the most widely used sampling methods in animal ecology is distance sampling (Buckland et al. 2001), which involves recording the distance to each individual detected at M sites (often referred to as 'transects') on a single occasion. Detection probability is modeled as a function of distance (d) to the observer, for example using the halfnormal detection function p = exp(−d 2 /(2σ 2 )) where σ is the half-normal shape parameter. In practice, the distance measurements are often binned into J distance intervals, which allows them to modeled as multinomial outcomes with cell probabilities π i computed as the product of the probability of detection and the probability of occurrence in each distance interval (Royle et al. 2004).
As currently implemented, the general form of the distance sampling model is identical to the multinomial-Poisson mixture described above, with the sole difference being that instead of Here, the log link is required because σ is a positive shape parameter of the detection function. In addition, distance sampling data are associated with many unique attributes such as distance interval cut-points and survey method (line-transect vs point count); therefore, we created the specialized function distsamp to fit the multinomial-Poisson model to distance sampling data.
Generalization The flexibility of the multinomPois function is extended even further in the gmultmix function, which allows for a negative binomial mixture (Dorazio et al. 2005) and relaxes the population closure assumption (Chandler et al. 2011).

unmarked usage
unmarked provides data structures, fitting syntax, and post-processing that form a cohesive framework for the analysis of ecological data collected using a metapopulation design. In order to achieve these goals, unmarked uses the S4 class system (Chambers 2008). As R's most modern system of class-based programming, S4 allows customization of functions, referred to as methods, to specific object classes and superclasses. For example, when the generic predict method is called with any unmarked model fit object as an argument, the actual predict implementation depends on the specific model that was fit. Use of class-based programming can provide more reliable and maintainable software while also making the program more user-friendly (Chambers 2008).

Preparing data
unmarked uses a custom S4 data structure called the unmarkedFrame to store all data and metadata related to a sampling design. Although this at first appears to add an extra layer of work for the user, there are several reasons for this design choice. The multilevel structure of the models means that standard rectangular data structures such as data.frames or matrices are not suitable for storing the data. For example, covariates might have been measured separately at the site level and at the visit level. Furthermore, the length of the response vector Y i at site i might differ from the number of observations at the site as in the multinomial Poisson model. In some cases, metadata of arbitrary dimensions may need to be associated with the data. For example, in distance sampling it is necessary to store the units of measure and the survey design type. Aside from these technical reasons, Gentleman (2009) pointed out that the use of such portable custom data objects can simplify future reference to previous analyses, an often neglected aspect of research. Repeated fitting calls using the same set of data require less code repetition if all data are contained in a single object. Finally, calls to fitting functions have a cleaner appearance with a more obvious purpose when the call is not buried in data arguments.
The parent S4 data class is called an unmarkedFrame and each unmarked fitting function has its own data type that extends the unmarkedFrame. To ease data importing and conversion, unmarked includes several helper functions to automatically convert data into an unmarkedFrame: csvToUMF which imports data directly from a comma-separated value text file, formatWide and formatLong which convert data from data frames, and the family of unmarkedFrame constructor functions.
An unmarkedFrame object contains components, referred to as slots, which hold the data and metadata. All unmarkedFrame objects contain a slot for the observation matrix y, a data.frame of site-level covariates siteCovs, and a data.frame of observation-level covariates obsCovs. The y matrix is the only required slot. Each row of y contains either the observed counts or detection/non-detection data at each of the M sites. siteCovs is an M -row data.frame with a column for each site-level covariate. obsCovs is an M J-row data.frame with a column for each observation-level covariate. Thus each row of obsCovs corresponds to a particular observation, with the order corresponding to site varying slower and observation within site varying faster. Both siteCovs and obsCovs can contain NA values corresponding to unbalanced or missing data. If a site-level covariate is missing, unmarked automatically removes all data for that site prior to fitting the model. Missing values in the obsCovs are handled by removing the corresponding occurrence or count observations such that the missing values make no contribution to the likelihood. unmarked provides constructor functions to make creating unmarkedFrames straightforward. For each specific data type, specific types of unmarkedFrames extend the basic unmarkedFrame to handle model-specific nuances.

Importing repeated count data
Here is an example of creating an unmarkedFrame for repeated count data (Section 2.2). First, load Mallard (Anas platyrhynchos) point count dataset described in Kéry et al. (2005).

R> library("unmarked") R> data("mallard")
Loading the mallard data makes three objects available within the R workspace. The matrix mallard.y contains the number of mallards counted at each of M = 239 sites (rows) on J = 3 visits (columns The site-level covariates are columns of the mallard.site data.frame, which also has M = 239 rows. The site-level covariates are elevation (elev), transect length (length), and the proportion of forest covering the site (forest).
The observation-level covariates are a list named mallard.obs with separate M ×J matrices for each observation-level covariate. Here, the two observation-level covariates are a measure of survey effort (ivel) and the date of the survey (date). Both have been standardized to a mean of zero and unit variance.  The unmarkedFrame constructors can accept obsCovs in this list format or as a data.frame in the format described in Section 3.1.
The following call to unmarkedFramePCount organizes the observations and covariates into an object that can be passed to the data argument of the fitting function pcount. The summary reveals that only 40 sites have at least one detection. Note that the "number of observations" in this case refers to number of visits. The term "observation" is used rather than "visit" because the definition of the J replicate surveys depends upon the sampling method. The tabulation of y observations provides additional evidence of sparse counts, with no mallards being detected during 576 of the surveys. The 58 NA values correspond to missing data.

Importing removal sampling data
To illustrate the slightly different syntax for removal sampling data example, we will import data from a removal survey of Ovenbirds (Seiurus aurocapillus) described by Royle (2004a).
The data consist of a list named ovendata.list with a matrix named data containing the removal counts for 4 visits and a data.frame called covariates containing site-level covariates information.

R> data("ovendata")
This loads a list named ovendata.list with components for the observed count data and covariates. The only additional specification required when creating an unmarkedFrame for the multinomial-Poisson model is to specify the particular type of data as either removal, double, or userDefined.

Preparing multi-season occupancy data
The data structures are more complex when surveys were conducted over more than one season. The following example demonstrates how to prepare a toy dataset for the colext fitting function. This detection/non-detection dataset has dimensions M = 4 sites, T = 2 seasons, and J = 3 visits per season. The function unmarkedMultFrame accepts siteCovs and obsCovs like all other unmarked-Frame classes, and it also has an argument yearlySiteCovs that can accept a list of M × T data.frames containing season-level covariates. The numPrimary argument specifies the number of seasons.

Manipulating unmarkedFrames
The various components of unmarkedFrames can be extracted and subsetted in a manner similar to the methods used to manipulate standard R objects. Subsetting can be accomplished using the 'bracket' notation. For example, the first five rows of data and all three removal occasions can be extracted from the Ovenbird removal study using In some cases, the covariate data may need to be manipulated after the unmarkedFrame has been created. The following code demonstrates how to extract the site-level covariates, standardize those that are continuous (columns 2 and 3), and reinsert them back into the unmarkedFrame.

Fitting models
As introduced in Section 2, each type of survey design has a corresponding fitting function. For example, to fit a repeated count model, we call pcount. Table 2 provides a summary of all models that unmarked currently fits. With the exception of the 'open' population models (colext, gmultmix, and pcountOpen), all fitting functions use a double right-hand sided formula syntax representing the hierarchical model structure. Specifically, covariates affecting the detection process are specified following the first tilde, and covariates of the state process follow the second tilde. No left-hand side of the formula is specified because the unmarkedFrame defines the response variable uniquely as the y slot.

Fitting a repeated count model
Continuing the Mallard example, the following call to pcount fits a binomial-Poisson mixture model (Section 2.2). The following code specifies that detection probability p should be modeled by day of year, including a quadratic term. We also wish to model abundance using elevation and proportion of area forested. As described in Section 2.2, covariates of detection are on the logit-scale and covariates of abundance are on the log-scale for the repeated count model.
R> fm.mallard.1 <-pcount(~date + I(date^2)~elev + forest, + data = mallardUMF, K = 50) R> fm.mallard.1 This initial fit suggests that Mallard abundance decreases with increasing elevation and forest. It also looks like a linear model might suffice for the detection model, so we subsequently fit the linear detection model as follows: R> fm.mallard.2 <-pcount(~date~elev + forest, data = mallardUMF, + K = 50) R> fm.mallard.2 This seems to be a better model according to both the Wald p value and AIC. The result suggests that detection probability decreases during the course of a year.

Fitting a multinomial-Poisson model
Here we demonstrate fitting a multinomial-Poisson mixture model to removal sampling data. The Ovenbird data has no observation-level covariates, so detection probability is assumed constant across visits. It is not necessary to specify that removal sampling was used when fitting the model because this information is already stored in the ovenFrame data. We model abundance as a function of understory forest coverage (ufp) and average basal tree area (trba).

Examining fitted models
Objects returned by unmarked's fitting functions also make use of the S4 class system. All fitted model objects belong to the unmarkedFit parent class. Thus, common operations such as extracting coefficient estimates, covariance matrices for estimates, and confidence intervals have been adapted to behave similar to R's base functions.
For example, we can extract estimated coefficients either from the entire model, or from the state or detection levels by specifying the type argument. Similarly, the vcov function extracts the covariance matrix of the estimates, using the observed Fisher information by default. vcov also accepts a type argument, as does the convenience method SE, which returns standard errors from the square root of the diagonal of the covariance matrix.
R> vcov(fm.mallard.2, type = "det") p(Int) p(date) p(Int) 0.03549024 0.00344853 p(date) 0.00344853 0.01302865 Extracting confidence intervals proceeds in a similar fashion. By default, the asymptotic normal approximation is used. Profile confidence intervals are also available upon request. This can take some time, however, because for each parameter, a nested optimization within a root-finding algorithm is being used to find the profile limit. The profile confidence intervals and normal approximations are quite similar here.

The nonparametric bootstrap
Nonparametric bootstrapping can also be used to estimate the covariance matrix. unmarked implements a two-stage bootstrap in which the sites are first drawn with replacement, and then within each site, the observations are drawn with replacement. First, bootstrap draws must be taken using the nonparboot function, which returns a new version of the unmarkedFit object with additional bootstrap sampling information. Thus, this new fit must be stored, either in a new fit object or the same one, and then subsequently queried for bootstrap summaries. To illustrate, we use the removal sampling data instead of the repeated count data because computations are much faster; however, bootstrapping is available for any of the models in unmarked.

Linear combinations of estimates
Often, meaningful hypotheses can be addressed by estimating linear combinations of coefficient estimates. Linear combinations of coefficient estimates can be requested with linearComb.
Continuing the Ovenbird example, the following code estimates the log-abundance rate for a site with ufp = 0.5 and trba= 0.

Linear combination(s) of Abundance estimate(s)
Estimate SE (Intercept) ufp trba 0.153 0.13 1 0.5 0 Multiple sets of coefficients may be supplied as a design matrix. The following code requests the estimated log-abundance for sites with ufp = 0.5 and trba = 1.

Back-transforming linear combinations of coefficients
Estimates of linear combinations back-transformed to the native scale are likely to be more interesting than the direct linear combinations. For example, the logistic transformation is applied to estimates of detection rates, resulting in a probability bound between 0 and 1. This is accomplished with the backTransform. Standard errors of back-transformed estimates are estimated using the delta method. Confidence intervals are estimated by back-transforming the confidence interval of the original linear combination.

Model selection
unmarked performs AIC-based model selection for structured lists of unmarkedFit objects. To demonstrate, we fit a few more models to the Ovenbird removal data, including an interaction model, two models with single predictors, and a model with no predictors.
Next, we organize the fitted models with the fitList function and the use the modSel method to rank the models by AIC. predict functions much like linearComb except that new data can be passed to it as a data.frame rather than a design matrix. When the first argument given to predict is a list of models created by fitList, predict computes model-averaged predictions, which may be useful in the presence of high model selection uncertainty.

Goodness of fit and the parametric bootstrap
To conduct goodness of fit tests, unmarked provides a generic parametric bootstrapping function parboot. It simulates data from the fitted model and applies a user-defined function that returns a fit-statistic such as the Pearson's χ 2 . R> set.seed(1234) R> chisq <-function(fm) { + observed <-getY(fm@data) + expected <-fitted(fm) + sum((observed -expected)^2/expected) + } R> pb <-parboot(fm.oven.1, statistic = chisq, nsim = 200) R> plot(pb, main = "") The above call to plot with a parametric bootstrap object as the argument produces a useful graphic for assessing goodness of fit ( Figure 2). The plot suggests that the model adequately explains these data.
Beyond serving as a tool to evaluate goodness of fit, parboot can be used to characterize uncertainty in any derived quantity of interest.

Future directions for unmarked development
unmarked has become a stable and useful platform for the analysis of ecological data, but several areas of development could improve its utility. First, new models need to be added to cover the range of sampling techniques and population dynamics commonly encountered.

Sampling method Population dynamics Closed
Open to Open to temporary demographic emigration processess Occurrence sampling occu colext Repeated counts pcount -pcountOpen Removal sampling, multinomPois gmultmix double observer sampling, and other multinomial designs Distance sampling distsamp -- Table 3: Model fitting functions classified by sampling method and population dynamics.
Missing cells indicate models that have not been developed but are likely to be investigated in the future. Table 3 illustrates some of the gaps that need to be filled. In most cases, models to fill these gaps have not been developed so more research is needed.
Second, each of the models in unmarked assumes independence among sites. However, ecologists often use sampling methods such as cluster sampling that induce spatial dependence. Typically, this is done for logistical convenience, but because few methods are available to account for spatial correlation and imperfect detection probability, the spatial dependence is often ignored. Rather than this being a weakness of the sampling design, we envision that this dependence can be used as information regarding the spatial distribution of individuals.
Third, many of the likelihoods are written in pure R, which can be slow for large problems.
We are currently translating many of these functions into C++ with the help of the R package Rcpp (Eddelbuettel and François 2011).
Finally, Markov chain Monte Carlo methods could be implemented for all of these models allowing for Bayesian inference. An important advantage of Bayesian analysis over classical methods is that the latent abundance or occurrence variables can be treated as formal parameters. Thus posterior distributions could easily be calculated for derived parameters such as the proportion of sites occupied. Bayesian analysis also would provide a natural framework for incorporating additional sources of random variation. For example, one could model heterogeneity among sites not accounted for by covariates alone.