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.