Modelling Population Dynamics: Chapter 7

Survival estimation from mark-re-encounter data

Section 7.2.1 Capture-recapture data

WinBugs implementation of Cormack-Jolly-Seber model

model
{
for (i in 1:n)  # loop over individuals
    {

    # State-space model likelihood
    px[i,f[i]] <- 1
    z[i, f[i]] ~ dbern(px[i,f[i]])  # initial state is alive (= 1)

    for (j in (f[i]+1):K)  # loop over time (after first capture)
        {
        # state probabilities 'px' (= 0 if dead at t-1, = phi if alive at t-1)
        px[i,j] <- phi[j-1] * z[i,j-1]
        # state equation (sample state at t given state at t-1)
        z[i,j] ~ dbern(px[i,j])

        # observation probabilities 'po' (= 0 if dead at t, = p if alive at t)
        po[i,j] <- p[j-1] * z[i, j]
        # observation equation (sample observation at t given state at t)
        y[i, j] ~ dbern(po[i,j])
        }
    }

    # Priors
    for (j in 1:(T-1))
        {
        phi[j] ~ dunif(0, 1)  # survival probability
        p[j] ~ dunif(0,1)  # resighting probability
        }
}

Section 7.2.4 Capture-recapture of European Dipper data

Data themselves reside in file named ed.dat

WinBUGS model (named “phitpt.bug” in subsequent R code)

# State-space modelling of mark-recapture data
model

{

  for (i in 1:n)  # loop / individual

  {

    px[i,f[i]] <- 1
    z[i, f[i]] ~ dbern(px[i,f[i]])  # initial state is alive (= 1)

    for (j in (f[i]+1):K)  # loop / individual (after first capture)

    {

      # state probabilities 'px' (= 0 if dead at t-1, = phi if alive at t-1)
      px[i,j] <- phi[j-1] * z[i,j-1] 
      
      # state equation (sample state at t given state at t-1)
      z[i,j] ~ dbern(px[i,j])
      
      # observation probabilities 'po' (= 0 if dead at t, = p if alive at t)
      po[i,j] <- p[j-1] * z[i, j]
      
      # observation equation (sample observation at t given state at t)
      y[i, j] ~ dbern(po[i,j])

    }

  }

# Priors
for (j in 1:(K-1))  
{
phi[j] ~ dunif(0, 1)  # survival probability
p[j] ~ dunif(0,1)  # resighting probability
}
}

Finally the R code to calling R2WinBUGS

niter = 500
nburn = 50

#------------------------------#
#--- model (phi-t,p-t) --------#
#------------------------------#

# read encounter histories
setwd("appropriate for your local configuration")
data <- as.matrix(read.table("ed.dat")) # file of encounter histories

# number of individuals 
n <- dim(data)[[1]] 

# number of capture events
K <- dim(data)[[2]] 

# compute the date of first capture
fc <- NULL
for (i in 1:n)
{
temp <- 1:K
fc <- c(fc,min(temp[data[i,]==1]))
}

# data
datax <- list(n=n,K=K,y=data,f=fc)

# list of inits for states
dataIV = data
for (i in 1:n)
    {
    if (fc[i] > 1)
        {
       for (j in 1:(fc[i]-1))
	   {dataIV[i,j] = NA}
	}
     }

# first list of inits
init1 <- list(p=rep(0.1,6),phi=rep(0.1,6),z=as.matrix(dataIV))

# second list of inits
init2 <- list(p=rep(0.9,6),phi=rep(0.9,6),z=as.matrix(dataIV))

# concatenate initial values
inits <- list(init1,init2)

# parameters to be monitored
parameters <- c("phi","p")

# do the MCMC stuff calling WinBUGS
library(R2WinBUGS)
deb = Sys.time()
sim.phitpt <- bugs(datax, inits, parameters,"phitpt.bug", n.chains=2, n.iter=niter,n.burnin = nburn)
fin = Sys.time()
duree = fin-deb

save(sim.phitpt,duree,file = "CJSmodel.Rdata")
round(sim.phitpt$summary,2)