Survival, recruitment and population size of Darwin’s leaf-eared mouse (Phyllotis darwini)
Capture-Recapture Data, using the Jolly–Seber (JS) model
Paola Correa, Diego Lizcano, Salvador Mandujano, Grabriel Andrade-Ponce et al.
05/05/2022
Model from the book “Bayesian Population Analysis using WinBUGS — A Hierarchical Perspective” (2012) by Marc Kéry and Michael Schaub.
10.1 Intro to JS
Capture–recapture data contain information not only about the disappearance of marked individuals from a study population (i.e., mortality and permanent emigration) but also about their arrival into the population (i.e., recruitment by locally born individuals and immigration). Estimation of these recruitment-related parameters is possible when the complete capture histories of the marked individuals are analyzed, not just the part following first capture, as for the CJS model (Chapter 7).
All open-population capture–recapture types of models described in the Chapters 7–9 condition on first capture. This means that only the part of a capture-history after first capture is modeled; the information about entry of individuals is discarded. In this chapter, we introduce the Jolly–Seber (JS) model (Jolly, 1965; Seber, 1965), which allows estimation of gains to and losses from the population by not conditioning on first capture. Additionally, population size and the number of recruits per occasion can be estimated. The JS model is an extension of the closed-population models in Chapter 6, which do not condition on first capture either but assume demographic closure (i.e., the absence of mortality, recruitment, or dispersal), and of the CJS model.
Capture-recapture models rely on assumptions:
Design: No mark lost, Identity of individuals recorded without error (no false positives), Captured individuals are a random sample.
Model: Homogeneity of survival and recapture probabilities, Independence between individuals (overdispersion).
Since the JS model uses the complete capture-history, additional assumptions are necessary. All assumptions listed for the CJS model (Section 7.1) also apply to the JS model. In addition, the JS model requires that all individuals (marked and unmarked) in the population have the same capture probability. In other words, newly captured individuals must be a random sample of all unmarked individuals in a population (Williams et al., 2002). This assumption is likely to be met when the same sampling protocol is applied for capture and recapture, but not when protocols differ.
What does survival actually mean in capture-recapture ?
Survival refers to the study area.
Mortality and permanent emigration are confounded.
Therefore we estimate apparent survival, not true survival.
Apparent survival probability = true survival × study area fidelity.
Consequently, apparent survival < true survival unless study area fidelity = 1.
Use caution with interpretation. If possible, combine with ring-recovery data, or go spatial to get closer to true survival.
10.7. Analysis of a real data set: Survival, Recruitment and Population size
We are interested in population size and in the variability of local recruitment over time.
Initial modeling suggested that adult survival was subject to strong temporal variation, whereas recapture probability was constant over time. We therefore fit a model with temporal random effects in survival, fixed time effects in recruitment, and a constant capture rate. We denote this model as (ϕt, bt, p)
Jolly-Seber models give population size and the number of recruits per occasion.
Load data and run jags model
library(ggplot2)
library(readxl)
library(AHMbook)
library(jagsUI)
library(MCMCvis)
Pdarwini_Males <- read_excel("D:/BoxFiles/Box Sync/CodigoR/Chinchilla/data/Pdarwini_Male_history_fixed.xls")
Pdarwini_Females <- read_excel("D:/BoxFiles/Box Sync/CodigoR/Chinchilla/data/Pdarwini_Female_history_fixed.xls")
# 1st model using females only didnt work
Pdarwini <- rbind(Pdarwini_Males, Pdarwini_Females)
y1 <- as.matrix( Pdarwini[,1:34])
n_occasions1 <- 34
g1 <- 2
group1 <- as.integer(Pdarwini$group)
nind1 <- 3069
df1 <- 3
##################
cat("
model {
# Priors and constraints
for (i in 1:M){
for (t in 1:(n.occasions-1)){
logit(phi[i,t]) <- mean.lphi + epsilon[t]
} #t
for (t in 1:n.occasions){
p[i,t] <- mean.p
} #t
} #i
mean.p ~ dunif(0.3, 0.9) # Prior for mean capture
mean.phi ~ dunif(0.1, 0.5) # Prior for mean survival
mean.lphi <- log(mean.phi / (1-mean.phi))
for (t in 1:(n.occasions-1)){
epsilon[t] ~ dnorm(0, tau)
}
tau <- pow(sigma, -2)
sigma ~ dunif(0, 0.6) # Prior for sd of indv. variation of phi
sigma2 <- pow(sigma, 2)
for (t in 1:n.occasions){
gamma[t] ~ dunif(0, 1)
} #t
# Likelihood
for (i in 1:M){
# First occasion
# State process
z[i,1] ~ dbern(gamma[1])
mu1[i] <- z[i,1] * p[i,1]
# Observation process
y[i,1] ~ dbern(mu1[i])
# Subsequent occasions
for (t in 2:n.occasions){
# State process
q[i,t-1] <- 1-z[i,t-1]
mu2[i,t] <- phi[i,t-1] * z[i,t-1] + gamma[t] * prod(q[i,1:(t-1)])
z[i,t] ~ dbern(mu2[i,t])
# Observation process
mu3[i,t] <- z[i,t] * p[i,t]
y[i,t] ~ dbern(mu3[i,t])
} #t
} #i
# Calculate derived population parameters
for (t in 1:n.occasions){
qgamma[t] <- 1-gamma[t]
}
cprob[1] <- gamma[1]
for (t in 2:n.occasions){
cprob[t] <- gamma[t] * prod(qgamma[1:(t-1)])
} #t
psi <- sum(cprob[]) # Inclusion probability
for (t in 1:n.occasions){
b[t] <- cprob[t] / psi # Entry probability
} #t
for (i in 1:M){
recruit[i,1] <- z[i,1]
for (t in 2:n.occasions){
recruit[i,t] <- (1-z[i,t-1]) * z[i,t]
} #t
} #i
for (t in 1:n.occasions){
N[t] <- sum(z[1:M,t]) # Actual population size
B[t] <- sum(recruit[1:M,t]) # Number of entries
} #t
for (i in 1:M){
Nind[i] <- sum(z[i,1:n.occasions])
Nalive[i] <- 1-equals(Nind[i], 0)
} #i
Nsuper <- sum(Nalive[]) # Size of superpopulation
}",fill=TRUE)##
## model {
## # Priors and constraints
## for (i in 1:M){
## for (t in 1:(n.occasions-1)){
## logit(phi[i,t]) <- mean.lphi + epsilon[t]
## } #t
## for (t in 1:n.occasions){
## p[i,t] <- mean.p
## } #t
## } #i
## mean.p ~ dunif(0.3, 0.9) # Prior for mean capture
## mean.phi ~ dunif(0.1, 0.5) # Prior for mean survival
## mean.lphi <- log(mean.phi / (1-mean.phi))
## for (t in 1:(n.occasions-1)){
## epsilon[t] ~ dnorm(0, tau)
## }
## tau <- pow(sigma, -2)
## sigma ~ dunif(0, 0.6) # Prior for sd of indv. variation of phi
## sigma2 <- pow(sigma, 2)
##
## for (t in 1:n.occasions){
## gamma[t] ~ dunif(0, 1)
## } #t
##
## # Likelihood
## for (i in 1:M){
## # First occasion
## # State process
## z[i,1] ~ dbern(gamma[1])
## mu1[i] <- z[i,1] * p[i,1]
## # Observation process
## y[i,1] ~ dbern(mu1[i])
##
## # Subsequent occasions
## for (t in 2:n.occasions){
## # State process
## q[i,t-1] <- 1-z[i,t-1]
## mu2[i,t] <- phi[i,t-1] * z[i,t-1] + gamma[t] * prod(q[i,1:(t-1)])
## z[i,t] ~ dbern(mu2[i,t])
## # Observation process
## mu3[i,t] <- z[i,t] * p[i,t]
## y[i,t] ~ dbern(mu3[i,t])
## } #t
## } #i
##
## # Calculate derived population parameters
## for (t in 1:n.occasions){
## qgamma[t] <- 1-gamma[t]
## }
## cprob[1] <- gamma[1]
## for (t in 2:n.occasions){
## cprob[t] <- gamma[t] * prod(qgamma[1:(t-1)])
## } #t
## psi <- sum(cprob[]) # Inclusion probability
## for (t in 1:n.occasions){
## b[t] <- cprob[t] / psi # Entry probability
## } #t
## for (i in 1:M){
## recruit[i,1] <- z[i,1]
## for (t in 2:n.occasions){
## recruit[i,t] <- (1-z[i,t-1]) * z[i,t]
## } #t
## } #i
## for (t in 1:n.occasions){
## N[t] <- sum(z[1:M,t]) # Actual population size
## B[t] <- sum(recruit[1:M,t]) # Number of entries
## } #t
## for (i in 1:M){
## Nind[i] <- sum(z[i,1:n.occasions])
## Nalive[i] <- 1-equals(Nind[i], 0)
## } #i
## Nsuper <- sum(Nalive[]) # Size of superpopulation
## }leis <- y1
nz <- 500
CH.aug <- rbind(leis, matrix(0, ncol = dim(leis)[2], nrow = nz))
# Bundle data
jags.data <- list(y = CH.aug, n.occasions = dim(CH.aug)[2], M = dim(CH.aug)[1])
# Initial values
# Good initial values for the latent state z are needed to run the model in JAGS. The simplest option that works is to give just a matrix with a 1 at all places.
z.init <- CH.aug
z.init[z.init==0] <- 1
inits <- function(){list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), sigma = runif(1, 0, 1), z = z.init)}
# Parameters monitored
parameters <- c("psi", "mean.p", "sigma2", "mean.phi", "N", "Nsuper", "b", "B")
# MCMC settings
ni <- 8000
nt <- 6
nb <- 3000
nc <- 3
# Call JAGS from R (BRT 127 min)
js <- jags(jags.data,
inits,
parameters,
"D:/BoxFiles/Box Sync/CodigoR/Chinchilla/R/js-tempran.jags",
n.chains = nc,
n.thin = nt,
n.iter = ni,
n.burnin = nb)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 121346
## Unobserved stochastic nodes: 121416
## Total graph size: 10304333
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
##
## Burn-in phase, 3000 iterations x 3 chains
##
##
## Sampling from joint posterior, 5000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.summary(js, digits = 3)## Summary for model 'D:/BoxFiles/Box Sync/CodigoR/Chinchilla/R/js-tempran.jags'
## Saved parameters: psi mean.p sigma2 mean.phi N Nsuper b B deviance
## MCMC ran for 1238.596 minutes at time 2022-05-06 08:14:37.
##
## For each of 3 chains:
## Adaptation: 100 iterations (sufficient)
## Burn-in: 3000 iterations
## Thin rate: 6 iterations
## Total chain length: 8100 iterations
## Posterior sample size: 833 draws
##
## **WARNING** Rhat values indicate convergence failure.
##
## DIC info: (pD = var(deviance)/2)
## pD = 204.2 and DIC = 30268.71saveRDS(js, "D:/BoxFiles/Box Sync/CodigoR/Chinchilla/model/js_tempran.RData")
#a <- readRDS("model/cjs_add1.RData")Check model parameter and convergence
Notice: All Plots are in logit scale.
# nice plots
MCMCtrace(js,
params = c("psi", "mean.p", "sigma2"),
ISB = FALSE,
exact = TRUE,
pdf = FALSE,
main_den= c("Inclusion probability",
"Mean capture",
"Variation of phi"))
MCMCtrace(js,
params = c("mean.phi", "Nsuper"),
ISB = FALSE,
exact = TRUE,
pdf = FALSE,
main_den=c("mean survival",
"Size of superpopulation"))
MCMCplot(js,
params = c("psi", "mean.p", "sigma2", "mean.phi" ),
ci = c(50, 90),
guide_lines = TRUE,
main= "Parameters")
nl <- js
# Code to produce Fig. 10-6
# Calculate per-capita recruitment
T <- dim(leis)[2]
f <- matrix(NA, ncol = T, nrow = length(nl$sims.list$B[,1]))
for (t in 1:(T-1)){
f[,t] <- nl$sims.list$B[,t+1] / nl$sims.list$N[,t+1]
}
n.lower <- n.upper <- f.lower <- f.upper <- f.mean <- numeric()
for (t in 1:T){
n.lower[t] <- quantile(nl$sims.list$N[,t], 0.025)
n.upper[t] <- quantile(nl$sims.list$N[,t], 0.975)
}
for (t in 1:(T-1)){
f.lower[t] <- quantile(f[,t], 0.025)
f.upper[t] <- quantile(f[,t], 0.975)
f.mean[t] <- mean(f[,t])
}
par(mfrow = c(1, 2))
plot(nl$mean$N, type = "b", pch = 19, ylab = "Population size", xlab = "", axes = F, cex = 1.5, ylim = c(10, max(n.upper)))
axis(1, at = seq(1, T, 2), labels = seq(1987, 2020, 2))
axis(1, at = 1:T, labels = rep("", T), tcl = -0.25)
axis(2, las = 1)
segments(1:T, n.lower, 1:T, n.upper)
plot(f.mean, type = "b", pch = 19, ylab = "per capita recruitment", xlab = "", axes = F, cex = 1.5, ylim = c(0, 0.4))
axis(1, at = seq(1, (T-1), 2), labels = seq(1987, 2020, 2))
axis(1, at = 1:(T-1), labels = rep("", T-1), tcl = -0.25)
axis(2, las = 1)
segments(1:(T-1), f.lower, 1:(T-1), f.upper)
Check if estimates of mean survival and its temporal variability are close to those obtained under the CJS model.And also if the capture probabilities differ from CJS.
print(sessionInfo(), locale = FALSE)## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
##
## Matrix products: default
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] MCMCvis_0.15.3 jagsUI_1.5.2 AHMbook_0.2.3 readxl_1.3.1 ggplot2_3.3.5
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.1 terra_1.5-21 xfun_0.26
## [4] bslib_0.2.5.1 purrr_0.3.4 lattice_0.20-44
## [7] colorspace_2.0-2 vctrs_0.3.8 generics_0.1.0
## [10] htmltools_0.5.1.1 yaml_2.3.5 utf8_1.2.2
## [13] rlang_1.0.2 jquerylib_0.1.4 pillar_1.7.0
## [16] glue_1.6.2 withr_2.2.0 DBI_1.1.2
## [19] sp_1.4-5 lifecycle_1.0.1 plyr_1.8.6
## [22] stringr_1.4.0 munsell_0.5.0 gtable_0.3.0
## [25] cellranger_1.1.0 raster_3.5-15 mvtnorm_1.1-2
## [28] coda_0.19-4 codetools_0.2-18 evaluate_0.14
## [31] knitr_1.36 parallel_4.0.3 fansi_1.0.2
## [34] highr_0.9 Rcpp_1.0.8 scales_1.1.1
## [37] plotrix_3.8-1 unmarked_1.1.1 jsonlite_1.8.0
## [40] rjags_4-10 digest_0.6.27 stringi_1.7.6
## [43] dplyr_1.0.7 grid_4.0.3 cli_3.0.1
## [46] tools_4.0.3 magrittr_2.0.2 sass_0.4.0
## [49] tibble_3.1.6 crayon_1.5.0 pkgconfig_2.0.3
## [52] ellipsis_0.3.2 MASS_7.3-54 RandomFieldsUtils_0.5.3
## [55] RandomFields_3.3.8 assertthat_0.2.1 rmarkdown_2.11
## [58] rstudioapi_0.13 R6_2.5.1 compiler_4.0.3