Bayesian model selection for spatial capture-recapture models.

A vast amount of ecological knowledge generated over the past two decades has hinged upon the ability of model selection methods to discriminate among various ecological hypotheses. The last decade has seen the rise of Bayesian hierarchical models in ecology. Consequently, commonly used tools, such as the AIC, become largely inapplicable and there appears to be no consensus about a particular model selection tool that can be universally applied. We focus on a specific class of competing Bayesian spatial capture-recapture (SCR) models and apply and evaluate some of the recommended Bayesian model selection tools: (1) Bayes Factor-using (a) Gelfand-Dey and (b) harmonic mean methods, (2) Deviance Information Criterion (DIC), (3) Watanabe-Akaike's Information Criterion (WAIC) and (4) posterior predictive loss criterion. In all, we evaluate 25 variants of model selection tools in our study. We evaluate these model selection tools from the standpoint of selecting the "true" model and parameter estimation. In all, we generate 120 simulated data sets using the true model and assess the frequency with which the true model is selected and how well the tool estimates N (population size), a parameter of much importance to ecologists. We find that when information content is low in the data, no particular model selection tool can be recommended to help realize, simultaneously, both the goals of model selection and parameter estimation. But, in general (when we consider both the objectives together), we recommend the use of our application of the Bayes Factor (Gelfand-Dey with MAP approximation) for Bayesian SCR models. Our study highlights the point that although new model selection tools are emerging (e.g., WAIC) in the applied statistics literature, those tools based on sound theory even under approximation may still perform much better.

INTRODUCTION

Following the highly influential text, Burnham and Anderson (2002), on model selection, ecologists and conservation biologists have drastically shifted their inferential practice from the “hypothesis testing” approach to the more appropriate “hypothesis discrimination” approach (Johnson & Omland, 2004). To date, Burnham and Anderson (2002) has been cited 44,168 times (as on January 3, 2019—Google Scholar), demonstrating the impact of this text. One may argue that this contribution has helped increase the pace of growth in ecological knowledge because it has paved the way for researchers to draw inferences more robustly because of the ability to now assess the influence of various competing a priori hypotheses (models) without altering the study question to suit the restrictive hypothesis testing paradigm (Bolker, 2008).

Since a large amount of ecological data are based on field observations, it called for ecologists to take the approach of “detectives” rather than “hypothetico‐deductive” scientists by formulating models using likelihood functions to confront various a priori hypotheses using observational data (Hilborn & Mangel, 1997). And by maximizing the likelihood and using a model selection tool, such as the Akaike's information criterion (Burnham & Anderson, 2002), researchers found a way to place increased faith on models favored by such criteria. Thus, a vast amount of ecological knowledge generated has relied on the robustness of such model selection tools in accurately discriminating hypotheses.

Recently, there has been an increased use of hierarchical models in ecology since they appear to address two important issues: (a) ecological scales are naturally hierarchical in structure and (b) hierarchical models form a natural way of incorporating the observation process (Royle & Dorazio, 2008). With powerful tools such as the MCMC, it is now possible to confront complex ecological models with data in a Bayesian inferential framework (Bolker, 2008). However, it remains unclear as to how to discriminate among competing hypotheses (models) because popular model selection criteria (such as AIC, BIC, or DIC; Burnham & Anderson, 2002) are not easy to apply or work poorly for complex hierarchical models (Millar, 2009).

It is well known that, asymptotically, the Bayes factor is the preferred model selection tool due to its consistency property, that is, to identify the true data‐generating model if (and only if) the true data‐generating model is included in the model set and if the data tend to the limit of infinite informativeness (Ghosh, Delampady, & Samanta, 2006; Robert, 2007). However, this property holds only under certain regularity conditions that are often difficult to verify for complex models (Berger, Ghosh, & Mukhopadhyay, 2003; Dass & Lee, 2004; Ghosh & Samanta, 2001). More prominently, there is vast literature expressing the difficulties in computing the marginal likelihood in applied statistical problems (Chan & Eisenstat, 2015; Wang, 2018). Hence, it becomes necessary to also consider alternatives to the Bayes factor or to find novel ways of applying them in practice.

Recently, Hooten and Hobbs (2015) summarized a wide array of Bayesian model selection methods that are available to ecologists. However, the generality of the recommendations provided by them remains unknown. Given such innate difficulties involved in discovering the “ideal” model selection tool both from the standpoint of theory and its application to a broad class of models, it appears to be prudent to explore the model selection issue by conditioning, at least, on a particular class of models.

Here, we evaluate various Bayesian model selection tools on a class of Bayesian spatial capture–recapture (SCR) models that are now used frequently for animal density estimation (Royle, Chandler, Sollmann, & Gardner, 2013). Although, previously, Goldberg et al. (2015) has attempted to apply the Bayes Factor (Gelfand‐Dey estimator) in an abundance estimation problem for leopards (Panthera pardus), their approach of computing the ratio term in the estimator seems inaccurate in the context of how the denominator has to be computed according to Gelfand and Dey (1994). Thus, we evaluate various Bayesian model selection tools by: (a) defining a class of competing models (in our case, these include the model developed in Dey, Delampady, Karanth, and Gopalaswamy (2019) along with simplified alternatives) that vary both in terms of structural and model complexity (b) simulating data sets from a “true” model (c) practically implementing a variety of Bayesian model selection tools, and in specific cases, also proposing alternatives previously not defined (d) assessing the efficacy of these implementations from the standpoint of model selection and parameter estimation and (e) providing recommendations to practitioners based on our results.

METHODS

We describe here the sampling design and development of the competing models in the candidate model set.

The candidate model set

Capture–recapture surveys are conducted by placing an array of detectors (e.g., camera traps, hair snares) to sample the species of interest within a bounded region over a fixed period of time. As an extension, spatial capture–recapture (SCR) models draw inference on the spatial distribution of animals using their spatial locations from the recorded capture–recapture samples (Royle et al., 2013).

In photographic capture–recapture surveys, an array of camera traps are placed over the study area. Each camera trap station consists of two cameras placed opposite to one another to photograph the flanks of animals passing by. Naturally marked species can be identified by the unique patterns on their flanks. However, an animal passing through a camera trap station does not necessarily result in identifiable flank photographs from both the sides. This is because there are many known and unknown factors that can influence camera firing rates. Such a detection process results in uncertain or partially identified individuals in the spatial capture—recapture sample and thus provided the motivation for the development of Dey et al. (2019). We summarize the description of the model below.

Sampling situation

The notation used in this article is described in Tables 1 and 2. However, we describe a few variables and parameters for ease in the model description below. Consider a capture–recapture survey of a species with naturally marked individuals in which two detectors are collocated at J trap stations (within a bounded geographic region ) and kept active for K sampling occasions. An individual can be completely identified if both the detectors record the individual simultaneously at least once during the course of study (Royle, 2015). We assume that each detector captures some mutually exclusive attributes of an individual. These capture outcomes are recorded as binary observations and for an individual i at trap station x on sampling occasion k corresponding to detectors 1 and 2, respectively. The paired Bernoulli outcomes give rise to bilateral spatial capture–recapture data for each individual i at location x on occasion k. The array of a bilateral capture history for an individual i is denoted by , which is of dimension 2 × J × K. In Example 2.1, we provide an example of a sample data set coming out of a spatial capture–recapture survey with two detectors deployed at each station.

Table 1

Notations of variables and parameters used in this article

Variables and parameters	Definition
V	A bounded geographic region of scientific or operational relevance where a population of individuals of certain species reside
N ~ Binomial (M, ψ)	Population size of the superpopulation, that is, the number of individuals within V
M	Maximum number of individuals within the state space V This is a fixed quantity defined by the investigator
ψ	Proportion of individuals that are real and present within V
θ	Probability that an individual is male
J	Number of trap stations in V
K	Number of sampling occasions
R	Maximum permissible value of movement range for each individual during the survey
ω 0	Baseline trap entry probability in the models M 1 and M 2, that is, probability that an individual passes through a trap station assuming its center of activity is also located at that trap station
p 0	Baseline detection probability in the models M 2 and M 4, that is, probability that an individual is detected by a detector assuming its center of activity is also located at that trap station
σ	σ measures the spatial extent of movement around individual activity center. σ = σm for male individuals, σ = σf for female individuals
d_ij=d(s_i,x_j)=s_i-x_j	Euclidean distance between points s i and x j
η_j(s_i,u_i)=ω_0exp-d(s_i,x_j)^22σ(u_i)^2	Probability that an individual i passes through a trap station x j on some occasion k and σ is modeled as a function of individual covariate on sex category ui
η_j(s_i)=ω_0exp-d(s_i,x_j)^22σ^2	Probability that an individual i passes through a trap station x j on some occasion k
ϕ	Probability that an individual i is detected by a detector on some occasion k given that it is present at that trap

Bold symbols represent collections (vectors).

Table 2

Notations of latent variables and data used in this article

Latent variables	Definition
S	Locations of the activity centers of N animals within V
s i = (si 1, s i2)′	Location of individual i's activity center.
Z = (z 1, z 2, …, zM)′	A vector of Bernoulli variables, z i = 1 if individual i is present
u = (u1, …, uM)′	A vector of Bernoulli variables, ui = 1 if individual i is male in the population and ui = 0 if it is a female
u_0(⊂u)	Vector of “missing” binary observations on sexes of the list of M individuals
L=(L 1, L 2, …, LM)′	Each Li takes value in {1,2,…,M} and denotes the true index of ith detector 2 individual

Bold symbols represent collections (vectors).

Maximum number of individuals within the state space

Suppose a capture–recapture survey is conducted where a pair of detectors (1 and 2) are deployed at each of the 3 ( = J) trap stations and kept active for 4 ( = K) sampling occasions. Two individuals get fully identified based on their obtained capture histories (captured in both cameras at least once during the survey). The capture history for each of these two fully identified individuals is of dimension 2 × 3 × 4. The detection histories are tabulated in Table 3. Here, individual 1 is fully identified owing to the capture event at trap 2 on occasion 4. Individual 1 is also fully identified as it is captured at trap 2 on occasion 4. We note that in practical applications, full identity of individuals is obtained only when there are metadata, such as time of capture, to ascertain simultaneous captures. Due to the absence of simultaneous capture events in the detection histories of the partially identified individuals, we are uncertain about whether these histories correspond to two different individuals or to the same individual.

Table 3

An example of detection histories for two fully identified individuals and partially identified individuals is presented. The circled 1s indicate the simultaneous captures of an individual by the detectors 1 and 2

	Occasion trap	Detectors 1				Occasion trap	Detectors 2
		1	2	3	4		1	2	3	4
Fully identified individual 1	1	0	1	0	1	1	0	0	1	0
	2	1	0	0	①	2	0	0	0	①
	3	0	0	1	1	3	1	0	0	0
Fully identified individual 2	1	①	0	0	0	1	①	1	0	0
	2	0	0	0	1	2	0	0	0	0
	3	1	1	0	0	3	0	0	1	0
Partially identified individual	1	1	0	0	1	1	—	—	—	—
	2	0	0	1	0	2	—	—	—	—
	3	0	0	0	0	3	—	—	—	—
Partially identified individual	1	—	—	—	—	1	0	0	1	0
	2	—	—	—	—	2	1	0	0	0
	3	—	—	—	—	3	0	0	1	0

Model likelihoods

We have considered four models, denoted by M 1, M 2, M 3, and M 4, and the model likelihoods of the corresponding models follow in Equations ((1), (2), (3), (4)). Each of the four models M 1–M 4 are SECR models for partially identified capture–recapture history. The basic difference between M 3 and M 4 is that M 4 does not explain the latent hierarchy of the two events: trap entry of an animal and subsequent detection conditional on trap entry, whereas the former model M 3 does. M 1 uses the covariate information on sex category and allows the movement scale parameter (σ in our case) of M 3 to be gender specific. The differences between M 2 and M 4 are exactly the same (see Table 4). We describe the specific models in below.

Table 4

Specification differences in the four competing models

Model	Trap entry and detection parameter separated?	Sex‐specific σ (with sex covariate u)?
M 1	Yes	Yes
M 2	No	Yes
M 3	Yes	No
M 4	No	No

When detection rates in recorded samples are low due to failure or malfunction of detectors, capture–recapture data may be comprised of individuals with uncertain identities or “partially identified individuals.” Dey et al. (2019) separately accounts for the process of animal arrival within the detection region of a detector and detection process by conditioning on animal arrival—thus modeling the underlying mechanism by which we obtain different events leading to partial identification.

The probability of animal arrival η (s ) (termed as “trap entry probability”) is modeled as a decreasing function of Euclidean distance between individual activity center s and trap station x : . Here, ω 0 is regarded as the “baseline trap entry probability” and σ quantifies the rate of decline in trap entry probability as d(s , x ) increases. The observation process is parameterized in terms of detection probability ϕ which denotes the probability that any arbitrary individual i is detected by a detector on some occasion k given its arrival at that trap.

The obtained capture history observations and from the two detectors 1 and 2 during a spatial capture–recapture survey may not be synchronized as detectors often perform imperfectly. These two observed data arrays are then augmented with “all‐zero” capture histories. We denote the zero‐augmented data sets by Y ( ) and Y ( ); each of them is of dimension M × J × K, M being an upper bound of the population size. This also makes the dimension of the likelihood fixed in each iteration of the Markov Chain Monte Carlo algorithm which in turn eases computation. A vector of M latent binary variables z = (z 1,…,z)′ is introduced where z = 1 implies that individual i is a member of the population. We assume that each z is a Bernoulli random variable with parameter ψ and is independent of other z's. Here ψ is the proportion of individuals that are real and present within . Thus, the true population size N follows the Binomial distribution with parameters M and ψ. The individuals from the two lists obtained from detector 1 and detector 2, respectively, are linked probabilistically by introducing a latent identity variable L = (L 1, L 2,…, L)′ which is a one‐to‐one mapping from an index set of individuals captured by detector 2 to {1, 2,…,M} giving the true index of each of detector 2 individuals. Without loss of generality, the true identity of each individual in the population is defined to be in the row order of the capture histories of detector 1. Then, the rows of detector 2 data set Y are reordered as indicated by L to synchronize with the individuals of the detector 1 data set Y. We denote this newly ordered detector 2 data set as Y.

It is sometimes helpful to introduce a binary covariate on sex category u on spatial animal movement, σ, as in Sollmann et al. (2011). We define σ as a function of the latent structural vector : , if u= 1, that is, individual i is a male; σ(u) = σ, if u i = 0; that is, individual i is a female. u i's are independently and identically distributed Bernoulli random variables with parameter θ, θ being the probability that an arbitrary individual in the population is male. Let be a vector of binary observations on sex category of the captured individuals. The vector of latent missing observations in u is denoted by u 0. Assuming that covariate information on individual sex category is available, the joint density of Y * and u under Dey et al. (2019) is the following:

where denotes the probability that an individual i passes through a trap station x on some occasion k, is the number of times individual i got detected on at least one the detectors over K occasions and . The above Equation (1) can be regarded as zero‐inflated Bernoulli density with extra zeros coming from no trap entry. Prior to Dey et al. (2019), Royle (2015) proposed an SCR model for partially identified individuals coming from spatial capture–recapture surveys. The joint density of Y * and u under Royle (2015) is the following:where p (s , u i) = p 0 exp(−d(s i, x )2/(2σ(u)2)) denotes the probability that an individual i is detected at x on occasion k. Note that, unlike model (Equation 1), here movement through the detection region is considered inherently as a part of the observation process and p 0 is regarded as “baseline detection probability” and σ, although related to animal movement, is regarded as the rate of decline in detection probability. Qualitatively, the absence of ϕ in (Equation 2) distinguishes this model from (Equation 1) and can be regarded as a less general model. Recently, Augustine et al. (2018) extended Royle (2015) by introducing separate parameters to distinguish captures by both detectors and captures by only one of the detectors.

In the absence of the sex covariate u, the joint density of under Dey et al. (2019) and Royle (2015) is, respectively, as given below: where denotes the probability that an individual i passes through a trap station x on some occasion k and p(s ) = p 0 exp(−d(s i, x )2/(2σ 2)) denotes the probability that an individual i is detected at a trap station x on some occasion k.

The assumed prior distributions for the model parameters are as follows: a uniform distribution over the interval (0,1) for the probability parameters ϕ, ω 0, p 0, ψ, and θ; a uniform distribution over the interval (0,R) for parameters σ, σ, and σ where R is high enough to expect that it would be impossible for animals to exhibit movement as widely as this scale during sampling. The prior distributions for these model parameters ϕ, ω 0, p 0, ψ, and θ, σ, σ are assumed to be independent. To ensure that the marginal distribution of the data is well defined, we have assumed proper priors for each of these parameters (Gopalaswamy & Delampady, 2016). We assume a uniform prior over the entire state space for each location of the activity center s and that these s 's are independently distributed. The latent variable L is assumed to have a uniform prior distribution over the permutation space of {1, …, M}. The prior specifications remain the same for all the model fits. The MCMC algorithm used to sample from the respective posterior density under each model is detailed in Appendix S1D. Thus, we have the four models,

M 1: Model with density (Equation 1),

M 2: Model with density (Equation 2),

M 3: Model with density (Equation 3), and

M 4: Model with density (Equation 4).

Model selection tools are used to find the “best” fitted model among a set of competing models. To provide an understanding of what such a best model is, Shmueli (2010) discusses two broad modelling philosophies—explanatory and predictive. Explanatory modeling of the data set is only relevant to obtain the most accurate representation of the underlying theory, whereas predictive modelling seeks to minimize the combination of bias and estimation variance. Hence, a model that explains the data well may not have the best predictive ability and vice versa. This naturally leads to two performance aspects by which a model can be assessed—(a) how well does the model explain the observed data set, that is, to what degree of confidence can it be thought of as a model from which the data generated, (b) how good is its predictive ability.

On such criteria, a wide range of model selection tools have thus emerged in the applied statistical literature (see Gelman, Hwang, & Vehtari, 2014; Höge, Guthke, & Nowak, 2019; Höge, Wöhling, & Nowak, 2018; Hooten & Hobbs, 2015; Kass & Raftery, 1995; for examples). In all these cases, however, we are assuming that the “true” model is in the list of competing models. Such classes of problems are formally categorized as ‘‐closed’ (Bernardo & Smith, 1994; Clarke, Clarke, & Yu, 2013; Vehtari & Ojanen, 2012) and is generally considered to be well‐studied (Clarke, Clarke, & Yu, 2013).

However, for the practitioner, the applicability of such guidelines for model selection is not clear because often such guidelines are formed using specific examples (Hooten & Hobbs, 2015) or based on asymptotic arguments in a theoretic sense (for example, Shibata, 1989; Stone, 1977; Watanabe, 2010). Even when the purpose of model selection is clear (explanatory or predictive) and also the model selection criterion is fixed, what techniques practically work is not clear. Practitioners often wish to confront fairly complex models (e.g., hierarchical models) with their data, and it is often a challenge to directly apply a candidate set of model selection tools without incorporating some approximations. Given these approximations due to model complexity and with limited sample sizes, the validity of such general guidelines (Hooten & Hobbs, 2015) remains to be tested in practice. Therefore, in this study, we apply a range of model selection tools and test the validity of such general model selection guidelines in the context of a popular class of models, spatial capture—recapture models, used for estimating animal density of some of the world's most iconic species (Broekhuis & Gopalaswamy, 2016; Elliot & Gopalaswamy, 2017; Royle, Karanth, Gopalaswamy, & Kumar, 2009; Sollmann et al., 2011).

Candidate model selection tools

We have considered four different Bayesian model selection methods for application and evaluation: Bayes factors, Deviance Information Criterion (DIC), Watanabe‐Akaike information criterion (WAIC), and posterior predictive loss. Two popular model selection tools (AIC and BIC; Burnham & Anderson, 2002) are not used here because they impose restrictive assumptions on the parameter space as the sample size increases—situations often encountered in many hierarchical models (Royle & Dorazio, 2008). For example, in the SCR models we study here, the concept of “number of parameters” is unclear and we therefore cannot apply criteria such as the AIC and the BIC directly.

We clarify our aim here is to discriminate between hypotheses (models) and not to conduct hypothesis tests of the model parameters. In the hypothesis testing paradigm that applies in our setup, it usually requires for models to be nested so that statistical tests (such as likelihood ratio tests) between different parameter values can be performed. However, this limitation does not apply to model selection under the hypothesis discrimination paradigm.

Bayes factors

Model comparison using Bayes factors requires the computation of the marginal likelihood, which involves the integration where denotes the model likelihood and π(µ) denotes the prior density of the parameters µ under M. This integration is difficult to compute in practice unless the models are very simple in structure, which is often not the case in ecology. Therefore, computation of marginal likelihood continues to be an active area of research in Statistics (Wang & Meng, 2016; Wang et al., 2018).

Estimation of marginal likelihood of data

Under our model settings, Y = (Y *, u) for models M 1, M 2, and Y = Y * for the other models M 3, M 4. µ denotes the collection of all parameters and latent variables for each model as a generic notation. Specifically, let µ = (µ , µ s), where µ is the collection of scalar parameters and µ is the collection of all latent variables. The Gelfand‐Dey estimator of marginal likelihood of data m(Y) is expressed as:where {µ ( ): d = 1, …, N iter} is a set of MCMC draws from the posterior π(µ|Y) and g(µ) is a tuning density. It is to be noted that by specifying g(µ) = π(µ) in (Equation 5), we obtain the harmonic mean estimator of the marginal likelihood

Details on these estimators and their properties can be found in Gelfand and Dey (1994), Kass and Raftery (1995). The harmonic mean estimator is easy to compute by just calculating the model likelihood at each MCMC sample draws from the posterior distribution. Although it is known that the harmonic mean estimator is a consistent estimator of the marginal likelihood, as noted in Kass and Raftery (1995) and Newton and Raftery (1994), it can have a simulation pseudo‐bias (Lenk, 2009). This bias occurs due to the limitations in numerical computations. Indeed, our study here aims to assess how and whether these potential biases in computation influence the model selection process for the particular class of models we consider—partial identity SCR models, in our case.

For our problem, the computation of (Equation 5) requires us to obtain the integrated likelihoods (marginals) under the different models that we consider. This becomes particularly tricky in the presence of high‐dimensional latent variables such as u 0, z, S, and L, which are elements of µ .

We have developed a novel approach to compute an estimate of the marginal likelihood. Here we partition the parameters of a model into two sets (real‐valued scalar and high‐dimensional structural parameters, which are denoted by µ and µ , respectively) and then propose two approximating strategies to account for the high‐dimensional structural parameters, , before applying the Gelfand‐Dey estimator. This general approach of partitioning can be applied to many other classes of hierarchical models as well and not limited to SCR models.

We describe below two approximating approaches to compute the Gelfand‐Dey estimator: the maximum a posteriori (MAP) approximation approach and the integrated likelihood (IL) approach.

Approach 1: MAP approximation

In this approach, we fix the high‐dimensional variables at their MAP estimates , assuming that their posterior distributions are well summarized by these estimates, which are derived from the MCMC draws. The Gelfand‐Dey estimator is then computed using the formula,where is a set of MCMC draws from the posterior π(µ , µ |Y). We begin with as an initial estimate of (µ, µ) where

This estimate of the posterior mode of (µ, µ) may not be optimal since in our high‐dimensional parameter setting, an MCMC sample of a practical size may not be enough to extensively explore the posterior surface. We, therefore, fix one of the parameters and explore the posterior surface to find d 1 such that . In this way, we obtain an improved MAP estimate of (µ , µ ), , if

Similarly, we then fix and find . This procedure is continued iteratively to eventually give us the best MAP estimate of the posterior mode . Suitable transformations of the parameters ensure that all the points in belong to the posterior support. A single MCMC chain is to be used for the algorithm, and we have shown that the obtained MAP estimate using our algorithm will be better than the usual estimate (µ , µ ) (see Appendix S1A). In theory, we can also increase sample size by merging all the chain outputs after convergence.

Approach 2: Integrated likelihood (IL) approximation

Ideally, we would like to compute the marginal likelihood m(Y) by integrating out all the latent variables with respect to their corresponding prior distributions from the model likelihoods. However, in the case of model likelihoods (Equations (1), (2), (3), (4)), this integration is computationally difficult for the permutation vector L because of its very large support. The integration over the variables u 0 and z can be performed analytically. The integration over S is evaluated numerically by partitioning the region into a sufficiently fine grid and then evaluating a Riemann sum (as direct integration cannot be expressed in a closed form). This integrated likelihood can then be used in (Equation 5) for estimating the marginal likelihood

One downside of the IL approximation approach is the lack of clarity about the interdependencies between the latent variables u 0, z, S, and L after carrying out the integrations. The derivations of the integrated likelihoods for each of the four models M 1–M 4 are given in Appendix S1B by ignoring any possible interdependencies. We have also assessed the robustness of the Gelfand‐Dey estimator by computing the marginal likelihood estimates using different tuning densities, for example, multivariate normal density, multivariate‐t density with varying degrees of freedom, and the truncated normal density following the suggestion of Geweke (1999). These technical details are described in Appendix S1C.

Deviance information criterion

Deviance is defined as D(µ) = −2log f(Y|µ). The deviance information criterion (DIC) is then defined as , where is an estimate of µ. The term p DIC can be viewed as effective number of parameters, which is a bias correction in evaluating a model's predictive accuracy. Two versions of p DIC are generally used (Gelman, Carlin, et al., 2014; Spiegelhalter, Best, Carlin, & van der Linde, 2002):

Although p DIC1 has been considered to be numerically more stable, p DIC2 has the advantage of being always positive (Gelman, Hwang, & Vehtari, 2014). A model with smaller DIC value is preferred. Model comparison using DIC is not invariant to parameterization and depends on the components of the model likelihood to be considered as the likelihood. Spiegelhalter et al. (2002) suggests that practitioners carefully decide on the parameters of interest so that they can avoid this potential pitfall. This advice is often difficult to implement in practice, especially when there exists inherent ambiguity in the interpretation of latent parameters. We note that Celeux, Forbes, Robert, and Titterington (2006) suggests several forms for the DIC that can be used for different hierarchical models but does not recommend any particular form as the best.

For the computation of deviance, we have use the MAP estimate of µ to obtain instead of the posterior mean, due to the presence of binary latent variables and unknown permutation vectors in the likelihood. We have then computed two versions of p DIC (Gelman, Carlin, et al., 2014; Hooten & Hobbs, 2015) using the MCMC draws from π(µ|Y) as follows:

Watanabe‐Akaike information criterion

The Watanabe‐Akaike information criterion (WAIC) is a Bayesian version of the AIC as it uses the posterior predictive distribution of the data to estimate the out‐of‐sample predictive accuracy of the model. Watanabe (2010) introduced the WAIC criterion based on the assumption of independence between the data points and has shown its asymptotic equivalence with cross‐validation. In our model formulations, we have assumed that the different data points correspond to capture–recapture data set for each of the M individuals. The WAIC is then defined as . A model with smaller WAIC value is preferred. In computing WAIC, we partition data Y in terms of individuals (Y 1, Y 2,…, Y). We compute the two commonly used versions of p WAIC (Hooten & Hobbs, 2015) using MCMC draws{µ (d): d = 1,…,N iter} from π (µ|Y) as follows:

We propose another version for p WAIC based on absolute error loss, recognizing that large variability in µ ( ) and the magnitude of squared errors themselves may have an impact on the efficiency of p WAIC of the above forms:

From a decision‐theoretic perspective, we think it is interesting to compare WAIC2, which is based on a square error loss function with WAIC3 that is based on an absolute error loss function. The comparison based on a simulated data analysis is described later in Section 4.

Posterior predictive loss

Gelfand and Ghosh (1998) derived a model selection criterion, popularly known as the posterior predictive loss criterion, by adopting a decision‐theoretic approach for measuring predictive accuracy of a model. The posterior predictive loss D criterion (based on a square error loss function) is composed of squared bias, variance of a new prediction, and a parameter w. The parameter w indicates the relative weight given to the loss for departure of the new predicted value from the observed data against departure of the same from the new data. In practice, this weight is taken to be very large, and asymptotically as w →∞, the D ∞ criterion is obtained as follows under our model setting:where is a vector of length 2MJK obtained by vectorizing observed data array Y and Y rep = (y1,rep, y2,rep,…,y2MJK,rep)′ represents the vector of all these replicates. The first term in the D ∞ criterion (see Equation 13) is the goodness‐of‐fit term, while the second term can be interpreted as a penalty term for model complexity. The model with the smallest D ∞ is to be preferred. In our analysis, the two data sets are obtained by vectorizing the two binary data arrays.

We compute the above expectation and variance using the MCMC draws. Given an MCMC sample µ (d): d = 1,…,N iter from π(µ|Y), we simulate from f(Y|µ ( )) for each d = 1,…,N iter. For instance, in model M 1, µ denotes the collection of the parameters ψ, θ, ϕ, ω 0, σ, σf, u 0, z, S, L. Then and . We summarize the various model selection methods and their variants in Table 5. Considering all these model selection tools and their variants (from approximation approaches to setting tuning densities), our evaluation is carried out on 25 unique tools. By no means is this an exhaustive list of candidate tools because there are a large number of tools available. Our purpose here really is to assess some popular tools and assess their efficacies on a particular class of models.

Table 5

Bayesian model selection methods used in this study

Sl. no.	Model selection method	Variant	Approximation method	Choices of tuning density	Eq. No.
1.	Bayes factor	Gelfand‐Dey estimator	MAP	Multivariate normal density, multivar‐ iate‐t density with degrees of freedom 10, 100, 500, 1,000, 10,000 and truncated multivariate normal density with confidence coefficients 0.90, 0,95, 0.99	(7)
2.	Bayes factor	Gelfand‐Dey estimator	IL	‐Do‐	(8)
3.	Bayes factor	Harmonic mean estimator	—	—	(6)
4.	DIC	p DIC1	MAP	—	(10)
5.	DIC	p DIC2	MAP	—	(10)
6.	WAIC	p WAIC1	—	—	(11)
7.	WAIC	p WAIC2	—	—	(11)
8.	WAIC	p WAIC3	—	—	(12)
9.	Posterior predictive loss	—	—	—	(13)

EVALUATION OF THE PERFORMANCE OF MODEL SELECTION METHODS

Simulation design and simulation scenarios

We have conducted simulations for 12 scenarios (provided in Table 6) grouped into 2 equal sized sets, to assess the performance of the models proposed here. We set σ = 0.3 and σ = 0.15 for the first set of 6 scenarios, σ = 0.4 and σ = 0.2 for the second set of 6 scenarios. We set (ω 0, ϕ) = {(0.01, 0.3), (0.05, 0.3), (0.05, 0.5), (0.03, 0.8), (0.01, 0.9), (0.05, 0.9)}, which gives us 6 different scenarios for each of the two sets corresponding to the values taken by ω 0 and ϕ. We assume that a total of 100 individuals are residing inside the state space of which 40 are male. Each of the simulation experiments is conducted within a rectangular state space of dimension 5 unit × 7 unit (Figure 1), after setting a buffer of 1 unit in both horizontal and vertical directions, a 10 × 16 trapping array of total J = 160 trap stations has been set (trap spacing is 0.3 unit on X axis and 0.3125 unit on Y axis). This meets the requirement suggested in Karanth and Nichols (2017). Each of the traps remains active for K = 50 sampling occasions simultaneously. For parameter estimation, we set the maximum possible number of individuals present in the population (M) at 400 for all the scenarios. Since MCMC approaches are computationally time consuming, we set the data‐generating parameters to generate data sets of varying levels of information content. This strategy ensured that we capture trends in model selection performance by the different methods. The experiment is repeated n sim = 10 times. The MCMC samples for each of the parameters are obtained (each of length 30,000), and the estimates are computed using those chains after a burn‐in of 10,000. It takes approximately 3 days to run a single chain of 30,000 iterations using our R code for model M 1, approximately 2 days for model M 2 and M 3, and approximately 1 day for model M 4 using our R code on a DELL Precision Rack 7,910 Server with 40 cores and 512 GB RAM at a clock speed of 2.20 GHz.

Table 6

Parameter specifications corresponding to different simulation scenarios

Scenario	M	N	N Male	ω 0	ϕ	σm	σf
1	400	100	40	0.01	0.3	0.3	0.15
2	400	100	40	0.01	0.9	0.3	0.15
3	400	100	40	0.01	0.3	0.4	0.20
4	400	100	40	0.01	0.9	0.4	0.20
5	400	100	40	0.03	0.8	0.3	0.15
6	400	100	40	0.03	0.8	0.4	0.20
7	400	100	40	0.05	0.3	0.3	0.15
8	400	100	40	0.05	0.5	0.3	0.15
9	400	100	40	0.05	0.9	0.3	0.15
10	400	100	40	0.05	0.3	0.4	0.20
11	400	100	40	0.05	0.5	0.4	0.20
12	400	100	40	0.05	0.9	0.4	0.20

Figure 1

Array of trap locations (denoted by “+”) within the state space (0.5) × (0.7)

Capture–recapture data sets are simulated independently under each of the 12 simulation scenarios (Table 6) under model M 1. Recall that, model M 1 corresponds to the statistical model in (Equation 1) with σ parameter modeled in terms of individual covariate on sex category u (see Section 2.1). Then, each simulated data set is fitted with all the four competing model M 1, M 2, M 3, and M 4. For a better mixing during the MCMC run, we transform the components of so that the transformed parameter space is the entire Euclidean space (details are given in Appendix S1C).

Defining performance measures

Probability of selecting the true model

In our study, since all the data sets are simulated from model M 1, it is considered as the true model. We have computed the proportion of times a model selection method chooses M 1 as the best out of n = 10 simulations. This proportion will serve as an estimate for the probability of selecting the true model. Additionally, we have also computed the proportion of times a model selection method chooses M 2, M 3, and M 4, respectively. The computed proportions of selecting these models from the simulation study will indicate the efficacy of the model selection methods.

Assessing the quality of parameter estimation

The abundance parameter N carries a lot of significance in ecology and conservation. Due to its importance, ecologists place their interests in the robustness and accuracy of its estimate, and will therefore rely on a model selection method that will achieve this by selecting that model from the model set that will provide the most reliable estimate of N. The precision and accuracy of the parameter estimate indicates the quality of the model fit, and we assess this by computing the average root mean square error (average RMSE).

Suppose {µ ( ): d = 1, …, N iter} denotes a set of MCMC draws from the posterior distribution of an arbitrary parameter µ for the t‐th simulated data set, t = 1,…,n sim. Mean square error (MSE) of µ for the t‐th simulated data set is estimated as Average RMSE is calculated by first averaging the estimated MSEs of different simulations and then taking the square root of the average:

Let us now summarize our experiment. The quantities related to the various Bayesian model selection methods (Section 2.2) are computed for each of the four models for the simulated data and analysis sets. We compare the conclusions drawn from the computed proportions with the conclusions drawn from average RMSEs of the parameters to study the behavior of the model selection methods with respect to varying information content. We also generate pairwise correlation plots from the MCMC draws to study the extent of parameter redundancy between various pairs of parameters as a consequence of lack of information content in the data.

RESULTS OF THE SIMULATION STUDY

Figure 2 shows the proportion of times different estimators of Bayes factor favors any particular model. The subplots have different scenarios 1–12 in the x‐axis and the competing models M 1–M 4 in the y‐axis. Plot (a) and Plot (b) correspond to the Gelfand‐Dey estimator of the Bayes factor favoring any particular model using (i) the MAP approximation approach and a multivariate normal density for tuning density g, (ii) the integrated likelihood approximation approach and a multivariate normal density for tuning density g, respectively. Plot (c) corresponds to the harmonic mean estimator of the Bayes factor. Plots (d–i) correspond to WAIC1, WAIC2, WAIC3, DIC1, DIC2, and posterior predictive loss, respectively.

Figure 2

Plot (a): The proportion of times Gelfand‐Dey estimator of Bayes factor favors any particular model using the MAP approximation approach and a multivariate normal density for g. Plot (b): The proportion of times Gelfand‐Dey estimator of Bayes factor favors any particular model using the integrated likelihood approximation approach and a multivariate normal density for g. Plot (c): The proportion of times harmonic mean estimator of Bayes factor using the favors any particular model. Plots (d)–(i) correspond to WAIC1, WAIC2, WAIC3, DIC1, DIC2, and posterior predictive loss, respectively

Our results suggest that the choice of the tuning density in computation of Gelfand‐Dey estimator had no impact on model selection (Figure 2). Thus, we focus only on the results corresponding to the normal tuning density. We observe that the Bayes factor with the GD‐MAP approximation favors M 4 in 70% of the times under scenario 1 and favors M 3 in at least 70% of the times under scenarios 2, 3, 4, and 7. Bayes factor (GD‐MAP) is in favor of M 1 more than 80% of the times under all the remaining scenarios 5–6 and 8–12. As seen here, the GD‐IL approximation performs slightly worse than the GD‐MAP approximation (Figure 2). However, surprisingly, we observe that the model choices by the harmonic mean estimator of Bayes factor performs well and favors the true model M 1 in majority of the scenarios.

We have considered three forms of the WAIC and two forms of the DIC. The corresponding plots are shown in Figure 2(d–h). WAIC1, WAIC3, DIC1, and DIC2 exhibit very similar tendencies in their model choices by favoring M 3 under all the scenarios except scenarios 1 and 3. These four methods favor M 4 under scenario 1. Under scenario 3, WAIC1 and DIC1 favor M 3 in majority of the times, whereas WAIC3 and DIC2 favor M 4. Like Bayes factor, these model selection methods (DICs and WAICs) also tend to favor simpler models under scenarios 1 and 3. The WAIC2 largely agrees with the other WAICs (and DICs) but more often selects the true model when data sets are more informative (Figure 2d–f). In general, DICs and WAICs seem to discourage the presence of high‐dimensional latent variables.

The plot showing the proportion of different model choices by posterior predictive loss is given in Figure 2i. The posterior predictive loss criterion D ∞ favors models with individual sex covariates (M 1 and M 2). The posterior predictive loss also appears to select the true model some of the times, even when there is not sufficient information in the data (scenarios 1 and 3).

In Figure 3, we see that the average RMSE of N is substantially higher in scenarios 1–3 relative to scenarios 4–12. We also have observed a substantial amount of correlation between some pairs of parameters involving N in the MCMC chains under scenarios 1–3 (r(N, θ) ≈ −0.7, r(N, σ ≈ −0.42). However, these correlations decrease for the other scenarios, likely due to increased information content in the data. For interested readers, we have also provided additional plots for the model selection methods (Figures 1–6), the RMSE plots for all the parameters (Figures 7 and 8), and a representative set of scatter plots (Figures 9–14) in the Appendix S1E.

Figure 3

Plot of average RMSE estimates of N over different simulation scenarios

DISCUSSION AND CONCLUSIONS

Contemporary practice of ecology and conservation biology relies largely on the use of model selection for hypotheses discrimination (Johnson & Omland, 2004). Simultaneously, there has been major growth in the use of hierarchical models in ecology, especially within the realm of Bayesian inference (Kéry & Royle, 2015; Royle & Dorazio, 2008). These models have now enabled statistical ecologists to fairly easily formulate complex ecological models, elegantly deal with the sampling process and also fit these complex models using powerful tools such as the MCMC (Kéry & Royle, 2015). However, the lack of availability of ready‐made model selection tools when practicing the Bayesian inference has sometimes motivated ecologists to continue using the likelihood‐based inferences, merely because one can use well‐known model selection tools such as the AIC (Burnham & Anderson, 2002) for inference.

To provide a context for this argument, in the spatial capture–recapture literature, we have essentially seen the development of three important likelihood functions: (1) Borchers and Efford (2008), (2) Royle et al. (2009) and (3) Royle, Sutherland, Fuller, and Sun (2015). Inferences for the models (1) and (3) are by maximizing the likelihood, while the inference for (2) is Bayesian. We note with interest that one of the reasons for the development of (3) was motivated on the pretext that model selection is much easier (using known tools such as the AIC) for practitioners using the likelihood approach, in spite of the problem having been already solved in the Bayesian context (Royle et al., 2013). It is specially of concern in the context of the models we study here in that investigators may be forced to integrate out S (activity centers of individuals) in order to construct tractable likelihoods and thus oversimplifying ecological reality. For example, by retaining the activity centers it is possible for investigators to confront open SCR models (Gardner, Sollmann, Kumar, Jathanna, & Karanth, 2018). We therefore hope that our study motivates the continued use of Bayesian methods by investigators for all its advantages, instead of opting out of them merely for the sake of using simple model selection methods.

In this study, we have tried to implement some selected Bayesian model selection methods on a specific class of advanced Bayesian SCR models (Royle, 2015—with and without sex covariates; Dey et al. (2019)–with and without sex covariates) dealing with partially identified individuals. We have found our Bayes factor implementation using the Gelfand‐Dey estimator (using the MAP approximation approach) to be the preferred choice as a model selection method over a wide range of simulation scenarios. This approach appears to work particularly well when information content in the data is moderate to high. The IL approximation approach also worked well but not as well as the MAP approximation approach perhaps because there exists interdependency between some of the latent variables which are not accurately captured during marginalization.

However, implementing Bayes factors for model selection using the Gelfand‐Dey estimator (with MAP approximation) is not very straightforward and requires investigators to incorporate some code during model fitting by MCMC. Interestingly, our study demonstrates that obtaining Bayes factors using the harmonic mean approach for marginal likelihood computation is less demanding but yet serves as a very good model selection method. This finding deviates from the popular view among applied scientists that it is futile to estimate the marginal likelihood using the harmonic mean approach (Lartillot & Philippe, 2006; Xie, Lewis, Fan, Kuo, & Chen, 2011). We surmise that this finding may be attributed to the fact that when the priors are bounded (as we have done, but for other reasons) and does not permit extremely low probabilities to occur at the tails. Further, our emphasis here is to ask how the tool performs to select the true model and estimate parameters for prediction for a given analytical setup. And it is not about how accurately the tools estimate the actual value of the marginal likelihood. Thus, many of the criticisms (Lartillot & Philippe, 2006; Xie et al., 2011) may become irrelevant in practice, but this would require further enquiry.

As our simulation study shows, the two goals of model selection and parameter estimation cannot be simultaneously achieved under certain circumstances (scenarios 1–3), especially when information content is low (indicated by high RMSE values and correlation coefficients). Hence, researchers have to clearly prioritize their objectives prior to data analysis. If the goal, for example, is to find a model that best estimates population size N, then we recommend the use of Bayes factor (Gelfand‐Dey with MAP approximation) or the Bayes factor (harmonic mean, due to its simplicity) because these appear to provide the most reliable estimates of N over all the simulation scenarios. However, if researchers are only interested to select the true model, especially when data are less informative (scenarios 1–3), we recommend the posterior predictive loss approach since they favor the true model nearly a one‐third of the times only in situations with such low information content. Of course, the dual objectives of model selection and parameter estimation are met when information content is moderate or high (scenarios 4–12), and as stated previously, we recommend either the Bayes factor (Gelfand‐Dey with MAP approximation) or the Bayes factor (harmonic mean) in such cases. However, the posterior predictive loss (with the squared error loss function) as used here does not select the true model when information content is moderate to high.

We also do not recommend the use of the DICs or the WAICs, since they do not appear to outcompete other model selection tools, either from the standpoint of model selection or parameter estimation, in any of the simulation scenarios. This is an interesting finding, because tools such as the WAIC are much newer tools developed by applied Bayesians and recommended for hierarchical models (Gelman, Carlin, et al., 2014; Hooten & Hobbs, 2015). Thus, our study brings back focus on the need to assess the strength of inference from a model selection method by conditioning on a true model and consequently evaluating a competing set of model selection methods prior to selecting the most appropriate one for the problem on hand. Table 7 summarizes the similarities and differences in inferences between our specific evaluation of model selection tools and the generally accepted evaluations. We believe this summary will motivate investigators to look deeper into their own models and the model selection methods they are using.

Table 7

Contrasting the performance of model selection tools based on the intended purpose and perceived applicability with findings from our specific study. In the table, we provide answers to the following questions: (a) Does this approach select the true model? (b) Does this approach favor models providing reliable estimates of parameters (specifically for N)? (c) How difficult is the approach to implement in practice? Comments in bold draw attention to the noticeable differences between the expected performance of a tool and its performance in our particular study

Model selection tool	Intended purpose and applicability	Reference	Findings from our specific study
Bayes factor (by Gelfand‐Dey estimator)	(a) Yes (b) Yes (c) Difficult		(a) Yes, very often (MAP) Yes, quite often (IL) (b) Yes, very often (MAP) Yes, quite often (IL) (c) Moderately difficult (MAP) Moderately difficult (IL)
Bayes factor (by harmonic mean estimator)	(a) Yes, but unreliable (b) Yes (c) Easy	Newton and Raftery (1994), and Kass and Raftery (1995)	(a) Yes, very often (b) Yes, very often (c) Easy
DIC1	(a) No (b) Yes (c) Easy	Spiegelhalter et al. (2002), Gelman, Hwang, & Vehtari, 2014), and Hooten and Hobbs (2015)	(a) No (b) No (c) Moderately difficult
DIC2	(a) No (b) Yes (c) Easy	Spiegelhalter et al. (2002), Gelman, Hwang, & Vehtari, 2014), and Hooten and Hobbs (2015)	(a) No (b) No (c) Moderately difficult (required MAP estimate)
WAIC1	(a) No (b) Yes (c) Easy	Watanabe (2010), Gelman, Hwang, & Vehtari, 2014), and Hooten and Hobbs (2015)	(a) No (b) No (c) Easy
WAIC2	(a) No (b) Yes (c) Easy	Watanabe (2010), Gelman, Hwang, & Vehtari, 2014), Hooten and Hobbs (2015)	(a) No, not often (b) Yes, but suboptimally (c) Easy
WAIC3	(a) No (b) Yes (c) Easy		(a) No (b) No (c) Easy
Posterior predictive loss	(a) No (b) Yes (c) Moderately difficult	Gelfand and Ghosh (1998), and Hooten and Hobbs (2015)	(a) No, not often (b) Yes, but suboptimally (c) Moderately difficult

Our approach does not, strictly speaking, permit us to draw conclusions and make inferences on the most suitable model selection tools beyond the restrictive set of competing models and the settings we have used in this study. However, we demonstrate the applicability of the Bayes factor on a set of structurally very complex hierarchical models. This was possible by using our approach of partitioning the real‐valued scalar and high‐dimensional structural parameters. We anticipate that this general approach could address difficulties faced by ecologists attempting to use the Bayes factors for comparing very complex hierarchical models. We would expect that it would be much easier to use the Bayes factor for relatively simpler and a much larger class of hierarchical models (Kéry & Royle, 2015; Royle et al., 2013). We note with interest that it is also unclear whether the routinely used AIC works as an appropriate model selection tool for MLE‐based SCR models as discussed in Efford and Mowat (2014). In real‐world problems, we can almost always find proper priors for parameters to ensure that Bayes factors exist. Therefore, we believe it is best to depend on tools built on sound statistical theory (by finding suitable approximations) rather than seeking answers based on seemingly elegant but not well‐tested tools.

CONFLICT OF INTEREST

None declared.

AUTHORS CONTRIBUTION

SD, MD, and AMG conceived the ideas and designed methodology; SD conducted data analysis; SD, MD, and AMG drafted the article. All authors contributed critically to the drafts and gave final approval for publication.

Click here for additional data file.