CaptureRecaptureLab.knit

Motivation

We illustrate the PX-DA approach to implementing the closed-population capture-recapture model using data corresponding to encounter histories for two species of salamander in a 15 \(\times\) 15m fenced plot in Great Smoky Mountains National Park (GSMNP; Bailey et al., 2004). The majority of salamanders in GSMNP comprise two species: the red-cheeked salamander (Plethodon jordani) or the pygmy salamander (Desmognathus wrighti). It is thought that both species had fairly low detection probabilities so that researchers were unable to count all individuals, even within a truly closed population (in which closure involved a fenced area and short sampling period). Our goal is to estimate the abundance of both species.

R libraries and data needed for this lab

You need to load the following libraries and data. Set the seed to 10 to compare your answers to ours.

library(rjags)
library(MCMCvis)
library(HDInterval)
# library(BayesNSF)
sal.df <- read.table("sal_data.txt", header=TRUE)
set.seed(10)

Problem

The data frame sal.df consists of 4 columns (y1 through y4) indicating the encounter histories of 225 individuals. The fifth column spp indicates the species, where 0 represents the red-cheeked salamander and 1 represents the pygmy salamander. We will develop one model to be fit to either species separately.

Draw the DAG.

See whiteboard!

Write the mathematical expression for the joint posterior distribution.

\[ [p, \psi, \mathbf{z}|\mathbf{y}] \propto \prod_{i=1}^M([y_i|p, z_i]\times[z_i|\psi])\times[p]\times[\psi] \]

Implement the JAGS model in lecture. Augment the data with all zero encounter histories up to \(M=1500\) (a reasonable size of super-population). Use a uniform prior with support on \((0, 1)\) for the detection probability \(p\), and a beta prior with \(\alpha_\psi = 0.001\) and \(\beta_\psi = 1\) for the membership probability \(\psi\). Check chains for convergence.

y0 <- apply(sal.df[sal.df$spp==0, 1:4], 1, sum) # P. jordani
y1 <- apply(sal.df[sal.df$spp==1, 1:4], 1, sum) # D. wrighti

J <- 4
M <- 1500
n0 <- length(y0)
n1 <- length(y1)
y0.aug <- c(y0, rep(0, M-n0))
y1.aug <- c(y1, rep(0, M-n1))

m.jags <-"
  model{
    for(i in 1:M){
      y[i] ~ dbin(z[i]*p, J) 
      z[i] ~ dbern(psi)
    }
    p ~ dunif(0, 1)
    psi ~ dbeta(0.001, 1)
    N = sum(z)
}
"
n.mcmc <- 10000
n.burn <- round(.2*n.mcmc)

# P. jordani
mod<-textConnection(m.jags)
m0.out<-jags.model(mod, data=list('y'=y0.aug, 'J'=J, 'M'=M), inits=list(z=rep(1, M)), n.chains=1, n.adapt=0)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1500
##    Unobserved stochastic nodes: 1502
##    Total graph size: 4509
## 
## Initializing model

update(m0.out, n.mcmc)
m0.samples <- jags.samples(m0.out, c('p','psi','N'), n.mcmc)

## NOTE: Stopping adaptation

# D. wrighti
mod<-textConnection(m.jags)
m1.out<-jags.model(mod, data=list('y'=y1.aug, 'J'=J, 'M'=M), inits=list(z=rep(1, M)), n.chains=1, n.adapt=0)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1500
##    Unobserved stochastic nodes: 1502
##    Total graph size: 4509
## 
## Initializing model

update(m1.out, n.mcmc)
m1.samples <- jags.samples(m1.out, c('p','psi','N'), n.mcmc)

## NOTE: Stopping adaptation

Summarize the posterior detection probability and abundance for either species, respectively. Do detection probabilities differ by species? What is the posterior probability that there were more pygmy salamanders than red-cheeked salamanders?

layout(matrix(1:4, 2, 2))

hist(m0.samples$p[1,,1], prob=TRUE, ylab="Probability", xlab="Detection", main="P. jordani")

hist(m1.samples$p[1,,1], prob=TRUE, ylab="Probability", xlab="Detection", main="D. wrighti")

hist(m0.samples$N[1,,1], prob=TRUE, breaks=seq(0,M,1), xlim=c(0,M), ylab="Abundance", xlab="N", main="")

hist(m1.samples$N[1,,1], prob=TRUE, breaks=seq(0,M,1), xlim=c(0,M), ylab="Abundance", xlab="N", main="")

# posterior probability more pygmy salamanders
pp <- sum(m1.samples$N[1,,1]>m0.samples$N[1,,1])/n.mcmc

Modify the hyperprior parameters on \(\psi\) (e.g., increase \(\alpha_\psi\) while keeping \(\beta_\psi\) the same). How does it affect inference on abundance? Alternatively, modify \(M\) and explore how it affects posterior abundance.

Explore on your own.

Salamander Capture-Recapture Modeling Lab

June 14, 2023

Motivation

R libraries and data needed for this lab

Problem

References