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).