**1. Simulate data**

We use a modified version of the R code provided in Kéry and Schaub (2012) to simulate detection/non-detection data from a dynamic occupancy model with constant parameters. These parameters are: R = 250 sites, J = 3 visits within K = 3 primary sessions or years or seasons, the initial occupancy is psi1 = 0.6, detection is p = 0.7, extinction is epsilon = 0.5 and colonization is gamma = 0.3. Recall that the occupancy status remains the same within primary periods (closure assumption) but may change between primary occasions.

```
# R – Number of sites
# J – Number of replicate surveys or visits
# K – Number of years or seasons
# psi1 – occupancy probability in first year / season
# p – detection
# epsilon and gamma – extinction and colonization probabilities
R = 250
J = 3
K = 3
psi1 = 0.6
p = 0.7
epsilon = 0.5
gamma = 0.3
# Set up some required arrays
site <- 1:R # Sites
year <- 1:K # Years
psi <- rep(NA, K) # Occupancy probability
muZ <- z <- array(dim = c(R, K)) # Expected and realized occurrence
y <- array(NA, dim = c(R, J, K)) # Detection histories
# Determine initial occupancy and demographic parameters
# Generate latent states of occurrence
# First year
z[,1] <- rbinom(R, 1, psi1) # Initial occupancy state
# Later years
for(i in 1:R){ # Loop over sites
for(k in 2:K){ # Loop over years
muZ[k] <- z[i, k-1]*(1-epsilon) + (1-z[i, k-1])*gamma # Prob for occ.
z[i,k] <- rbinom(1, 1, muZ[k])
}
}
# Generate detection/nondetection data
for(i in 1:R){
for(k in 1:K){
prob <- z[i,k] * p
for(j in 1:J){
y[i,j,k] <- rbinom(1, 1, prob)
}
}
}
# Compute annual population occupancy
for (k in 2:K){
psi[k] <- psi[k-1]*(1-epsilon) + (1-psi[k-1])*gamma
}
# format data
yy <- matrix(y, R, J*K)
```

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

Copy and paste yy 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**

The states are unoccupied, occupied and dead for E-SURGE while the observation are simply nondetected (0) or detected (1) for the species.

The vector of initial state probabilities is [1-psi1 psi1] (the dead state is not displayed in E-SURGE).

The transition matrix is:

[1-gamma gamma 0

epsilon 1-epsilon 0

0 0 1].

The observation matrix is:

[1 0

1-p p

1 0].

- Start » New session
- 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*`* y -`

`y * -`

`- - *`*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:

*Init state*`i`or`to`equivalently*Transition*`to.t(1 2 4 5 7 8)+to.t(3 6)`*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 and fix at the Transition step Beta#1 (local extinction within primary periods) and Beta#2 (colonization within primary periods) to 0 by modifying the 0.5 value and ticking the box in front of it. 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# 2# IS( 1, 2)( 1, 1)( 1 1) | 0.591558944 0.527325918 0.652810183 0.032174630 **Initial occupancy**

Par# 30# T( 2, 1)( 3, 1)( 1 1) | 0.530782337 0.466147058 0.594401142 0.032898369 **Extinction**

Par# 31# T( 1, 2)( 3, 1)( 1 1) | 0.324049015 0.267153279 0.386668564 0.030608580 **Colonization**

Par# 62# E( 2, 2)( 1, 1)( 1 1) | 0.702004402 0.670170072 0.731994902 0.015786328 **Detection**

**4. Check results with R package Unmarked**

Backtransformed linear combination(s) of **Initial** estimate(s)

Estimate SE LinComb (Intercept)

0.592 0.0322 0.37 1

Backtransformed linear combination(s) of **Colonization** estimate(s)

Estimate SE LinComb (Intercept)

0.324 0.0306 -0.735 1

Backtransformed linear combination(s) of **Extinction** estimate(s)

Estimate SE LinComb (Intercept)

0.531 0.0329 0.123 1

Backtransformed linear combination(s) of **Detection** estimate(s)

Estimate SE LinComb (Intercept)

0.702 0.0158 0.857 1

**5. Further models in E-SURGE**

To have time-dependent colonization and extinction probabilities, modify the syntax in GEMACO at the Transition step by distinguishing the 3rd and 6th time intervals with a comma in the second term `to.t(1 2 4 5 7 8)+to.t(3,6)` and fix the same values as before to 0. To get survey-specific detection probabilities, modify the syntax in GEMACO at the Event step by putting together the corresponding parameters `t(1)&t(4)&t(7)+t(2)&t(5)&t(8)+t(3)&t(6)&t(9)`. Shortcuts can be useful to specify such models.