Analyses
We modeled detection-weighted, single-species occupancy (MacKenzie et al. 2002, Royle and Nichols 2003) in a Bayesian framework in the western and central Great Basin separately. We built models for species with >30 detections in ≥4 of the 9 survey years in the western Great Basin (24 species) and >30 detections in ≥10 of the 18 survey years in the central Great Basin (23 species).
We compared three occupancy models for each species. The first included data from all points sampled along the full elevational gradient. The second and third examined occupancy in the lowest and highest 25% (lower and upper edges ) of the gradient. The three models had the same formulation and covariates. All continuous covariates were scaled and centered around zero. Model selection techniques differed between the full-gradient model and the edge models. We applied indicator variable selection to the full-gradient model and the Watanabe-Akaike Information Criterion (WAIC) to the edge models, which did not converge when indicator variable selection was applied.
We modeled detection probability as
Cijk ~ \(\left\{\par \begin{matrix}\text{Bernoulli}\left(p_{\text{ijk}}\right)\text{\ \ if\ }Z_{\text{ik}}=1\\ 0\ \ \ if\ Z_{\text{ik}}=0\\ \end{matrix}\right.\ \)
\begin{equation} \text{logit}\left(p_{\text{ijk}}\right)=\beta 1+\beta 2*\text{jday}_{\text{ijk}}+\beta 3*\text{time}_{\text{ijk}}+\beta 4*{\text{time}^{2}}_{\text{ijk}}+\alpha 1*\text{observer}_{\text{ijk}}\nonumber \\ \end{equation}\begin{equation} \alpha 1\ \sim\ Normal\left(0,\ tau1\right)\nonumber \\ \end{equation}\begin{equation} tau1\ \sim\ Uniform\left(0,5\right),\nonumber \\ \end{equation}
where Cijk is the observed presence or absence of the species at point i during visit j in year k , and p is the probability of detecting a species given its presence. We used a logit link function to model four detection covariates: Julian date (jday ), time of day and its quadratic transformation, and a random, observer-level effect. β1 is the mean point-level detection probability for a given species, andα1 is a random effect of observer identity on detection, with a mean of 0 and a precision of tau1 .
We connected the detection process to the occupancy process throughZik , which we treated as a Bernoulli random variable governed by the success probability ψ :
\begin{equation} Z_{\text{ik}}\sim\ Bernoulli(\psi_{\text{ik}})\nonumber \\ \end{equation}\begin{equation} \text{logit}\left(\psi_{\text{ik}}\right)=\ \beta 5+\beta 6*\mathbf{X}_{\mathbf{\text{ik}}}+\alpha 2*\text{point}_{\text{ik}}+\alpha 3*\text{canyon}_{\text{ik}}\nonumber \\ \end{equation}\begin{equation} \alpha 2\ \sim\ Normal(0,\ tau2)\nonumber \\ \end{equation}\begin{equation} tau2\ \sim\ Uniform(0,1)\nonumber \\ \end{equation}\begin{equation} \alpha 3\ \sim\ Normal(0,\ tau3)\nonumber \\ \end{equation}\begin{equation} \text{tau}3\ \sim\ Uniform(0,1)\nonumber \\ \end{equation}
where Zik is the occupancy state (0 = absent, 1 = present) at point i in year k. We applied a logit link function to ψ to model occupancy covariates.Xik represents the vector of covariate values at point i in year k . Covariates in X included year, elevation, the interaction of year and elevation, spring temperature, winter precipitation, spring precipitation, and NDVI. If the interaction of year and elevation was included in the best model, and its posterior density estimate did not overlap zero, we concluded that the mean elevational distribution of the species had shifted upslope or downslope across years. We included random effects on the occupancy process to account for unmeasured differences among points. α2 is a point-level random effect with a mean of 0 and precision of tau2 , and α3 is a canyon-level random effect with a mean of 0 and precision oftau3 . We used vague prior distributions for intercepts, covariates, and random effects. All covariates were centered and scaled.
We conducted simulations to test whether species’ elevational distributions shifted linearly over time or were consistent with stochastic processes, such as random chance or population variability. For each species for which a significant interaction was included in the best model, we constructed a null distribution for the time trend by randomly permuting the year variable in 100 copies of the dataset. The model previously identified as best was estimated on these simulated data sets. We then conducted a two-tailed test of whether the estimated mean of (year x elevation) in each simulation was greater than the absolute value of the estimate from the original model. For example, if the estimated effect of year x elevation on occupancy was 0.65, we determined the percentage of simulations in which the mean estimate of the interaction term was greater than 0.65 or less than -0.65. We would expect >90% of simulations to include a mean estimate that was equal to or greater than the absolute value of the observed estimate if the effect was not a product of stochasticity.
We implemented all models in JAGS (Plummer 2003) with the jagsUI package (Kellner 2019) in R (R Core Team 2020). We based posteriors on three chains of 50,000 iterations after a 10,000 sample burn-in and adaptive phase. We classified convergence as Rhat <1.15 (Gelman and Hill 2007). No pairs of candidate covariates were collinear. We examined model fit on the basis of separate Bayesian p-values for the detection and occupancy processes, and we estimated mean occupancy and detection. We classified the fit of models as good if both of the Bayesian p-values were 0.05-0.95 and estimated mean detection and occupancy were 0.15.
We implemented a simple linear model of the effect of year on the mean elevation at which a species was observed. The resulting slope and intercept allowed us to calculate the average elevational distance that the species’ distribution shifted over the survey period, which was not possible to estimate from the results of the occupancy model. We used standard error estimates to calculate the 95% confidence interval of the elevational shift. To determine whether spring temperature, winter precipitation, spring precipitation, or NDVI changed over the survey period, we used generalized linear models (GLMs). We included a point-level random effect and modeled all variables as a function of elevation, year, and their interaction.