Running Head
Spatial dynamic multistate model
Title
A multistate dynamic site occupancy model for spatially aggregated sessile communities
Authors
Keiichi Fukaya1
The Institute of Statistical Mathematics, 10-3 Midoricho, Tachikawa, Tokyo 190-8562 Japan
J. Andrew Royle
USGS Patuxent Wildlife Research Center, 12100 Beech Forest Road, Laurel, Maryland 20708 USA
Takehiro Okuda
National Research Institute of Far Seas Fisheries, Japanese Fisheries Research and Education Agency, 2-12-4 Fukuura,
Kanazawa-ku, Yokohama, Kanagawa, 236-8648 Japan
Masahiro Nakaoka
Akkeshi Marine Station, Field Science Center for Northern Biosphere, Hokkaido University, Aikappu, Akkeshi, Hokkaido
088-1113, Japan
Takashi Noda
Faculty of Environmental Earth Science, Hokkaido University, N10W5, Kita-ku, Sapporo, Hokkaido 060-0810, Japan
Summary
1. Estimation of transition probabilities of sessile communities seems easy in principle
but may still be difficult in practice because resampling error (i.e., a failure to
resample exactly the same location at fixed points) may cause significant estimation
bias. Previous studies have developed novel analytical methods to correct for this
estimation bias. However, they did not consider the local structure of community
composition induced by the aggregated distribution of organisms that is typically
observed in sessile assemblages and is very likely to affect observations.
2. We developed a multistate dynamic site occupancy model to estimate transition
probabilities that accounts for resampling errors associated with local community
structure. The model applies a nonparametric multivariate kernel smoothing method-
ology to the latent occupancy component to estimate the local state composition
near each observation point, which is assumed to determine the probability distri-
bution of data conditional on the occurrence of resampling error.
3. By using computer simulations, we confirmed that an observation process that de-
pends on local community structure may bias inferences about transition probabil-
ities. By applying the proposed model to a real dataset of intertidal sessile commu-
nities, we also showed that estimates of transition probabilities and of the properties
of community dynamics may differ considerably when spatial dependence is taken
into account.
1kfukaya@ism.ac.jp
1
4. Results suggest the importance of accounting for resampling error and local com-
munity structure for developing management plans that are based on Markovian
models. Our approach provides a solution to this problem that is applicable to
broad sessile communities. It can even accommodate an anisotropic spatial corre-
lation of species composition, and may also serve as a basis for inferring complex
nonlinear ecological dynamics.
Key words
Classification error; Community dynamics; Hierarchical models; Intertidal; Kernel smooth-
ing; Site occupancy models; Spatial correlation; Transition probability
2
Introduction
Markov models are a general class of mathematical models that are widely used to
describe dynamics of ecological systems. In the context of community ecology, they offer
a simple and useful representation of community dynamics. For example, a very basic
Markov model summarizes the dynamics of a community (i.e., changes in community
composition or site occupancy dynamics among ecological states) with a transition
matrix, P, which consists of the transition probability from state k to j, pjk, in the
element in its j th row and k th column. Such linear Markov models may not provide a
fully realistic description of community dynamics in nature, which may be inherently
nonlinear (Spencer & Tanner 2008; Tanner et al. 2009). Nevertheless, they can still
provide a good approximation of observed community dynamics and composition
(Wootton 2001a,b; Hill et al. 2002). Linear Markov models also allow us to assess a
wide range of properties of species-level dynamics (e.g., rate of colonization, disturbance,
and replacement) as well as those of community-level dynamics (e.g., equilibrium
community composition, as well as its sensitivities, mean turn over time, and damping
ratio), which can be derived from the transition probability matrix (Hill et al. 2004).
Markov community models have been applied to communities of sessile organisms
(e.g., terrestrial plants, corals and intertidal/subtidal communities) because in those
systems transition probabilities can be estimated with field data by tracking species
occupancy at many fixed observation points over multiple periods of time (e.g., Wootton
2001a,b; Hill et al. 2002, 2004; Tanner et al. 2009; Fig. 1A). However, estimation of
transition probabilities may be biased when there are resampling errors, that is, failures
to resample exactly the same fixed-point locations (Conway-Cranos & Doak 2011).
When such resampling error occurs, researchers observe a point that is different from,
but probably close to, the correct location of the fixed point; the result may be a
“misclassification” of the occupancy state. Such an error may sometimes be inevitable
in field sampling, where somewhat crude tools are necessarily used and/or organisms are
small (Conway-Cranos & Doak 2011). New analytical methods have been proposed to
3
correct this estimation bias which assume that the probability distribution of data,
conditional on the occurrence of resampling error, is related to the relative state
frequency of ecological states (Conway-Cranos & Doak 2011; Fukaya & Royle 2013). On
the one hand, Conway-Cranos & Doak (2011) proposed a maximum likelihood approach,
in which transition probabilities and resampling error rates are estimated based on
multinomial likelihoods. On the other hand, Fukaya & Royle (2013) proposed to use the
dynamic site occupancy modeling framework to account for resampling error in the
estimation of transition probabilities.
In both of these studies, spatial aspects of community structure were not considered
explicitly. It is very likely, however, that the local structure of community composition
affects the observation process described above, because sessile organisms typically
represent locally aggregated spatial distributions. Therefore, accounting for such spatial
dependence will be necessary to make more reliable inferences about underlying
community dynamics. In this study, we developed an extension of the multistate
dynamic site occupancy model proposed by Fukaya & Royle (2013) to account for
spatial dependence in the estimation of transition probabilities. We used additional
information about the spatial coordinate of each observation point to estimate the state
composition near each point. Estimation of this local community structure was possible
by applying the nonparametric multivariate kernel smoothing methodology (Diggle
et al. 2005) to the latent occupancy component. Our approach can even account for an
anisotropic spatial correlation of species composition and may also lead to some
potential extensions of the model, which we describe in the Discussion. By using
computer simulations, we confirmed that an observation process that depends on local
community structure may cause some bias in the inference of transition probabilities.
By applying the proposed model to a real dataset of intertidal sessile communities, we
also showed that estimates of transition probabilities and of community dynamics may
be considerably different when spatial dependence is accounted for.
4
The model
Model description
We assumed that there were, in total, I fixed observation points within a permanent
quadrat that could be occupied by one of the possible S ecological states (Fig. 1A).
Hence, observation points and the quadrat corresponded to “sites” (local populations)
and the “meta-population”, respectively, in the meta-population design concept for the
site occupancy modeling framework (Kéry & Schaub 2012; Kéry & Royle 2016). The
ecological states are typically defined by dominant species or groups of species within the
community. “Free space” may also commonly be considered. The occupancy state was
assessed for each site over T periods of time. We define zit and yitn (zit, yitn ∈ 1, . . . , S)
as the occupancy state at site i (i = 1, . . . , I) at time t (t = 1, . . . , T ), and the state
observed by a researcher (i.e., data recorded) at the nth measurement
(n = 1, . . . , N(i, t)) at site i at time t, respectively. The N(i, t) observations consisted of
replicated samples at site i and time t; these samples may have been collected by
repeated surveys conducted within a sufficiently short period or by independent
observers. This sampling scheme formally resembles Pollock’s robust design, which is
classically considered in capture-recapture studies (Pollock 1982) in which secondary
resampling is performed within a single primary sampling period (t). We note that the
dynamic site occupancy modeling framework described below also permits missing data.
Changes in site occupancy states (zit) across primary sampling periods are described
with a transition probability matrix P. The element in the j th row and k th column
pjk(1 ≤ j, k ≤ S) represents the transition probability from state k to j. Resampling
error probability e is also considered to account for a certain type of observation error
that arises when there is failure to resample the exact location of the fixed observation
site (Conway-Cranos & Doak 2011). It is thus assumed that when resampling error does
not occur at a sampling occasion for site i at time t (occurring with a probability 1− e),
the occupancy state zit is recorded with probability one, whereas when resampling error
occurs, the “occupancy state” observed is a random variable that follows a certain
5
probably distribution (Fig. 1B). Fukaya & Royle (2013) assumed that when a
resampling error occurs, the probability of state s being observed is equal to the relative
dominance of that state within the quadrat, fts =
∑
i 1(zit = s), where 1(x) is an
indicator function. Conceptually, this assumption implicitly requires that organisms are
homogenously distributed within the quadrat so that the state composition is identical
over all sites. The same assumption was made by Conway-Cranos & Doak (2011) to
account for a similar type of resampling error in longitudinal observations of sessile
communities. To account for the spatial dependence of community structure, we here
consider the relative dominance of the state near each site, gits, and relate it, instead of
fts, to the data distribution conditional on the occurrence of the resampling error.
Although the model framework could accommodate more ecological realities, we here
restricted our model by using the simplest elements possible for notational simplicity
and clarity. We believe this approach will facilitate understanding the basic model
structure. Possible extensions of the model are described in the Discussion. Formally,
the model we propose is as follows.
Observation model — We assume that on each sampling occasion, with probability
(1− e), observers assess the exact location of site i and find the true occupancy state zit.
But otherwise (i.e., with probability e), they fail to observe the exact location of that
site. In the latter case, we assume that they record state s with a probability that
equals the local dominance of that state, gits. This observation distribution, conditional
on the occurrence of the resampling error, represents a categorical distribution with a
probability vector (git1, . . . , gitS). With a latent indicator variable mitn representing the
occurrence of resampling error for the nth measurement at site i at time t, the model for
this observation process is expressed as:
yitn = zit when mitn = 0
yitn ∼ Categorical(git1, . . . , gitS) when mitn = 1.
(1)
We assume that resampling errors occur independently, and thus that mitn simply
follows a Bernoulli distribution:
6
mitn ∼ Bernoulli(e). (2)
Process model — Assuming that the dynamics of the occupancy state is a Markov
process, the conditional probability of latent state zit for t = 2, 3, . . . , T can be expressed
as follows:
zit|(zi,t−1=s) ∼ Categorical(p1s, p2s, . . . , pSs), (3)
where pjs gives the transition probability from state s to j. The transition probability
matrix P is expressed as follows:
P =

p11 . . . p1s . . . p1S
...
. . .
...
...
ps1 . . . pss . . . psS
...
...
. . .
...
pS1 . . . pSs . . . pSS

(4)
where column s gives the vector of transition probabilities for zit|(zi,t−1=s). Note that
this transition probability matrix is column-stochastic (i.e., each column sums to unity).
For t = 1, the initial occupancy probability of each site is defined by a probability
vector ϕ = (ϕ1, . . . , ϕS). The model for initial occupancy states is thus expressed as:
zi1 ∼ Categorical(ϕ1, . . . , ϕS). (5)
A kernel estimator for the local state composition — The local state
dominance gits is unobserved, and we now describe a model for gits. We estimate gits
nonparametrically by using a multivariate kernel smoothing methodology to estimate
the spatial structure of multinomial probabilities (Diggle et al. 2005).
We let xi be a two-dimensional coordinate vector for site i. For each i, t, and s, the
kernel regression estimator for the local state composition is expressed as follows:
gits =
∑I
j=1KΣ(xi,xj)1(zjt = s)∑S
r=1
∑I
j=1KΣ(xi,xj)1(zjt = r)
, (6)
where zjt is the latent occupancy state at site j at time t, and KΣ(xi,xj) is a
7
two-dimensional Gaussian kernel function with a bandwidth matrix Σ:
KΣ(xi,xj) = exp
{
−1
2
(xi − xj)′Σ−1(xi − xj)
}
. (7)
The bandwidth matrix Σ is a 2× 2 positive definite covariance matrix:
Σ =
 σ21 ρσ1σ2
ρσ1σ2 σ
2
2
 , (8)
where σ1, σ2 > 0 is the scale parameter for each dimension and −1 < ρ < 1 is the
correlation parameter. Note that the assumption of a common kernel function across
ecological states simplifies the estimator (eqn 6) as follows:
gits =
∑I
j=1KΣ(xi,xj)1(zjt = s)∑I
j=1KΣ(xi,xj)
. (9)
With the kernel estimator of this form, gits becomes higher at site i as the number of
sites that are occupied by state s at time t increases. Although all the sites occupied by
s at time t contributes to gits, the weight decreases as a function of the distance from
site i: the scale of this distance dependence is controlled by the bandwidth matrix Σ.
We note that Σ accounts for an anisotropic spatial dependence (i.e., direction
dependence) of the local community structure: the spatial dependence is isotropic only
when σ1 = σ2 and ρ = 0.
As long as there are no significant sources of observation error other than the
resampling error we considered above, the estimated vector (git1, . . . , gitS) is interpreted
naturally as a representation of the composition of ecological states surrounding site i at
time t. Equivalent to fts in the model of Fukaya & Royle (2013), gits is a derived
parameter obtained as a function of latent occupancy state {zit}, whereas the function
defining gits now depends on an unknown bandwidth matrix with parameters to be
estimated. We also note that when the bandwidth of the kernel is sufficiently large, gits
will no longer vary over sites within the quadrat. The result is a spatially homogeneous
situation that is equivalent to the model considered by Fukaya & Royle (2013).
8
Estimation of states and parameters
In this model, elements of the transition probability matrix P, the resampling error rate
e, the initial occupancy probability vector ϕ, and the elements of the bandwidth matrix
Σ are parameters to be estimated. The model also involves occupancy state zit, error
occurrence indicator mitj, and the local state composition gits (a derived quantity) as
latent state variables. The model assumes that observations at a particular time point
depend on all the state variables at that time. Such a class of dynamic models is called
a factorial hidden Markov model in the machine learning community, and for such a
model it is known that an exact likelihood inference is computationally intractable
(Ghahramani & Jordan 1997). We thus adopt here a Bayesian approach for parameter
inference, in which the joint posterior distribution of parameters and latent state
variables is obtained by using Markov chain Monte Carlo (MCMC) methods. We specify
vague priors for parameters for the observation and process models:
e ∼ Beta(1, 1), (10)
(p1s, . . . , pSs) ∼ Dirichlet(1, . . . , 1), (11)
(ϕ1, . . . , ϕS) ∼ Dirichlet(1, . . . , 1), (12)
where the diffuse Dirichlet prior (eqn 11) is specified for each column of the transition
probability matrix. Vague priors would also be specified for elements of the bandwidth
matrix, Σ, as follows:
σ1 ∼ Uniform(0, U) (13)
σ2 ∼ Uniform(0, U) (14)
ρ ∼ Uniform(−1, 1) (15)
where U is a sufficiently large positive constant defining the upper limit of the uniform
prior.
The joint posterior distribution of these parameters and latent variables ({zit}, {mitn},
9
{gits}) can be obtained by using the MCMC methods. The model can be implemented
using the BUGS software. We have provided a model script in the supplemental
material that can be used to fit the model using JAGS software (Plummer 2003).
Simulation study
Materials and methods
To examine the effect of an observation process that is relevant to the local community
structure on the inference of community dynamics, we conducted a simulation study.
Using the system and observation models described above, with fixed model parameters
that were to be estimated, we simulated a number of replicate datasets that were
obtained from the hierarchical data-generating process described in the previous section.
Simulations and analyses described in what follows were conducted using R (versions
3.1.0 to 3.2.5). An R script for this simulation is provided in the supplemental material.
Three models were fitted to these simulated data sets. The first was a classical, naive
estimator for transition probabilities, defined as p̂jk = njk/
∑
l nlk, where njk is the
number of sites that were in state k at a certain point in time and in state j at the
following time (Spencer & Susko 2005). This model does not account for any type of
observation error. The second was a “non-spatial” multistate dynamic site occupancy
model proposed by Fukaya & Royle (2013), which does account for resampling errors.
The model assumes that the distribution of the observed state, conditional on the
occurrence of the resampling error, is proportional to the relative dominance of each
state in the quadrat that is homogeneous in space. Finally, the “spatial” multistate
dynamic site occupancy model described in the previous section was also fitted to
simulated datasets. Note that this model is equivalent to the data-generating model and
was thus expected to perform the best among the three models.
To generate replicated datasets, we set the number of sites I = 225, the length of time
series T = 5, and the number of ecological states S = 5. The sites were assumed to be
aligned on a 15× 15 grid, and their coordinates were represented as
10
x1 = (1, 1),x2 = (2, 1), . . . ,x225 = (15, 15). We used a transition probability matrix that
was randomly drawn from a diffuse Dirichlet distribution
(p1s, . . . , pSs) ∼ Dirichlet(1, . . . , 1). For the observation process, we assumed an isotropic
bandwidth kernel by setting σ1 = σ2 = 1 and ρ = 0 (an anisotropic case is also examined
and reported in Appendix S2). With these settings, we generated 96 replicate datasets
for each of 6 levels of error rate (e = 0, 0.15, 0.3, 0.45, 0.6, 0.75).
We considered two scenarios for the manner of data collection. The first scenario
considered a situation of limited information where no replicated observations were
obtained; that is, we set N(i, t) = 1 for all i and t. The second scenario simulated a
more ideal situation where three replicated observations were conducted for each survey
(i.e., the robust design); that is, we set N(i, t) = 3 for all i and t.
Among the three models and two sampling scenarios described above, the
performances of the models were evaluated based on the mean square error (MSE),
square of bias (Bias2) and variance (Var) of the point estimate of each parameter. For
transition probabilities, the above criteria were averaged over all elements in the
probability matrix to evaluate the magnitude of the estimation errors in a sequence of
transition probabilities. Details of the procedure for parameter estimation are described
in Appendix S1.
Results
The results of the simulation study are summarized in Fig. 2. We first note that when
the error rate is small, the posterior distribution for the bandwidth parameter σ in the
spatial model sometimes failed to converge. In such cases, the posterior distribution of σ
was typically distributed over a wide range, a reflection of the flat specified prior
distribution, the result being spatially homogeneous estimates of the local state
composition gits. This result suggests an identification issue for σ in these settings, an
issue that appears to become worse when replicated samples are lacking (Fig. 2). In the
spatial model, σ determines the probability distribution of data conditional on the
11
occurrence of resampling error. Hence, σ is not satisfactorily determined by data when
resampling errors occur infrequently. Even when the posterior sampling for σ was
unsuccessful, however, the convergence of other parameters in the spatial model (P and
e) was typically achieved. In these cases, estimates of transition probabilities were often
very similar in the spatial model and the non-spatial model. This similarity reflects the
fact that the estimated local state composition gits was spatially homogeneous, the
result being a spatial model that was virtually the same as the non-spatial model, in
which a spatially homogeneous state composition was assumed implicitly.
For this reason, we excluded estimates from the calculation of MSE, bias, and
variance of σ̂, if the difference between the point estimate and the true value of σ was
more than 10 (Fig. 2, lower panels). We omitted the results of MSE, square bias, and
the variance of σ̂ for e = 0 because this procedure excluded almost all estimates for
e = 0. We also note that in the spatial model with e = 0 and three replicated
observations, there was an estimate of e that was exceptionally large. We also excluded
that estimate from the calculation of MSE, bias, and variance of ê.
Not surprisingly, overall, the bias was the smallest in the spatial model compared to
other models (Fig. 2). When comparisons were made among the results of the
single-replicate case, estimates of transition probabilities from the naive estimator were
the least variable (Fig. 2, upper panels). However, the naive estimator became
profoundly biased as the error rate increased, the result being the worst MSE of the
transition probabilities among the examined models. The bias and variance of the
transition probability estimates in the non-spatial and spatial model increased with the
error rate. The bias of the transition probabilities was considerably smaller in these
models compared to the naive model, resulting in the better MSE as the error rate
increased. Although the difference between the non-spatial and spatial models in the
estimation of transition probabilities is not so large, the bias and MSE of the spatial
model was smaller than that of the non-spatial model when the error rate was high (i.e.,
e ≥ 0.45). Replicated observations improved the estimation of transition probabilities in
the spatial and non-spatial model. Under these conditions, the performance of the
12
non-spatial model was even better than that of the spatial model with a single replicate.
The MSE of ê tended to be higher for the non-spatial model than for the spatial
model (Fig. 2, middle panels) because the bias of the non-spatial model increased faster
than that of the spatial model, although the variance of ê was in general smaller in the
non-spatial model than in the spatial model. In the non-spatial model, replicated
observations did not mitigate the MSE and bias of ê. In the spatial model, the MSE and
bias of σ̂ tended to be higher when the error rate was low (Fig. 2, lower panels), where,
as noted above, the estimation of σ may be difficult.
In summary, these results suggest that, when the local community structure affects
observations, both the naive estimator and the non-spatial multistate dynamic site
occupancy model may provide worse estimates than the spatial model. When the error
rate is low, the bandwidth parameter of the spatial model may not be satisfactorily
determined by data, but the model may still perform as well as the non-spatial model by
estimating the spatially homogeneous relative state composition. As the error rate
increased, the spatial model on average outperformed the non-spatial model. We note
that these results hold even when the underlying spatial correlation is anisotropic
(Appendix S2).
Application to real data
Materials and methods
We applied the proposed spatial model, in addition to the non-spatial model (Fukaya &
Royle 2013) and the naive estimator for transition probabilities, to an intertidal sessile
community dataset collected at a total of 25 permanent quadrats within five shores —
Mochirippu (MC), Mabiro (MB), Aikappu (AP), Monshizu (MZ) and Nikomanai
(NN)— located on the Pacific coast of eastern Hokkaido, Japan (detailed descriptions of
the study sites and biogeographic features of the area can be found in Okuda et al. 2004,
Nakaoka et al. 2006 and Fukaya et al. 2014). Each quadrat contains 200 sites at which
occupancy states were assessed once per year from 2002 to 2012.
13
For each quadrat, models were fitted independently to data in a fashion analogous to
that used in the simulation study described in the previous section. To fit the spatial
model, we restricted the bandwidth matrix to Σ = diag(σ21, σ
2
2), where σ1 and σ2 were
the scale parameters along the vertical and horizontal directions, respectively.
As derived parameters, we also obtained estimates of two quantities that inform some
dynamical aspects of communities: the mean turnover time and the damping ratio (Hill
et al. 2004). The former quantifies the average turnover time of a site randomly selected
from the stationary community. It is calculated as
∑S
s=1ws(1− pss)−1 where ws is the
equilibrium relative dominance of state s obtained as the normalized dominant right
eigenvector of the transition probability matrix P. The latter quantifies the lower bound
of the rate of convergence to equilibrium. It is equated to 1/|λ2| where λ2 is the
second-largest eigenvalue of P.
More details of the sampling design and model fitting are described in Appendix S3.
Results
As an example, parameter estimates obtained at a particular quadrat (MC2) for each
model are shown in Table 1. The estimated diagonal elements of the transition
probability matrix (i.e., the probability of persistence of the occupied state) tended to
be higher for the non-spatial model than those obtained with the naive estimator. These
estimated probabilities were even higher for the spatial model. This tendency was
generally found in other quadrats and reflects a consistent pattern of estimated
properties of community dynamics (Fig. 3). The estimated mean turnover time tended
to be longest in the spatial model, intermediate in the non-spatial model, and shortest
in the naive method. Furthermore, the damping ratio was estimated to be fairly high for
the naive method, whereas it was lower for the spatial and non-spatial models. These
results suggest that ignoring the resampling error and the effect of local community
structure on observations may lead to an overestimation of the “velocity” of community
dynamics because of an underestimation of persistence probabilities.
Interestingly, in the spatial model, the scale of bandwidth for the vertical axis was
14
estimated to be approximately half that for the horizontal axis (Table 1). This result
coincides with the known observation that rocky intertidal communities may be more
variable along the vertical than the horizontal direction because environments tend to
vary more rapidly along the vertical axis (Benedetti-Cecchi 2001; Valdivia et al. 2011).
The estimated resampling error rate tended to be higher in the spatial model than in
the non-spatial model (Table 1). Observations on these parameters noted here were also
generally found in other quadrats. We show in Fig. 4 a snapshot of the estimated local
community composition for a particular quadrat (MC2) and point in time.
Discussion
In ecological studies, consideration of spatial effects has been thought to be fundamental
(Legendre 1993; Lichstein et al. 2002; Royle et al. 2007). It has also been widely
recognized that accounting for observation processes in data collection and analyses is
important because ignoring these processes may lead to biased estimates and thus
misleading conclusions (Williams et al. 2002; Kéry & Schaub 2012; Kéry & Royle 2016).
These two considerations have been incorporated into the dynamic site occupancy
modeling framework to study metapopulation dynamics while accounting for imperfect
detection (e.g., Bled et al. 2011; Risk et al. 2011; Yackulic et al. 2012). As far as we
know, however, such spatial dynamic site occupancy models have been applied only to
population dynamics of single species (but see Yackulic et al. (2014)). The spatial
multistate dynamic site occupancy model we developed here would be widely applicable
to studying the dynamics of sessile communities in various systems, including corals,
mussels and terrestrial plants, for which precise resampling of observation points may be
difficult in the field for practical reasons (Conway-Cranos & Doak 2011).
For spatially referenced site occupancy data, dynamic site occupancy models with
autologistic components have been proposed (Bled et al. 2011, 2013; Yackulic et al.
2012, 2014; Eaton et al. 2014; see also Broms et al. 2016). In these models, the
adjacency of sites must be defined a priori, because probabilities of occupancy, local
15
colonization and extinction are assumed to depend on the occupancy states of
neighboring sites. It has also been typical in the application of these models that sites
are aligned on grids. However, as has been done in this study, another approach for
modeling spatial effects is to use kernel-weighting functions. Several dynamic site
occupancy models that were recently developed based on the metapopulation theory use
this approach to account for spatial dependency of local colonization and extinction
probabilities (Risk et al. 2011; Heard et al. 2013; Sutherland et al. 2014; Chandler et al.
2015). Because the scale of correlation is determined by the data, this modeling
approach does not require a priori determination of site adjacency. The observation
sites may also not need to be aligned in a reticular pattern. A related weighting function
approach has also been used to identify the spatial scale at which landscape variables
affect abundance (Chandler & Hepinstall-Cymerman 2016). We note, however, that a
possible difficulty in using kernel weighting may be its associated computational burden.
Evaluating the kernel weights for every combination of two sites may be
time-consuming, and the amount of computation time required increases very rapidly
with the number of sites being considered.
The multivariate Gaussian kernel used in the proposed model can account for an
anisotropic spatial correlation of community structure. This characteristic of the
multivariate Gaussian kernel can be useful when a known, or even unknown,
environmental gradient exists in the quadrat and affects the distribution of sessile
species. Although other types of kernels may be used to model the spatial structure, the
multivariate Gaussian kernel allows a straightforward implementation of anisotropic
correlation by specifying a bandwidth matrix with scale and correlation parameters.
However, results of our simulation study suggest that, in the proposed model,
bandwidth parameters may be difficult to estimate when the resampling error rate is
quite low. In such cases, the proposed spatial model would not be a reasonable option
for the inference of transition probabilities. In this regard, adoption of a robust design is
recommended because obtaining replicated samples may mitigate this problem (Fig. 2).
We note that dynamic occupancy models in general do not require equal numbers of
16
replications for every site and time.
Results from a dataset of intertidal communities on the Pacific coast of eastern
Hokkaido highlight the fact that different models may produce considerably different
estimates of transition probabilities as well as of the properties of community dynamics,
which are obtained from the transition probability matrix (Table 1, Fig. 3). We have
shown that, compared to the non-spatial model (Fukaya & Royle 2013), the spatial
model tends to yield a higher persistence probability of each state, longer mean turnover
time and lower damping ratio. These differences were even greater when the spatial
model was compared to the naive model. In the assessment of their intertidal sessile
community data from California, Conway-Cranos & Doak (2011) also found that the
estimated damping ratio was higher when resampling errors were not corrected. These
results have implications for the development of management plans that are based on
Markovian models. They suggest that ignoring the existence of resampling errors, as
well as local community structure, may lead to overestimation of the recovery rate of
ecological communities from natural or anthropogenic disturbance, the result being
poorly informed management strategies. We note, however, that little is known about
how the inference of ecological quantities, which are obtained as a function of transition
probabilities, may be affected by these factors. A further study is therefore required.
In the spatial model, resampling error rate was estimated considerably higher than
the nonspatial model (Table 1). This was a general tendency in the quadrats we
examined (mean and SD of point estimates of e: 0.74 ± 0.15). Although we do not
know if such a high level of error rate is typical, the existence of resampling error
appears to be the rule rather than the exception, at least in the field observation of
rocky intertidal communities. Conway-Cranos & Doak (2011) reported that, by using
their proposed method, the estimated error rate was 0.353 in the surveys of intertidal
assemblage in California. Different sampling protocols or types of communities should
result in different levels of error rate. In any case, however, we recommend accounting
for resampling error in some way rather than ignoring it, because the latter may lead to
biased inference.
17
The hierarchical formulation of the model allows us to readily extend the proposed
model, at least conceptually, to add more ecological realism (Royle & Dorazio 2008;
Kéry & Schaub 2012; Kéry & Royle 2016). For example, if some site- or time-specific
environmental covariates are available, they may be incorporated into the model to take
account of variations in transition probabilities and error rates. A possible extension of
the model that may have a particular ecological significance is to link estimated local
state dominance gits values to transition probabilities. Because interactions between
sessile organisms are limited to neighboring individuals, the rate of local colonization,
persistence, and replacement are expected to depend on surrounding local community
structure (Wootton 2001; Kawai & Tokeshi 2006). Although the dependency of
transition probabilities on species frequency has been considered in previous studies
about Markovian community dynamics (i.e., nonlinear Markov community model;
Spencer & Tanner 2008; Tanner et al. 2009), only the connection between transition
probabilities and the overall (i.e., global) species frequency has been explored. However,
by using the underlying occupancy component explicitly, dynamic site occupancy
models can model a relationship between local occupancy density and the rates of
colonization and/or extinction (Bled et al. 2011; Risk et al. 2011; Yackulic et al. 2012),
although such dependency has rarely been modeled in a multispecies context (but see
Yackulic et al. (2014) who explored intra- and inter-specific effects on site occupancy
dynamics of two owl species). Because the proposed model explicitly involves local
community structure as a model component, it will naturally provide a basis for
inferring complex nonlinear ecological dynamics.
Acknowledgements
We thank I. K. Shimatani and S. Eguchi for their helpful comments and discussion. We
also thank two anonymous referees for helpful reviews of the manuscript. For providing
access to laboratory facilities, we are grateful to the staff and students of the Akkeshi
Marine Stations of Hokkaido University. We acknowledge many researchers and students
18
who helped with our fieldwork. This research was supported by an allocation of
computing resources of the SGI ICE X and SGI UV 2000 supercomputers from the
Institute of Statistical Mathematics. Funding was provided by the Japan Society for the
Promotion of Science (Grants-in-Aid for Scientific Research No 15K18617 and
Grant-in-Aid for JSPS Fellows No 16J07614 to KF, Grants-in-Aid for Scientific Research
Nos 14340242, 18201043 and 21241055 to MN, and Nos 20570012, 24570012 and
15K07208 to TN). The authors declare no conflict of interest. Any use of trade, firm, or
product names is for descriptive purposes only and does not imply endorsement by the
U.S. Government.
Data Accessibility
Data from intertidal rocky shore communities, scripts for data analyses and results are
available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.30mb4.
Author Contributions Statement
KF and JAR conceived the ideas and designed methodology; KF, TO, MN and TN
collected the data; KF conducted simulations and data analyses; KF and JAR led the
writing of the manuscript. All authors contributed critically to the drafts and gave final
approval for publication.
References
Benedetti-Cecchi, L. (2001) Variability in abundance of algae and invertebrates at
different spatial scales on rocky sea shores. Marine Ecology Progress Series, 215,
79–92.
Bled, F., Nichols, J.D. & Altwegg, R. (2013) Dynamic occupancy models for analyzing
species’ range dynamics across large geographic scales. Ecology and Evolution, 3,
4896–4909.
19
Bled, F., Royle, J.A. & Cam, E. (2011) Hierarchical modeling of an invasive spread: the
Eurasian Collared-Dove Streptopelia decaocto in the United States. Ecological
Applications, 21, 290–302.
Broms, K.M., Hooten, M.B., Johnson, D.S., Altwegg, R. & Conquest, L.L. (2016)
Dynamic occupancy models for explicit colonization processes. Ecology, 97, 194–204.
Chandler, R. & Hepinstall-Cymerman, J. (2016) Estimating the spatial scales of
landscape effects on abundance. Landscape Ecology, 31, 1383–1394.
Chandler, R.B., Muths, E., Sigafus, B.H., Schwalbe, C.R., Jarchow, C.J. & Hossack,
B.R. (2015) Spatial occupancy models for predicting metapopulation dynamics and
viability following reintroduction. Journal of Applied Ecology, 52, 1325–1333.
Conway-Cranos, L.L. & Doak, D.F. (2011) Sampling errors create bias in Markov
models for community dynamics: the problem and a method for its solution.
Oecologia, 167, 199–207.
Diggle, P., Zheng, P. & Durr, P. (2005) Nonparametric estimation of spatial segregation
in a multivariate point process: bovine tuberculosis in Cornwall, UK. Journal of the
Royal Statistical Society: Series C (Applied Statistics), 54, 645–658.
Eaton, M.J., Hughes, P.T., Hines, J.E. & Nichols, J.D. (2014) Testing metapopulation
concepts: effects of patch characteristics and neighborhood occupancy on the
dynamics of an endangered lagomorph. Oikos, 123, 662–676.
Fukaya, K., Okuda, T., Nakaoka, M. & Noda, T. (2014) Effects of spatial structure of
population size on the population dynamics of barnacles across their elevational
range. Journal of Animal Ecology, 83, 1334–1343.
Fukaya, K. & Royle, J.A. (2013) Markov models for community dynamics allowing for
observation error. Ecology, 94, 2670–2677.
Ghahramani, Z. & Jordan, M.I. (1997) Factorial hidden Markov models. Machine
learning, 29, 245–273.
20
Heard, G.W., McCarthy, M.A., Scroggie, M.P., Baumgartner, J.B. & Parris, K.M.
(2013) A bayesian model of metapopulation viability, with application to an
endangered amphibian. Diversity and Distributions, 19, 555–566.
Hill, M.F., Witman, J.D. & Caswell, H. (2002) Spatio-temporal variation in Markov
chain models of subtidal community succession. Ecology Letters, 5, 665–675.
Hill, M.F., Witman, J.D. & Caswell, H. (2004) Markov chain analysis of succession in a
rocky subtidal community. American Naturalist, 164, E46–61.
Kawai, T. & Tokeshi, M. (2006) Asymmetric coexistence: bidirectional abiotic and biotic
effects between goose barnacles and mussels. Journal of Animal Ecology, 75, 928–941.
Kéry, M. & Royle, J.A. (2016) Applied Hierarchical Modeling in Ecology: Analysis of
Distribution, Abundance and Species Richness in R and BUGS. Volume 1: Prelude
and Static Models. Academic Press, London.
Kéry, M. & Schaub, M. (2012) Bayesian Population Analysis using WinBUGS: A
Hierarchical Perspective. Academic Press, London.
Legendre, P. (1993) Spatial autocorrelation: trouble or new paradigm? Ecology, 74,
1659–1673.
Lichstein, J.W., Simons, T.R., Shriner, S.A. & Franzreb, K.E. (2002) Spatial
autocorrelation and autoregressive models in ecology. Ecological Monographs, 72,
445–463.
Nakaoka, M., Ito, N., Yamamoto, T., Okuda, T. & Noda, T. (2006) Similarity of rocky
intertidal assemblages along the Pacific coast of Japan: effects of spatial scales and
geographic distance. Ecological Research, 21, 425–435.
Okuda, T., Noda, T., Yamamoto, T., Ito, N. & Nakaoka, M. (2004) Latitudinal gradient
of species diversity: multi-scale variability in rocky intertidal sessile assemblages along
the Northwestern Pacific coast. Population Ecology, 46, 159–170.
21
Plummer, M. (2003) JAGS: A program for analysis of Bayesian graphical models using
Gibbs sampling. Proceedings of the 3rd international workshop on distributed
statistical computing (DSC 2003), volume 124, p. 125. Technische Universit at Wien,
Austria.
Pollock, K.H. (1982) A capture-recapture design robust to unequal probability of
capture. Journal of Wildlife Management, 46, 752–757.
Risk, B.B., de Valpine, P. & Beissinger, S.R. (2011) A robust-design formulation of the
incidence function model of metapopulation dynamics applied to two species of rails.
Ecology, 92, 462–474.
Royle, J.A. & Dorazio, R.M. (2008) Hierarchical Modeling and Inference in Ecology: the
Analysis of Data from Populations, Metapopulations and Communities. Academic
Press, London.
Royle, J.A., Kéry, M., Gautier, R. & Schmid, H. (2007) Hierarchical spatial models of
abundance and occurrence from imperfect survey data. Ecological Monographs, 77,
465–481.
Spencer, M. & Susko, E. (2005) Continuous-time Markov models for species
interactions. Ecology, 86, 3272–3278.
Spencer, M. & Tanner, J.E. (2008) Lotka-Volterra competition models for sessile
organisms. Ecology, 89, 1134–1143.
Sutherland, C.S., Elston, D.A. & Lambin, X. (2014) A demographic, spatially explicit
patch occupancy model of metapopulation dynamics and persistence. Ecology, 95,
3149–3160.
Tanner, J.E., Hughes, T.P. & Connell, J.H. (2009) Community-level density dependence:
an example from a shallow coral assemblage. Ecology, 90, 506–516.
22
Valdivia, N., Scrosati, R.A., Molis, M. & Knox, A.S. (2011) Variation in community
structure across vertical intertidal stress gradients: how does it compare with
horizontal variation at different scales? PLoS ONE, 6, e24062.
Williams, B.K., Nichols, J.D. & Conroy, M.J. (2002) Analysis and Management of
Animal Populations. Academic Press, London.
Wootton, J.T. (2001a) Local interactions predict large-scale pattern in empirically
derived cellular automata. Nature, 413, 841–844.
Wootton, J.T. (2001b) Prediction in complex communities: analysis of empirically
derived Markov models. Ecology, 82, 580–598.
Yackulic, C.B., Reid, J., Davis, R., Hines, J.E., Nichols, J.D. & Forsman, E. (2012)
Neighborhood and habitat effects on vital rates: expansion of the Barred Owl in the
Oregon Coast Ranges. Ecology, 93, 1953–1966.
Yackulic, C.B., Reid, J., Nichols, J.D., Hines, J.E., Davis, R. & Forsman, E. (2014) The
roles of competition and habitat in the dynamics of populations and species
distributions. Ecology, 95, 265–279.
Supporting Information
Additional Supporting Informationmay be found online in the supporting information
tab for this article:
Appendix S1. Parameter estimation in the simulation study.
Appendix S2. Additional simulation study: anisotropic spatial correlation.
Appendix S3. Sampling design and model fitting to the rocky intertidal data.
Code S1. JAGS model code for the proposed model.
Code S2. An R script for the simulation study.
23
T
ab
le
1.
P
ar
am
et
er
es
ti
m
at
es
(t
ra
n
si
ti
on
p
ro
b
ab
il
it
y
m
at
ri
x
P
,
re
sa
m
p
li
n
g
er
ro
r
p
ro
b
ab
il
it
y
e,
an
d
b
an
d
w
id
th
m
at
ri
x
Σ
)
ob
ta
in
ed
b
y
fi
tt
in
g
ea
ch
of
th
e
th
re
e
d
iff
er
en
t
m
o
d
el
s
(n
ai
ve
es
ti
m
at
or
,
n
on
-s
p
at
ia
l
m
o
d
el
,
an
d
sp
at
ia
l
m
o
d
el
)
to
th
e
in
te
rt
id
al
se
ss
il
e
co
m
m
u
n
it
y
d
at
as
et
of
a
p
ar
ti
cu
la
r
q
u
ad
ra
t
(M
C
2)
.
F
or
th
e
sp
at
ia
l
an
d
n
on
-s
p
at
ia
l
m
o
d
el
s,
fo
r
w
h
ic
h
p
os
te
ri
or
sa
m
p
li
n
g
w
as
co
n
d
u
ct
ed
u
si
n
g
M
C
M
C
,
p
os
te
ri
or
m
o
d
es
(o
f
e
an
d
Σ
)
an
d
sp
at
ia
l
m
ed
ia
n
s
(o
f
P
)
ar
e
re
p
or
te
d
.
T
h
e
ro
w
s
an
d
co
lu
m
n
s
of
th
e
tr
an
si
ti
on
p
ro
b
ab
il
it
y
m
at
ri
x
co
rr
es
p
on
d
to
:
(1
)
fr
ee
sp
ac
e,
(2
)
th
e
b
ar
n
ac
le
C
ht
ha
m
al
u
s
da
ll
i,
(3
)
th
e
re
d
al
ga
G
lo
io
pe
lt
is
fu
rc
at
a,
(4
)
th
e
ar
ti
cu
la
te
d
ca
lc
ar
eo
u
s
re
d
al
ga
C
or
al
li
n
a
pi
lu
li
fe
ra
,
an
d
(5
)
ot
h
er
or
ga
n
is
m
s.
M
o
d
el
P
ar
am
et
er
N
a
iv
e
N
o
n
-s
p
a
ti
a
l
S
p
a
ti
al
P
     0.5
62
0
.3
3
0
0
.3
7
7
0
.0
4
3
0
.1
6
4
0.
0
85
0
.2
6
8
0
.1
0
0
0
.0
1
1
0
.0
8
0
0.
3
21
0
.3
4
4
0
.4
4
7
0
.0
1
1
0
.1
4
1
0.
0
05
0
.0
0
9
0
.0
0
1
0
.4
8
4
0
.1
5
5
0.
0
27
0
.0
4
9
0
.0
7
5
0
.4
5
1
0
.4
6
0     
     0.6
7
0
0
.2
2
1
0
.3
0
7
0
.0
1
3
0
.0
28
0
.0
6
4
0
.4
1
6
0
.1
0
0
0
.0
1
4
0
.0
1
6
0
.2
5
8
0
.3
2
8
0
.5
8
5
0
.0
1
1
0
.0
20
0
.0
0
3
0
.0
0
5
0
.0
0
2
0
.5
0
1
0
.1
9
6
0
.0
0
5
0
.0
2
9
0
.0
0
6
0
.4
6
1
0
.7
40
          
0.
77
2
0
.1
6
0
0
.2
08
0
.0
15
0
.0
40
0.
05
3
0
.5
7
4
0
.0
60
0
.0
16
0
.0
29
0.
16
5
0
.2
3
9
0
.7
13
0
.0
15
0
.0
37
0.
00
3
0
.0
0
8
0
.0
02
0
.6
84
0
.0
98
0.
00
7
0
.0
2
0
0
.0
17
0
.2
70
0
.7
96
     
e
0.
2
3
9
0.
71
5
Σ
d
ia
g(
0
.6
4
4,
1.
27
8)
24
(A) (B) 
Correct site location 
 (mitn = 0) 
Incorrect site location 
 (mitn = 1) 
yitn = zit 
Occurrence of 
resampling error 
(Equation 2) 
Conditional data 
distribution 
(Equation 1) 
Quadrat 
Sites (fixed points) 
Permanent marks 
yitn ~ Categorical(git1, ..., gitS) 
Observation 
Fig. 1. (A) A diagram of sampling scheme for sessile communities considered in this
study. The occupancy state of each site is identified in the field, or in a room if the
quadrat is photographed. The location of each site is not marked, but referred by using
some sampling device that indicate the position of the site relative to the permanent
marks which are established on the periphery of the quadrat. Resampling errors, failures
to find the correct location of the site, may occur when the sampling device is incorrectly
positioned and/or an observation is made with an inappropriate angle. (B) A hierarchical
representation of the observation process assumed in the proposed model.
25
Mean Square Error Square of Bias Variance
0.000
0.005
0.010
0.00
0.01
0.02
0.03
0.00
0.05
0.10
0.15
P
e
sigm
a
0.0 0.2 0.4 0.6 0.0 0.2 0.4 0.6 0.0 0.2 0.4 0.6
Error rate
V
al
ue
Model
Spatial
Non−spatial
Naive
Replicates
1
3
Fig. 2. Results of simulation study. Mean square error (left column), square of bias
(middle column), and variance (right column) of parameter estimates obtained from 96
replicated datasets using three different models (spatial model, non-spatial model, and
naive estimator) are shown. Upper, middle, and lower panels indicate results for the tran-
sition probability matrix P̂, resampling error rate ê and the bandwidth scale parameter
σ̂, respectively. Results of σ̂ for e = 0 in the spatial model are not shown because of a
convergence issue (see text for details).
26
2.5
5.0
7.5
10.0
12.5
1
2
3
4
5
M
ean Turnover T
im
e
D
am
ping R
atio
MC1 MC2 MC3 MC4 MC5 MB1 MB2 MB3 MB4 MB5 AP1 AP2 AP3 AP4 AP5 MZ1 MZ2 MZ3 MZ4 MZ5 NN1 NN2 NN3 NN4 NN5
Quadrat
V
al
ue
Model
Spatial
Non−spatial
Naive
Fig. 3. Mean turnover time (upper panel) and damping ratio (lower panel) estimated us-
ing three models (spatial model, non-spatial model, and naive estimator) for each quadrat.
27
Free space
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
Chthamalus dalli
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
Gloiopeltis fucata
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
Corallina pilulifera
●
●
● ●
●
●
●
●
●
●
●
●
●
●
Other organisms
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
Fig. 4. Observed data, occupancy, and local state dominance for a particular quadrat
(MC2) and point in time estimated using the spatial model, mapped on the spatial co-
ordinates of each observation point (20 horizontal rows × 10 vertical columns). Filled
circles (•) and crosses (+) represent data (the state observed at each observation point)
and estimated site occupancy state (the state that had the highest posterior probability
of occupancy at each observation point), respectively. Note that for each site, both filled
circles and crosses appear in only one of the five panels, which represent different ecolog-
ical states. The estimated local state dominance, gits, is represented by colored squares.
A brighter square indicates greater dominance of the state at that location.
28