Statistical analysis
A statistical model was developed to assess whether the public-sourced photographs provide evidence of long-term and/or wide-spread variation in patterns of mountain goat shedding. We assumed that the progression of the amount of coat shed within a season could be described by the logistic equation. Let f (t ) denote the mean fraction of coat shed on day of year (DOY), t . In its simplest form our model is given by:
\begin{equation} f\left(t\right)=\frac{e^{\alpha(t-\tau)}}{1+e^{\alpha(t-\tau)}}\nonumber \\ \end{equation}
where τ is the DOY when 50% of the coat has been shed, andα describes the rate of shedding. Rate of change in f has maximum value α /4 and occurs when t = τ . We refer to τ and α as the shedding date and shedding rate, respectively. These two parameters may be affected by animal state or environmental variables. Here, the state of an animal is defined by its sex and, if female, the presence of a kid. Although animal state is known for the YWP study it is often not clear from community photographs. Animal state is described by a letter pairing; the first letter describes the sex of the animal (F = female, M = male, X = unknown) and the second letter describes the presence of kid (Y = yes, N = no, X = unknown). Assuming only females may be associated with kids, there are six animal states, three of which are unambiguous: FN (female without a kid), FY (female with a kid), and MN (male without a kid). We assume that animal state may affect the timing of shedding. Letτ 0 be the shedding date of a female without kid, and suppose the shedding date of males differs to females byτ M. Shedding date of females when with kid differs by τ K.
Shedding date and shedding rate may also correlate with elevation (E) and latitude (L), due to their relation with temperature and photoperiod. Specifically, for an animal photographed in year yat elevation x E and latitudex L, we assume
\(\alpha=\alpha_{0}+\alpha_{E}x_{E}+\alpha_{L}x_{L}+\alpha_{Y}y+\alpha_{[y]}\),
and
\(\tau=\tau_{0}+\tau_{M}x_{M}+\tau_{K}x_{K}+\tau_{E}x_{E}+\tau_{L}x_{L}+\tau_{Y}y+\tau_{[y]}\),
where x M and x K are binary predictors indicating whether the animal is male, or with kid (e.g., if animal state is FY then x M = 0 andx K = 1). Year factors into the model as both a continuous predictor and a random effect. Parametersα Y and τ Y describe long-term, smooth trends in shedding date and rate, which we might expect to differ from zero under climate change. Alternatively, theα [y ] and τ [y ] are random effects associated with year y , drawn from normal distributions with mean zero and standard deviation σα andστ , respectively, and represent any year-specific stochastic effects on the timing of shedding that are common to all locations (e.g., large-scale weather fluctuations).
We fit our model to the observed estimates of shedding fractions by first converting the fractions into integers. If fraction f had been shed then the response variable is set to n = round(fN ), where N is the number of shedding bins, implying 0 < n < N . We chose N= 25 bins; our statistical findings are robust to the choice ofN . The probability of observing n bins shed at timet when the expected fraction shed is f (t ), is given by the beta-binomial distribution, denotedP BB(n ). The beta-component accounts for over-dispersion in the observations relative to the binomial distribution (e.g., due to observation error, animal differences, or unknown environmental variation). We formulate the beta-binomial using the parameter φ (common to all observations) so that its variance inflation factor, relative to the binomial distribution, is\(v=1+(N-1)\phi/(1+\phi)\) (Richards, 2008).
The model described above can be fit to observations of shedding where animal state is unambiguous. We can also fit this model to observations when animal sex or presence of kid is unknown. Suppose at any time proportion p of animals are female and proportion q of the females are associated with a kid. These assumptions imply that the average proportion of animals in states: FN, FY, and MN, arep (1-q ), pq , and (1-p ), respectively. When animal state is ambiguous (i.e., FX, XN, or XX) the probability of observing n shed bins is a weighted sum of the three beta-binomial distributions associated with the unambiguous states, where the weights are calculated using p and q . Specifically,
\begin{equation} \Pr\left(n\middle|\text{FX}\right)=\left(1-q\right)P_{\text{BB}}\left(n\middle|\text{FN}\right)+qP_{\text{BB}}\left(n\middle|\text{FY}\right),\nonumber \\ \end{equation}\begin{equation} \Pr\left(n\middle|\text{XN}\right)=\frac{p\left(1-q\right)P_{\text{BB}}\left(n\middle|\text{FN}\right)+(1-p)P_{\text{BB}}\left(n\middle|\text{MN}\right)}{p\left(1-q\right)+(1-p)},\nonumber \\ \end{equation}
and
\begin{equation} \Pr\left(n\middle|\text{XX}\right)=p\left(1-q\right)P_{\text{BB}}\left(n\middle|\text{FN}\right)+pqP_{\text{BB}}\left(n\middle|\text{FY}\right)+(1-p)P_{\text{BB}}\left(n\middle|\text{MN}\right),\nonumber \\ \end{equation}
where P BB(n |j ) is the probability of observing n bins shed when the animal is in an unambiguous state j , which is calculated according to the beta-binomial distribution, as described above.
This model has 15 parameters (see Table 1) and we estimate them using Bayesian methods based on Monte Carlo sampling. Specifically, we implemented the model using the R programming language and the package rstan (Supplementary Materials). We specified relatively uninformative priors for all parameters so that the posterior distributions were strongly dependent on the data (Supplementary Materials). We used three sampling chains to visually check that our model formulation converged and posterior parameter distributions were based on 1000 samples after a 1000 sample burn-in.
For the YWP portion of our study we fit a simpler form of the model. As this study was conducted at a single site, and within a single season, we did not estimate effects of elevation, latitude or year (i.e.α E = α L =α Y = τ E =τ L = τ Y = 0). An important difference with sampling design between both studies is that animals were identified and repeatedly surveyed in the YWP study. To account for repeated measures, we used a random effect term associated with each animal; however, given the limited number of animals in the study (14) we only included a single random effect and associated it with shedding date. These random effect terms were drawn from a normal distribution with mean zero and standard deviation, σ ID. The model of shedding that we fit to the YWP study had six parameters (Table 2).