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.