# Fitting occupancy models with E-SURGE: Hidden Markov modelling of presence-absence data

Olivier Gimenez, Laetitia Blanc, Aurélien Besnard, Roger Pradel, Paul F. Doherty, Eric Marboutin, Rémi Choquet

► To cite this version:

Olivier Gimenez, Laetitia Blanc, Aurélien Besnard, Roger Pradel, Paul F. Doherty, et al. Fitting occupancy models with E-SURGE: Hidden Markov modelling of presence-absence data. Methods in Ecology and Evolution, 2014, 5 (6), pp.592-597. 10.1111/2041-210X.12191. hal-01911593

HAL Id: hal-01911593

https://hal.science/hal-01911593v1

Submitted on 17 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# FITTING OCCUPANCY MODELS WITH E-SURGE:
## HIDDEN MARKOV MODELLING OF PRESENCE-ABSENCE DATA

Olivier Gimenez¹, Laetitia Blanc¹, Aurélien Besnard¹, Roger Pradel¹, Paul F. Doherty Jr², Eric Marboutin³, and Rémi Choquet¹

¹Centre d'Ecologie Fonctionnelle et Evolutive, UMR 5175, campus CNRS, 1919 Route de Mende, 34293 Montpellier Cedex 5, France.

²Department of Fish, Wildlife and Conservation Biology, Colorado State University, Fort Collins, CO 80523-1474, USA

³Office National de la Chasse et de la Faune Sauvage, ZI Mayencin, 5 allée de Béthléem, 38610 Gières, France

### Abstract.

1. Occupancy – the proportion of area occupied by a species – is a key notion for addressing important questions in ecology, biogeography and conservation biology. Occupancy models allow estimating and inferring about species occurrence while accounting for false absences (or imperfect species detection).

This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process, which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1111/2041-210X.12191

This article is protected by copyright. All rights reserved.

Accepted Article
2. Occupancy models can be formulated as hidden Markov models (HMM) in which the state process captures the Markovian dynamic of the actual but latent states while the observation process consists of observations that are made from these underlying states.

3. We show how occupancy models can be implemented in program E-SURGE, which was initially developed to analyse capture-recapture data in the HMM framework. Replacing individuals by sites provides the user with access to several features of E-SURGE that are *not available altogether* or *just not available* in standard occupancy software: i) flexible model specification through a user-friendly syntax without having to write custom code, ii) decomposition of the observation and state processes in several steps to provide flexible parameterisation, iii) up-to-date diagnostics of model identifiability and iv) advanced numerical algorithms to produce fast and reliable results (including site random effects).

4. To illustrate E-SURGE features, we provide implementation and analysis details for several occupancy models. We also provide simulated and real-world examples as well as further specifications and information in a companion wiki platform http://occupancyinesurge.wikidot.com/.

*Key words:* capture-recapture; detectability; detection-nondetection; E-SURGE; hidden Markov models; presence-absence; species occurrence.

## INTRODUCTION

Occupancy models allow estimating and inferring about species occurrence while accounting for false absences or imperfect species detection (MacKenzie *et al.* 2006). These models have been extensively used to address important questions in fields as diverse as conservation biology, biogeography, wildlife epidemiology, metapopulation dynamics and
community ecology (review in Bailey, MacKenzie, & Nichols 2013).

Following the seminal work of MacKenzie and colleagues (MacKenzie *et al.* 2006), it was soon realized that occupancy models could be formulated as hidden Markov models (HMMs) in which two time series run in parallel: the state process captures the Markovian dynamic of the actual but latent states (e.g., site occupied vs. unoccupied) while the observation process consists of observations that are made from these underlying states (e.g., species detected vs. undetected; e.g., Royle & Kéry 2007). Interestingly, the formulation of capture-recapture models as HMMs was also witnessed in the capture-recapture literature (Pradel 2005; Gimenez *et al.* 2012). A close connection exists between capture-recapture and occupancy models that can be realized by interchanging individuals in a capture-recapture context and sites in an occupancy context.

Several software are available to fit occupancy models, either in the Frequentist framework with programs PRESENCE (Hines 2013), MARK (White & Burnham 1999) and the R package unmarked (Fiske & Chandler 2011), or in the Bayesian framework using the BUGS language (WinBUGS, OpenBUGS or JAGS; e.g., Kéry & Schaub 2012). WinBUGS requires writing custom code and knowledge about the Bayesian theory. Programs PRESENCE and MARK often require the construction of design matrices to specify models, a process that can be prone to user's errors. PRESENCE and unmarked do not incorporate random effects and MARK does so only in limited ways.

Here, by exploiting the equivalence between occupancy and capture-recapture models, we illustrate how occupancy models can be implemented in program E-SURGE (Choquet, Rouan, & Pradel 2009) which was initially developed to analyse capture-recapture data in the HMM framework. We believe that some users will benefit greatly by realizing that E-SURGE has features that are not implemented together in available occupancy software, namely i) flexible model specification through a user-friendly syntax without having to write custom

Accepted Article
code, ii) advanced numerical algorithms to produce fast and reliable results, including the incorporation of site random effects, and several other features that are simply *not implemented* in these software, namely iii) decomposition of the observation and state processes in several steps to provide flexible parameterisation and, iv) up-to-date diagnostics of model identifiability, in other words a reliable way of counting the number of parameters entering the calculation of the Akaike Information Criterion. We believe that E-SURGE is a powerful option in the suite of software that can be used for occupancy problems; however few users may recognize this option, even if they are current E-SURGE users. Our goal in this paper is to make E-SURGE more accessible to analysts and field biologists working on occupancy, and related, problems.

## HIDDEN MARKOV MODELLING OF OCCUPANCY DATA

In Figure 1, we provide the HMM formulation of the general dynamic occupancy models to carry out inference about occurrence and how extinction and colonization drive changes in occurrence.

The parameters of interest are the probability of local extinction $\varepsilon$ and of colonization $\gamma$ as well as the detection probability $p$ and the probability of initial occupancy $\psi_1$ where we have assumed, for simplicity, all parameters are constant across periods and sites. A HMM is built around three pieces of information: the vector of initial state probabilities, the matrix of transition probabilities linking states in successive sampling occasions, and the matrix of observation probabilities linking observations and states at a given occasion. At the first sampling occasion $t = 1$, with the first state being ‘unoccupied’ and the second ‘occupied’, the vector of initial state probabilities is:

Accepted Article
$$[1-\psi_1 \ \psi_1] \qquad (1)$$

Then, the states are distributed as a first-order Markov chain governed by the transition matrix with states unoccupied and occupied at t in rows and states unoccupied and occupied at t + 1 in columns:

$$\begin{bmatrix} 1-\gamma & \gamma \\ \varepsilon & 1-\varepsilon \end{bmatrix} \qquad (2)$$

The observation process conditional on underlying occupancy states is summarized by a matrix with unoccupied and occupied states at t (rows) and undetected and detected observations at visits j (columns):

$$\begin{bmatrix} 1 & 0 \\ 1-p & p \end{bmatrix} \qquad (3)$$

Single-season occupancy models can be reformulated as HMMs and fitted in E-SURGE by imposing no extinction ($\varepsilon=0$) and no colonisation ($\gamma=0$) in the dynamic model. The extension to multiple states with uncertainty (Nichols et al. 2007; MacKenzie et al. 2009) is illustrated with breeding states and we consider three possible states. A site is ‘unoccupied’, ‘occupied by non-breeders’ or ‘occupied by breeders’, while the observations are species undetected (coded 0), species detected without young (coded 1) and species detected with young (coded 2). We use $\psi^1$ (resp. $\psi^2$) the probability that the site is occupied by non-breeders (resp. by breeders), $p^1$ (resp. $p^2$) the detection probability of non-breeders (resp. of breeders). The model can also accommodate state uncertainty. For example,
reproduction may occur on a site, but the young may not be detected. We introduce $\delta$ the probability of detecting evidence of reproduction, given the site is occupied with young. Thus, the vector of initial state probabilities is:

$$[1 - \psi^1 - \psi^2 \quad \psi^1 \quad \psi^2] \qquad (4)$$

while the transition matrix is the identity matrix. The primary modifications are in the observation matrix which can be written as

$$\begin{bmatrix} 1 & 0 & 0 \\ 1 - p^1 & p^1 & 0 \\ 1 - p^2 & p^2(1 - \delta) & p^2\delta \end{bmatrix} \qquad (5)$$

which can be split in two steps to highlight the successive processes of detection and breeding state ascertainment:

$$\begin{bmatrix} 1 & 0 & 0 \\ 1 - p^1 & p^1 & 0 \\ 1 - p^2 & 0 & p^2 \end{bmatrix} \times \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 - \delta & \delta \end{bmatrix} \qquad (6)$$

We introduce an intermediate set of observations: undetected, detection of non-breeders, and detection of breeders. The first matrix in the product (6) specifies the probabilities of these intermediate observations conditional on the states, while the second matrix specifies the probabilities of reproduction conditional on the intermediate observations.
FITTING MODELS IN E-SURGE

A typical session in E-SURGE follows the workflow provided in Figure 2. Here, we
provide a brief description of the main features of E-SURGE. We refer to the E-SURGE
manual (Choquet & Nogué 2013) as well as Choquet (2008) and Choquet, Rouan, & Pradel
(2009) for additional details.

There are a few standard preliminary steps that we do not detail here. Most
importantly, the user needs to select the Occupancy option. The number of events to specify is
the number of observations (e.g., undetected, detection of nonbreeders and detection of
breeders), while the number of states should always be the number of states to be used (e.g.,
unoccupied, occupied by non-breeders, occupied by breeders) plus one. This additional state
is needed because E-SURGE was initially developed to estimate demographic parameters and
uses the absorbing state ‘dead’ to handle mortality as a transition from ‘alive’ states to the
‘dead’ state.

The model construction is served in E-SURGE by two modules called GEPAT
(GEnerator of PATtern of elementary matrices) and GEMACO (Generator of Matrices of
Constraints; Choquet 2008). In GEPAT, the user specifies the structure of the vector of initial
state probabilities, the matrix of transitions governing the state process and the matrix of
observations (conditional on the states) governing the observation process. To specify a
parameter that will be estimated (i.e., that will be assigned an effect in GEMACO, see below),
any letter can be used. The minus sign `-` means that the parameter corresponding to this cell
will always be set to 0 while the asterisk '*' means the complementary of the sum of all the
other parameters on the same row.

Of practical interest, the three elements of a HMM can be specified through a
multistep process that proves very useful in accommodating state uncertainty (e.g., eqn. 4).

<figure><img src="image_2.png" /></figure>
In GEMACO, we specify the effects (*sensu* the design matrix in programs MARK and PRESENCE) using a self-explanatory syntax: for example, a season effect will be specified by 't' for time, a group effect by 'g' while a constant effect will be 'i' for intercept. If the effect of a site covariate needs to be investigated, we use 'i + xind' where xind specifies the slope of the relationship. A site random effect can also be specified in the same spirit using the syntax 'i + ind' (see below for an example). The matrices defined at the GEPAT step can be manipulated using the syntax 'from' for rows and 'to' for columns. For example, f(1).to(2) will pick the element in row 1 and column 2 of the corresponding matrix, and if a time effect is required on this element, then the syntax will be f(1).to(2).t.

## CASE STUDIES IN E-SURGE USING SIMULATED DATA

To illustrate the use of E-SURGE for fitting occupancy models, we used simulated data.

Results from fitting the models to these data are provided in Table 1. Estimates for all parameters were close to the true values and the 95% confidence intervals covered the true values in all cases.

Below, we detail the most important steps of implementing these models in E-SURGE. We refer the reader to the companion wiki website http://occupancyinesurge.wikidot.com/ for full details on how to simulate the data and implement the model.

### *Dynamic models*

We start with the dynamic occupancy model described in formulas (1), (2) and (3) above.

This mathematical formulation of the model can be translated for E-SURGE as follows. In GEPAT, the vector of initial state probabilities is specified as

Accepted Article
* Ψ

which corresponds to (1). The additional dead state is not displayed at this step. Then, the
transition matrix in (2) is specified as

* γ -

ε * -

- - *

where the last row and column are for the dead state. Finally, the observation matrix in (3) is
coded as:

* -

* p

* -

where the last row corresponds to the dead state. In GEMACO, we write *i* at the Initial state and Event steps to impose a constant parameter. Regarding the Transition step, the data were simulated with 3 seasons (primary sessions) and 3 visits within each season (secondary session), therefore we use `to.t(1 2 4 5 7 8)+to.t(3 6)` (or `to.t(1 2 4 5 7 8,3 6)`).

The `to` specifies the columns of the matrix to be different, resulting in distinguishing the colonization and extinction parameters (note that `from` would produce exactly the same result by differentiating the rows). To handle the robust design, `t(1 2 4 5 7 8)` makes the intervals between secondary occasions identical, and the corresponding parameters will be fixed to 0 to impose closure within a primary session (i.e., neither extinction nor colonization). The term `t(3 6)` makes the intervals between the last visit in a season and the first visit of the next one (between primary sessions) identical.
Single-season models

We view a single-season model as a constrained dynamic model. Thus we start with the
dynamic model as described above. In GEPAT, the specification for the Initial state and the
Event steps is exactly the same as for the dynamic model. To satisfy the closure assumption,
we impose neither colonization nor extinction between sampling occasions by specifying at
the Transition step:

* - -
- * -
- - *

In GEMACO, we use i for all steps to obtain constant parameters.

The extension of this model to multiple states with uncertainty is obtained as follows
in E-SURGE. In GEPAT, the initial state probabilities in (4) are:

* Ψ Ψ

while the Transition step is the same as before. The difference lies in the specification of the
observation process that is accomplished in two steps to match (6) in GEPAT: step 1, the
detection matrix, is:

* - - -
* p - -
* - p -
* - - -

while step 2, the reproduction assignment matrix, is:
* - -
- * -
- * δ
* - -

In GEMACO, we use *to* for the Initial state step to distinguish the occupancy probabilities according to states, *from* for step 1 of the Event step to distinguish the detection probabilities according to states and *i* for step 2 to have a constant assignment probability.

A REAL-WORLD EXAMPLE ILLUSTRATING HETEROGENEITY IN SPECIES DETECTION

We use a real-world example to further illustrate E-SURGE capabilities. In particular, we show how E-SURGE can easily accommodate site-specific covariates, deal with heterogeneity in the species detection process, and identify parameter redundancy issues.

Based on signs of presence collected between 2002 and 2006 by a network of trained volunteers, Rolland *et al.* (2011) modelled the occurrence of the Eurasian lynx (*Lynx lynx*) in the French Jura mountains. Briefly, a grid with cells of 9km x 9km was used, resulting in 197 sites, and the study period was split into five 1-year periods. Observers searched for signs of lynx presence and each cell was assumed to be visited at least once each year.

The authors used a single-season occupancy model as described above. They, via model selection, found that the probability that a site was used by lynx was best explained by the proportion of forest cover $x_i$ in site *i*. Thus on the logit scale, we have:

$$ \operatorname{logit}(\psi_i) = \beta_0 + \beta_1 x_i \quad (7) $$

where the $\beta$'s are regression parameters to be estimated.

To further incorporate site heterogeneity in the species detection process, we
implemented a model with a site random effect and another one with a 2-class finite mixture.

A feature of E-SURGE is that it allows accommodating site random effects. Let us consider site-to-site variation in the detection probability as follows:

$$ \operatorname{logit}(p_i) = \mu + \eta_i \qquad (8) $$

where $\mu$ is the mean detection probability and $\eta_i$ are normally distributed with mean 0 and variance $\sigma^2$. Maximum likelihood estimates of $\mu$ and $\sigma^2$ are obtained through numerical integration via the Gauss-Hermite quadrature (Gimenez & Choquet 2010). More complex random structures can also be accommodated (Choquet & Gimenez 2012; Choquet et al. 2013).

As an alternative, we consider a 2-class finite mixture model (Royle 2006). We considered 4 states: unoccupied in class 1, unoccupied in class 2, occupied in class 1, occupied in class 2. There are two observations: species undetected (0) or detected (1). The vector of initial state probabilities is:

$$ \begin{bmatrix} \pi(1-\psi_1) & (1-\pi)(1-\psi_2) & \pi\psi_1 & (1-\pi)\psi_2 \end{bmatrix} \qquad (9) $$

where $\pi$ is the proportion of sites that belong to class 1 and $\psi_1$ (resp. $\psi_2$) is the occupancy probability for sites that belong to class 1 (resp. 2). We obtain $\pi$ by summing the first and third components in (9), while the ratio of the fourth component over $\pi$ gives $1-\psi_2$ hence $\psi_2$ and the ratio of the third component over $\pi$ gives $\psi_1$.

However, when this mixture model is fitted, E-SURGE detects that the model is parameter-redundant. More precisely, parameters $\pi$, $\psi_1$ and $\psi_2$ cannot be estimated
separately, and the number of estimable parameters is 4 while the model has 5. To make these parameters estimable, the constraint $\psi_1 = \psi_2$ can be imposed to ensure that the occupancy probability is no longer redundant. This means disentangling i) the assignment of sites to a class of heterogeneity and ii) the occupancy process. Here again, as in the multistate occupancy model with uncertainty in (6), the ability in E-SURGE to specify the state or the observation process in several steps proves useful. More precisely, we write the vector of initial state probabilities (9) as the matrix product of X and Y in GEPAT where:

$$X = \begin{bmatrix} \pi & 1-\pi \end{bmatrix} \qquad (10)$$

$$Y = \begin{bmatrix} 1-\psi_1 & 0 & \psi_1 & 0 \\ 0 & 1-\psi_2 & 0 & \psi_2 \end{bmatrix} \qquad (11)$$

then force the $\psi$'s to be the same in GEMACO. The matrix of observation probabilities is:

$$\begin{bmatrix} 1 & 0 \\ 1 & 0 \\ 1-p_1 & p_1 \\ 1-p_2 & p_2 \end{bmatrix} \qquad (12)$$

where the states are in rows and the observations in columns.

We fitted 3 models to the lynx data: $M_0 = [\psi(\text{forest}), p]$,
$M_{\text{ran}} = [\psi(\text{forest}), p(\text{random})]$ and $M_{\text{mix}} = [\psi(\text{forest}), p(\text{mixture})]$ with self-explanatory notation. The covariate forest cover had a similar positive effect on the occupancy probability whatever the structure on the detection probability

$(\beta_1 = 7.4 (\text{SE}=1.3) \text{ in } M_0, \beta_1 = 9.3 (\text{SE}=0.5) \text{ in } M_{\text{ran}} \text{ and } \beta_1 = 9.1 (\text{SE}=1.8) \text{ in } M_{\text{mix}})$. Site-to-
site variation in the detection process was equally captured by both models with heterogeneity (AICc[$M_0$]=1000.1 vs. AICc[$M_{ran}$]=978.6 and AICc[$M_{mix}$]=980.2). In $M_{ran}$, the mean detection probability was estimated at $\hat{\mu}$ = 0.10 (SE = 0.20) (hence $\text{logit}^{-1}(\hat{\mu})$ = 0.53 (SE = 0.05)) while the estimated standard deviation of the random effect was $\hat{\sigma}$ = 1.25 (SE = 0.22). In $M_{mix}$, the detection probability in state 1 was $\hat{p}_1$ = 0.25 (SE = 0.05) in proportion $\hat{\pi}$ = 0.54 (SE = 0.10) while the detection probability in state 2 was $\hat{p}_2$ = 0.74 (SE = 0.08). By having the ability to easily construct heterogeneity models for detection probability in conjunction with occupancy we have more confidence in our results. The relationship between lynx occupancy and forest cover is stronger with the highest ranked models ($M_{ran}$ and $M_{mix}$) than with a lower ranking model that does not include any additional structure on detection ($M_0$). By modelling heterogeneity in detection we are better able to uncover and describe important covariates of occupancy, which hopefully leads to better management.

## FURTHER E-SURGE CAPABILITIES

From a numerical perspective, E-SURGE provides a reliable number of parameters via an algorithm described in Choquet & Cole (2012), which is crucial in particular with parameter-redundant models or boundary estimates, and is one of the key steps for correct model selection using AIC. To date, we are aware of a single study investigating parameter redundancy in occupancy models which, actually, used E-SURGE (Fréville *et al.* 2013). With the ongoing development of complex models, the issue of parameter redundancy is likely to occur more often and this algorithm in E-SURGE should be useful.

From a modelling perspective, E-SURGE offers the possibility to fit multiple species models (MacKenzie, Bailey, & Nichols 2004; Waddle *et al.* 2010). Additionally, users may often be interested in deriving the probability that a site is occupied for sites without

This article is protected by copyright. All rights reserved.

Accepted Article
detections, either to build maps of presence or to calculate the so-called finite-population
occupancy (MacKenzie *et al.* 2006). To do so, the latent states of HMMs need to be estimated
and this can be achieved in a frequentist framework via the Viterbi algorithm (Rouan *et al.*
2009; Fiske *et al.* 2013), which is available in E-SURGE (Rouan *et al.* 2009).

The companion wiki website http://occupancyinesurge.wikidot.com/ presents such
examples.

CONCLUSIONS

Although initially developed for capture-recapture data, E-SURGE can be efficiently
used to build and analyse a variety of occupancy models via the HMM framework. E-SURGE
includes a user-friendly syntax for specifying models without having to write custom code
and uses advanced numerical algorithms to produce fast and reliable results. By strengthening
the link between the two fast growing user communities of capture-recapture and occupancy,
E-SURGE has the potential to provide a unified framework for the construction and analysis
of hidden-Markov models in ecology.

ACKNOWLEDGEMENTS

We thank the associate editor and three referees for comments on our work. We are indebted
to Marc Kéry for a careful and detailed review of our manuscript. L. Blanc was supported by
a PhD fellowship from Labex CeMeB

DATA ACCESSIBILITY

We provide the simulated data and the lynx data in the companion wiki platform
http://occupancyinesurge.wikidot.com/.

Accepted Article
Bailey, L.L., MacKenzie, D.I. & Nichols, J.D. (2013) Advances and applications of occupancy models. *Methods in Ecology and Evolution*.

Choquet, R. (2008) Automatic generation of multistate capture recapture models. *The Canadian Journal of Statistics*, **36**, 43–57.

Choquet, R. & Cole, D.J. (2012) A Hybrid Symbolic-Numerical Method for Determining Model Structure. *Mathematical Biosciences*, **236**, 117–125.

Choquet, R. and O. Gimenez (2012). Towards built-in capture-recapture mixed models in program E-SURGE. *Journal of Ornithology*, **152**, 625-639.

Choquet, R. & Nogué, E. (2013) E-SURGE 1.8 User's Manual. Montpellier.

Choquet, R., Sanz-Aguilar, A., Doligez, B., Nogué, E., Pradel, R., Gustafsson, L. & Gimenez, O. (2013) Estimating Demographic Parameters in the Wild with Dependence among Individuals within Cluster. Methods in Ecology and Evolution, 4, 474–482

Choquet, R., Rouan, L. & Pradel, R. (2009) Program E - SURGE : A Software Application for Fitting Multievent Models. *Environmental and Ecological Statistics* (eds D.L. Thomson, E.G. Cooch & M.J. Conroy), pp. 845–865. Springer US.

Fiske, I.J. & Chandler, R.B. (2011) unmarked : An R Package for Fitting Hierarchical Models of Wildlife Occurrence and Abundance. *Journal Of Statistical Software*, **43**, 1–23.

Fiske, I.J., Royle, J.A. & Gross, K. (2013) Inference for finite-sample trajectories in dynamic multi-state site-occupancy models using hidden Markov model smoothing.
*Environmental and Ecological Statistics*, DOI: 10.1007/s10651-013-0256-1.

Fréville, H., Choquet, R., Pradel, R. & Cheptou, P.-O. (2013). Inferring seed bank from hidden Markov models: new insights into metapopulation dynamics in plants. *Journal of Ecology*, **101**, 1572–1580.

Gimenez, O. & Choquet, R. (2010) Incorporating individual heterogeneity in studies on

<figure><img src="image_3.png" /></figure>
marked animals using numerical integration: capture-recapture mixed models. *Ecology*, **91**, 951–957.

Gimenez, O., Lebreton, J.-D., Gaillard, J.-M., Choquet, R. & Pradel, R. (2012) Estimating demographic parameters using hidden process dynamic models. *Theoretical Population Biology*, **82**, 307–316.

Hines, J.E. (2013) PRESENCE 5.9 – Software to estimate patch occupancy and related parameters.

Kéry, M. & Schaub, M. (2012) *Bayesian Population Analysis Using WinBUGS: A Hierarchical Perspective*. Academic Press.

MacKenzie, D.I., Bailey, L.L. & Nichols, J.D. (2004) Investigating species co-occurrence patterns when species are detected imperfectly. *Journal of Animal Ecology*, **73**, 546–555.

MacKenzie, D.I., Nichols, J.D., Royle, J.A., Pollock, K.H., Bailey, L.L. & Hines, J.E. (2006)
*Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurrence*. Academic Press.

MacKenzie, D.I., Nichols, J.D., Seamans, M.E. & Gutiérrez, R.J. (2009) Modeling species occurrence dynamics with multiple states and imperfect detection. *Ecology*, **90**, 823–835.

Nichols, J.D., Hines, A.J.E., Mackenzie, D.I., Seamans, M.E. & Gutiérrez, R.J. (2007)
Occupancy estimation and modeling with multiple states and state uncertainty. *Ecology*,
**88**, 1395–1400.

Pollock, K.H. (1982) A capture-recapture design robust to unequal probability of capture. *The Journal of Wildlife Management*, **46**, 752–757.

Pradel, R. (2005) Multievent: an extension of multistate capture-recapture models to uncertain states. *Biometrics*, **61**, 442–7.

Rolland, J., Basille, M., Marboutin, E. & Gaillard, J.-M. (2011) Comparing profile methods and site-occupancy modelling for the study of occurrence of an elusive species.

Accepted Article
Rouan, L., Gaillard, J.-M., Guédon, Y., Pradel, R., 2009b. Estimation of lifetime reproductive success when reproductive status cannot always be assessed. In: Thomson, D.L., Cooch, E.G., Conroy, M.J. (Eds.), *Modeling Demographic Processes in Marked Populations*. Springer, Berlin, Germany, pp.867–879.

Royle, J.A. (2006) Site occupancy models with heterogeneous detection probabilities. *Biometrics*, **62**, 97–102.

Royle, J.A. & Kéry, M. (2007) A Bayesian state-space formulation of dynamic occupancy models. *Ecology*, **88**, 1813–23.

White, G.C. & Burnham, K.P. (1999) Program MARK: survival estimation from populations of marked animals. *Bird Study*, **46**, 120–139.

<table><caption>**Table 1.** Estimates of the parameters in the occupancy models fitted to the simulated data set in E-SURGE. A single dataset was generated from the model under study with 250 sites, 3 primary occasions and 3 secondary occasions whenever required. For each parameter, the true value, the maximum likelihood estimate and the 95% confidence interval are provided.</caption><thead><tr><th>Model</th><th>Parameter</th><th>True value</th><th>Estimate</th><th>95% confidence interval</th></tr></thead><tbody><tr><td rowspan="4">Dynamic model</td><td>ψ<sub>1</sub></td><td>0.6</td><td>0.59</td><td>(0.53, 0.65)</td></tr><tr><td>γ</td><td>0.3</td><td>0.32</td><td>(0.27, 0.39)</td></tr><tr><td>ε</td><td>0.5</td><td>0.53</td><td>(0.47, 0.59)</td></tr><tr><td>p</td><td>0.7</td><td>0.70</td><td>(0.67, 0.73)</td></tr><tr><td rowspan="2">Single-season</td><td>ψ</td><td>0.8</td><td>0.8</td><td>(0.70, 0.87)</td></tr><tr><td>p</td><td>0.5</td><td>0.5</td><td>(0.44, 0.56)</td></tr></tbody></table>
<table><tr><td rowspan="4">Multistate with uncertainty</td><td>ψ<sup>1</sup></td><td>0.3</td><td>0.30</td><td>(0.23, 0.37)</td></tr><tr><td>ψ<sup>2</sup></td><td>0.5</td><td>0.50</td><td>(0.43, 0.57)</td></tr><tr><td>p<sup>1</sup></td><td>0.5</td><td>0.52</td><td>(0.43, 0.61)</td></tr><tr><td>p<sup>2</sup></td><td>0.7</td><td>0.70</td><td>(0.65, 0.76)</td></tr><tr><td></td><td>δ</td><td>0.8</td><td>0.79</td><td>(0.73, 0.85)</td></tr></table>

Figures legends

**Figure 1.** Schematic hidden Markov representation of a dynamic occupancy model with imperfect detection. We consider site *i* between seasons (primary periods) *t* and *t* + 1 and three visits (secondary periods) within each primary period. Each site is closed (i.e., no change in its occupancy status) within primary periods but open (i.e., allowing for changes in occupancy status) between primary periods (robust design sensu; Pollock 1982). The first layer (circles) is a succession of hidden states or latent states (0 for unoccupied vs. 1 for occupied) of site *i* at season *t* + 1 (*X*<sub>*it+1*</sub>) depending on its states at time *t* (*X*<sub>*it*</sub>). The dynamic of the states is driven by transition probabilities, here probabilities of colonization

γ = Pr(*X*<sub>*it+1*</sub> = 1 | *X*<sub>*it*</sub> = 0) and local extinction ε = Pr(*X*<sub>*it+1*</sub> = 0 | *X*<sub>*it*</sub> = 1). The second layer (squares) corresponds to the detection (*Y*<sub>*ij*</sub> = 1) or not (*Y*<sub>*ij*</sub> = 0) of the target species on site *i* at visit *j* = 1, 2 or 3 conditional on site *i* being in state *X*<sub>*it*</sub>. These events are driven by the species detection probability *p* = Pr(*Y*<sub>*ij*</sub> = 1 | *X*<sub>*it*</sub> = 1). The initial state probabilities ψ<sub>1</sub> = Pr(*X*<sub>*i1*</sub> = 1) and 1 − ψ<sub>1</sub> are not represented.
**Figure 2.** Workflow diagram for E-SURGE. We describe the successive steps of a typical analysis in E-SURGE, from data input to model fitting through model building and effects specification. Steps that need to be accomplished through pull-down menus are in white boxes; other steps can be done directly from the main interface. We provide details on key steps in using E-SURGE in the text or boxes.

<figure><img src="image_4.png" /></figure>

Accepted Article
Workflow for E-SURGE

<figure><img src="image_5.png" /></figure>