**1. Simulate data**

We modified the R codes provided in Kéry and Schaub (2012) to simulate detection/non-detection data from a dynamic occupancy model with constant parameters except extinction which was written as a function of a site-specific covariate. More precisely, we used: R = 250 sites, J = 3 visits within K = 3 primary sessions or years or seasons, initial occupancy psi1 = 0.6, detection p = 0.7, colonization gamma = 0.3 and extinction logit(epsilon)=0.1 + 0.5 * covariate.

```
# R – Number of sites
# J – Number of replicate surveys
# K – Number of years
# psi1 – occupancy probability in first year
# p – detection
# epsilon and gamma – extinction and colonization probabilities
R = 250
J = 3
K = 3
psi1 = 0.6
p = 0.7
gamma = 0.3
# Relationship survival – covariate
X <- runif(n = R, min = 0, max = 1)
epsilon <- plogis(0.1 + 0.5 * X) # Apply inverse logit
mean(epsilon)
# 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[i]) + (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)
}
}
}
# format data
yy <- matrix(y, R, J*K)
yy <- cbind(yy,X)
```

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

- 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,
*1 individual covariate*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*`* -`

`* p`

`* -`

- 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*`f.primary+[f(2).primary]>xind+secondary`where`primary`and`secondary`are shortcuts for`[t(3 6)]`and`[t(1 2 4 5 7 8)]`respectively; the operator larger than '>' allows to restrict the effect of a covariate to a set of occasions, groups, etc… For more details, see section 4.5 of the E-SURGE manual.*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#3 (colonization and extinction 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 (on the logit scale), the limits of the 95% confidence interval and the SE:

Beta# 1# | 0.420051128 -0.685237895 -0.154864361 +0.135299371 **initial occupancy**

#####################################################################

Beta# 2# | -0.932845751 -1.210953873 -0.654737629 +0.141891899 **colonization**

Beta# 3# | +0.233789021 -0.330620003 +0.798198044 +0.287963788 **intercept extinction**

Beta# 4# | +0.580983611 -0.365449389 +1.527416610 +0.482873979 **slope extinction**

#####################################################################

Beta# 5# | +0.743071199 +0.588339955 +0.897802443 +0.078944512 **detection**

#####################################################################

**4. Validation in Unmarked and Presence**

- with UNMARKED

Call:

colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~X,

pformula = ~1, data = countUMF)

Initial:

Estimate SE z P(>|z|)

0.42 0.135 3.1 0.00191

Colonization:

Estimate SE z P(>|z|)

-0.933 0.142 -6.57 4.89e-11

Extinction:

Estimate SE z P(>|z|)

(Intercept) 0.234 0.288 0.812 0.417

X 0.581 0.483 1.203 0.229

Detection:

Estimate SE z P(>|z|)

0.743 0.0789 9.41 4.84e-21

AIC: 2052.154

- with PRESENCE

Untransformed Estimates of coefficients for covariates (Beta's)

estimate std.error

A1 psi1 : 0.420051 0.135302

B1 gam1 : -0.932846 0.141895

C1 eps1 : 0.233789 0.287903

C2 eps1.C2 : 0.580983 0.482787

D1 P[1-1] : 0.743071 0.078943

**5. Further models**

In what precedes, site-specific covariates correspond to individual covariates in the capture-recapture terminology. If a covariate is measured at the year / season level (field effort for instance), then it can be incorporated as an external covariate using a text file that gathers the covariate values (see sections 4.4 and 5.4 of the E-SURGE user's manual) and is required while writing the syntax in GEMACO.