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.