Site Level Covariate On Extinction

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.

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