State-space modeling of population dynamics
We used a state-space model that includes a data model that represents
an observational process and a process model that represents the
underlying population dynamics. This approach allowed us to separate
observational and process errors.
Data model. The observed data vector
Yt, j contains the number of observed
grid points for C. dalli (YC, t, j ) andG. furcata (YG, t, j ), and those that were
bare or occupied by other species (Yother, t, j )
at quadrat j and census year t :
Yt, j = (YC, t,
j , YG, t, j , Yother, t,
j ). (1)
Population size (cm2 coverage) of C. dalli andG. furcata and the remaining plot area, including bare area and
that occupied by other species, are denoted as NG,
t, j , NC, t, j , andNother, t, j , respectively. Note that the sum of
these terms equals the total area of each quadrat (1000
cm2). We assumed that the observed data
Y follow a multinomial distribution (Multi), in which
the cell probabilities are the proportional coverage of each element:
Yt,j ~ Multi(40,
Nt, j / 1000), (2)
where Nt, j =(NC,t,j , NG,t,j ,Nother,t,j ) and the number of trials is equal to
the number of grid points (40).
Process model . For a logarithmic Gompertz population model
(Ives et al. 2003), a normally distributed stochastic term with variance
σ2 was applied as the population dynamics for each
species in each quadrat:
ln(NC,t,j ) ~N (rC,j + (1 −
αCC,j ) × ln(NC, t-1, j ) −
αCG,j × ln(NG, t-1, j ),
σC2) (3a)
ln(NG,t,j ) ~N (rG,j + (1 −
αGG,j ) × ln(NG, t-1, j ) −
αGC,j × ln(NC, t-1, j ),
σG2), (3b)
where rC and rG represent
the intrinsic growth rate of C. dalli and G. furcata ,
respectively. αCC and αGGare the strength of intraspecific density dependence for C. dalliand G. furcata , respectively; αCG is the strength
of interspecific density dependence of G. furcata on C.
dalli , and vice versa for αGC.
Parameter estimation. The Markov chain Monte Carlo (MCMC)
method was used to jointly estimate all parameters in all quadrats for
the census period 2002–2019, using the software JAGS through the R
package rjags and R2jags (Plummer 2015) in R v3.6.2 (R
Core Team 2013). We used independent and uninformative prior
distributions for the parameters. Priors for σC and
σG were specified as the uniform distribution (U), U(0,
100). Priors for intrinsic growth rates and strengths of density
dependence were specified as the normal distribution (N), N(μ,s ), where μ is the mean of the parameter that was specified as
the normal distribution N(0, 0.0001) and s is the variance of the
parameter that was specified as the uniform distribution U(0, 100).
We used four independent chains with 100,000 iterations each, and the
first 50,000 iterations were used as burn-in iterations to ensure that
the chains had converged. In addition, we thinned the chains by a factor
of 10 to reduce autocorrelation in the posterior samples and to produce
a reasonable amount of output. We used the Gelman and Rubin Ȓ
convergence diagnostics (Gelman & Rubin 1992) and visual inspection of
the chains to ensure convergence. The value of Ȓ of all parameters was
between 1.001 and 1.036, indicating that the posterior distribution of
all parameters in all quadrats adequately converged (1.000 < Ȓ
< 1.050; Gelman & Rubin 1992).
The JAGS code of the model is given in the Supporting Information.