Simplest Dynamic Model

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.

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License