Estimating the Effects of Habitat and Biological Interactions in an Avian Community

We used repeated sightings of individual birds encountered in community-level surveys to investigate the relative roles of habitat and biological interactions in determining the distribution and abundance of each species. To analyze these data, we developed a multispecies N-mixture model that allowed estimation of both positive and negative correlations between abundances of different species while also estimating the effects of habitat and the effects of errors in detection of each species. Using a combination of single- and multispecies N-mixture modeling, we examined for each species whether our measures of habitat were sufficient to account for the variation in encounter histories of individual birds or whether other habitat variables or interactions with other species needed to be considered. In the community that we studied, habitat appeared to be more influential than biological interactions in determining the distribution and abundance of most avian species. Our results lend support to the hypothesis that abundances of forest specialists are negatively affected by forest fragmentation. Our results also suggest that many species were associated with particular types of vegetation as measured by structural attributes of the forests. The abundances of 6 of the 73 species observed in our study were strongly correlated. These species included large birds (American Crow (Corvus brachyrhynchos) and Red-winged Blackbird (Agelaius phoeniceus)) that forage on the ground in open habitats and small birds (Red-eyed Vireo (Vireo olivaceus), House Wren (Troglodytes aedon), Hooded Warbler (Setophaga citrina), and Prairie Warbler (Setophaga discolor)) that are associated with dense shrub cover. Species abundances were positively correlated within each size group and negatively correlated between groups. Except for the American Crow, which preys on eggs and nestlings of small song birds, none of the other 5 species is known to display direct interactions, so we suspect that the correlations may have been associated with species-specific responses to habitat components not adequately measured by our covariates.

S2 Appendix: MCMC algorithm used to fit multispecies N-mixture models In this supplement we describe the MCMC algorithm that was used to compute summaries of the posterior distribution of model parameters. We use bracket notation [1] to specify probability density functions; thus, [x, y] denotes the joint density of random variables X and Y , [x|y] denotes the conditional density of X given Y = y, and [x] denotes the unconditional (marginal) density of X.
We used the MCMC algorithm to generate a Markov chain whose stationary distribution is equivalent to a posterior with density function proportional to where β = (β 1 , . . . , β K ), α = (α 1 , . . . , α K ), ε = (ε 1 , . . . , ε I ), and n = (n 1 , . . . , n I ) are parameters of the multispecies N-mixture model. The vector b = (b 1 , . . . , b K ) contains auxiliary parameters that were used to specify the prior distribution of Σ hierarchically [2] as described below. The vector y ik = (y ik1 , . . . , y ikJ i ) contains the detection frequencies of birds of species k observed during the ith survey, and y = {y ik , ∀ik} denotes the entire set of observations. The MCMC algorithm uses either a Gibbs or a Metropolis-Hastings (M-H) sampler depending on the parameter. Each of the following full-conditional distributions is sampled in one iteration of the algorithm: 1. The full conditional for N ik has a familiar form: where y ik· = J i j=1 y ikj is the total number of birds of species k detected during the ith survey. In other words, this is the conditional distribution of the number of birds that were present but not detected during the survey.
2. We assumed a hierarchical prior distribution for Σ [2] that allowed us to specify marginally noninformative priors for its elements. Specifically, the hyperparameters of this prior can be chosen so that each standard deviation parameter σ k has a Halft prior of arbitrarily high noninformity [3] and each correlation parameter ρ kl has a uniform prior on (-1,1). The hierarchical prior distribution for Σ is a mixture of Inverse-Wishart and Inverse-Gamma distributions: . [2] showed that the marginal prior density for ρ kl is proportional to (1 − ρ kl ) ν/2−1 ; therefore, we let ν = 2 to specify a marginally uniform prior for each correlation parameter. [2] also showed that the marginal prior distribution for σ k is a Half-t distribution with ν degrees of freedom and scale parameter s k ; therefore, we specified a noninformative prior for σ k by choosing s k to be arbitrarily high.
The conditional conjugacy of this prior for Σ leads to full conditional distributions of familiar form that are relatively easy to sample. Specifically, where (Σ −1 ) kk is the kth diagonal element of Σ −1 . Once Σ −1 is drawn from its full conditional, it is a simple matter to compute the matrix of correlation parameters as 3. The full conditional density of ε factors into a product of I independent terms. To sample the ith target density, which is proportional to we used the M-H algorithm. Following [4], we used a multivariate-t distribution as a proposal and selected its parameters to approximate the target distribution.
Specifically, let f (ε i ) = log([ε i |·]). We assigned the mean of the proposal distribution to equalε i , the value of ε i that maximized f (ε i ). This maximimization was done numerically using an analytical gradient g(ε i ) = −Σ −1 ε i + n i − λ i and hessian 4. The full conditional density of β factors into a product of K independent terms. For each term we assumed a normal prior distribution with mean zero and diagonal covariance matrix V β k , using an arbitrarily high value for the nonzero elements of V β k to specify prior ignorance. To sample the kth target density, which is proportional to we used the M-H algorithm. Following [4], we used a multivariate-t distribution as a proposal and selected its parameters to approximate the target distribution. Specif- . We assigned the mean of the proposal distribution to equalβ k , the value of β k that maximized f (β k ). This maximimization was done numerically, and the covariance of the proposal distribution was then computed by inverting the negative of the hessian matrix {−H(β k )} −1 . The degrees of freedom parameter of the proposal may be used as a tuning parameter; however, in practice we simply assigned this parameter to be constant (2.0).

5.
The full conditional density of α factors into a product of K independent terms. Let Q k equal the number of parameters in α k . We assumed a noninformative prior for α k that comprises Q k independent t-distributions, each with mean zero, scale parameter σ α = 1.566, and degrees of freedom parameter ν α = 7.763. This distribution approximates a Uniform(0, 1) prior on the inverse-logit scale and assigns low probabilities to values outside the interval (-5,5) [5]. To sample the kth target density, which is proportional we used the M-H algorithm. We used a multivariate-t distribution as a proposal and selected its parameters to approximate the target distribution. Specifically, let f (α k ) = log([α k |·]). We assigned the mean of the proposal distribution to equalα k , the value of α k that maximized f (α k ). This maximimization was done numerically, and the