Data analysis
We used multinomial logistic regression to estimate the probability of a
hare being in moult category i (white, moulting, or brown) on
each ordinal day d . We included year k specific intercepts
to test if moult timing varied between years. The model used the
following equation:
\begin{equation}
p\left(y=i\right)=\frac{e^{\alpha_{\text{ik}}+\beta_{1i}\times d+\beta_{2i}\times a_{j}{+\beta}_{3i}\times l_{j}{+\beta}_{4i}\times c_{j}+s_{\text{ij}}}}{1+\sum_{m=1}^{i-1}e^{\alpha_{i}+\beta_{1i}\times d+\beta_{2i}\times a_{j}{+\beta}_{3i}\times l_{j}{+\beta}_{4i}\times c_{j}+s_{\text{ij}}}}\nonumber \\
\end{equation}Where d is ordinal day, and aj ,lj , and cj , respectively,
are the altitude, latitude, and climate zone at camera site j (as
adapted from Zimova et al 2019 and 2020b). Parametersβ1-4j represent the slopes for the different
covariates. We included a category i and year k specific
intercept αik as well as a site-specific random
intercept sij to correct for multiple
observations per camera site. The brown category was set to 0 in both
spring and autumn models to provide a baseline for comparison with
moulting and white categories.
We implemented the multinomial logistic regression models in a Bayesian
framework using JAGS (Denwood 2016) called from R V4.1.3 (R Core Team
2022) with the jagsUI package (version 1.5.2) (Kellner 2021). We
standardized altitude, latitude, and climate zone (mean = 0, SD = 1)
before running the models. Ordinal day was included as an explanatory
covariate to enable estimation of the probability of a hare being white,
brown, or moulting between coats on specific days of the year. We
checked for collinearity between the covariates using both the variance
inflation factor (VIF) and Pearson correlation coefficient. For every
covariate combination the Pearson values were below 0.6 (Suppl 6) and
VIF values were below 2.0 (Suppl 7).
We used uninformative, normally distributed priors with a mean of 0 and
precision of 0.01 for all slopes and the year k specific
intercepts. For the site s random intercept, we used a mean of 0
and a standard deviation defined as a vague prior with a uniform
distribution between 0 and 100. We ran the model with three chains all
thinned by 100 for 120,000 iterations, with a burn-in of 60,000
iterations. We confirmed model convergence using traceplots and the
Gelman-Rubin convergence statistic (R-hat) (Brooks and Gelman 1998) with
all variables used in the final models having R-hat values of 1.10 or
less. Additional models in which the dataset was subset to only include
cameras located south of 61° North were used to test if model
performance was affected by camera trap placement north of this
latitude. The results were consistent with those obtained using the full
dataset indicating that camera placement north of 61° did not affect
model performance. We produced all figures using the ggplot2(Wickham 2016), raster (Hijmans 2022), and cowplot (Wilke
2020) packages.