**1. Simulate data**

We use R code provided in Kéry and Schaub (2012) to simulate detection - nondetection data with psi = 0.8, p = 0.5, R = 200 sites and T = 3 seasons (or years).

```
# Sample sizes (spatial and temporal replication)
R <- 200
T <- 3
# Fix process parameters
psi <- 0.8 # Occupancy probability
p <- 0.5 # Detection probability
# Create structure to contain counts
y <- matrix(NA, nrow = R, ncol = T)
# Ecological process: Sample true occurrence (z, yes/no) from a Bernoulli (occurrence probability = psi)
z <- rbinom(n = R, size = 1, prob = psi) # Latent occurrence state
# Observation process: Sample detection/nondetection observations from a Bernoulli(with p) if z=1
for (j in 1:T){
y[,j] <- rbinom(n = R, size = 1, prob = z * p)
}
```

**2. Create the dataset to be analysed**

Copy and paste y created at the previous step in R in a text file. Then, format it as an .inp file.

The simulated dataset we will use in the analyses below can be found here.

**3. Analysis in E-SURGE**

We considered 2 states: unoccupied and occupied plus the the state dead in E-SURGE. The observations are species undetected (0) or detected (1).

The vector of initial state probabilities is (the dead state is not displayed):

[1-pi pi]

The transition matrix is:

[1 0 0

0 1 0

0 0 1]

The observation (or event) matrix is:

[1 0

1-p p

1 0]

- Start » New session (preferably in the directory where your dataset is)
- Data » Load data (Mark); click OK in the window that pops up asking 'How many columns do we extract from the data?'
- In the 'DATA' section in the main window, click the 'Modify' button and use 3 states and 1 age class
- Models » Markovian states only » Occupancy
- In the 'Advanced Numerical' section in the main window, tick the 'Compute C-I (Hessian)' box to get confidence intervals
- In the 'COMPUTE A MODEL' section in the main window, click on the 'Gepat' yellow button and use the following Matrix Patterns:

*Initial State*`* p`*Transition*`* - -`

`- * -`

`- - *`*Event*`* -`

`* b`

`* -`

- Click Exit to go to the next step.
- In the 'COMPUTE A MODEL' section in the main window, click on the 'Gemaco' green button and use the following syntax in the 'Model definition' dialog box (in the right bottom corner):

*Init state*`i`*Transition**Event*`i`

- Gemaco » Call Gemaco (all phrases) or Ctr+G, then click Exit
- In the 'COMPUTE A MODEL' section in the main window, click on the 'IVFV' pink button. Exit.
- In the 'COMPUTE A MODEL' section in the main window, click on the 'RUN' red button to fit the model to the simulated dataset.
- When the dialog box pops up, modify the model name if needed, then click OK
- In the 'Output' section of the main window, click on 'Selected Model Results (.out)' to get the results. More precisely, check out the 'Reduced set of parameters' section in the output file. The three lines below are organised as follows: initial state, transition and event parameters, with the maximum likelihood estimates, the limits of the 95% confidence interval and the SE:

Par# 1# IS( 1, 2)( 1, 1)( 1 1) | **0.799997178** 0.704026332 0.870570024 0.042427507 **(occupancy)**

Par# 18# E( 2, 2)( 1, 1)( 1 1) | **0.499997662** 0.441099929 0.558895459 0.030190038 **(detection)**

**4. Check results with program MARK**

Real Function Parameters of {p(.) Psi(.) PIM}

95% Confidence Interval

Parameter Estimate Standard Error Lower Upper

------ --— ---- --— ----

1:p **0.5000000** 0.0301904 0.4411016 0.5588984 **(detection)**

2:Psi **0.8000000** 0.0424264 0.7040317 0.8705711 **(occupancy)**