Capture-Recapture CJS of Darwin’s leaf-eared mouse (Phyllotis darwini)
Fixed_group_and_time_effects
Paola Correa, Diego Lizcano, Salvador Mandujano, Grabriel Andrade-Ponce et al.
22/4/2022
7. Estimation of survival probabilities using capture-recapture data
7.6. Models with time and group effects
7.6.1. Fixed group and time effects
Model from the book “Bayesian Population Analysis using WinBUGS — A Hierarchical Perspective” (2012) by Marc Kéry and Michael Schaub.
Adapted to stan by https://github.com/stan-dev/example-models/tree/master/BPA
State process: the latent alive/dead state
Observation process: relates true state to observed state
phi ~ dunif(0, 1) # Apparent survival.
p ~ dunif(0, 1) # Recapture.
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).
Cormack-Jolly-Seber models give the probability of a capture history for each of many individuals, conditional on first capture, based on parameters for survival and detection probability
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.
Load data and run jags model
library(ggplot2)
library(readxl)
library(AHMbook)
library(jagsUI)
library(MCMCvis)
## Read data
## The data generation code is in bpa-code.txt, available at
## http://www.vogelwarte.ch/de/projekte/publikationen/bpa/complete-code-and-data-files-of-the-book.html
# stan_data <- read_rdump("D:/BoxFiles/Box Sync/CodigoR/Chinchilla/R/cjs_add.data.R")
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 males because group 1 is males
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
# Pdarwini_data <- list(y=y1,
# n_occasions=n_occasions1,
# g=g1,
# group=group1,
# nind=nind1,
# df=df1,
# R=stan_data$R)
# Create vector with occasion of marking
get.first <- function(x) min(which(x!=0))
f <- apply(y1, 1, get.first)
# Specify model in BUGS language
cat("
model {
# Priors and constraints
for (i in 1:nind){
for (t in f[i]:(n.occasions-1)){
logit(phi[i,t]) <- beta[group[i]] + gamma[t]
p[i,t] <- p.g[group[i]]
} #t
} #i
# for survival parameters
for (t in 1:(n.occasions-1)){
gamma[t] ~ dnorm(0, 0.01)I(-10, 10) # Priors for time effects
phi.g1[t] <- 1 / (1+exp(-gamma[t])) # Back-transformed survival of males
phi.g2[t] <- 1 / (1+exp(-gamma[t]-beta[2])) # Back-transformed survival of females
}
beta[1] <- 0 # Corner constraint
beta[2] ~ dnorm(0, 0.01)I(-10, 10) # Prior for difference in male and female survival
# for recapture parameters
for (u in 1:g){
p.g[u] ~ dunif(0, 1) # Priors for group-spec. recapture
}
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,f[i]] <- 1
for (t in (f[i]+1):n.occasions){
# State process
z[i,t] ~ dbern(mu1[i,t])
mu1[i,t] <- phi[i,t-1] * z[i,t-1]
# Observation process
y[i,t] ~ dbern(mu2[i,t])
mu2[i,t] <- p[i,t-1] * z[i,t]
} #t
} #i
}
",fill = TRUE)##
## model {
##
## # Priors and constraints
## for (i in 1:nind){
## for (t in f[i]:(n.occasions-1)){
## logit(phi[i,t]) <- beta[group[i]] + gamma[t]
## p[i,t] <- p.g[group[i]]
## } #t
## } #i
## # for survival parameters
## for (t in 1:(n.occasions-1)){
## gamma[t] ~ dnorm(0, 0.01)I(-10, 10) # Priors for time effects
## phi.g1[t] <- 1 / (1+exp(-gamma[t])) # Back-transformed survival of males
## phi.g2[t] <- 1 / (1+exp(-gamma[t]-beta[2])) # Back-transformed survival of females
## }
## beta[1] <- 0 # Corner constraint
## beta[2] ~ dnorm(0, 0.01)I(-10, 10) # Prior for difference in male and female survival
## # for recapture parameters
## for (u in 1:g){
## p.g[u] ~ dunif(0, 1) # Priors for group-spec. recapture
## }
##
## # Likelihood
## for (i in 1:nind){
## # Define latent state at first capture
## z[i,f[i]] <- 1
## for (t in (f[i]+1):n.occasions){
## # State process
## z[i,t] ~ dbern(mu1[i,t])
## mu1[i,t] <- phi[i,t-1] * z[i,t-1]
## # Observation process
## y[i,t] ~ dbern(mu2[i,t])
## mu2[i,t] <- p[i,t-1] * z[i,t]
## } #t
## } #i
## }# Initial values
# In JAGS we have to give good initial values for the latent state z. At all occasions when an individual was observed, its state is z = 1 for sure. In addition, if an individual was not observed at an occasion, but was alive for sure, because it was observed before and thereafter (i.e. has a capture history of e.g. {101} or {10001}), then we know that the individual was alive at all of these occasions, and thus z = 1. Therefore, we should provide initial values of z = 1 at these positions as well. The following function provides such initial values from the observed capture histories:
known.state.cjs <- function(ch){
state <- ch
for (i in 1:dim(ch)[1]){
n1 <- min(which(ch[i,]==1))
n2 <- max(which(ch[i,]==1))
state[i,n1:n2] <- 1
state[i,n1] <- NA
}
state[state==0] <- NA
return(state)
}
# (Note that the function known.state.cjs is used in section 7.3.1 as well for another purpose)
# Bundle data
jags.data <- list(y = y1, f = f, nind = dim(y1)[1], n.occasions = dim(y1)[2], z = known.state.cjs(y1), g = length(unique(group1)), group = group1)
# Function to create a matrix of initial values for latent state z
cjs.init.z <- function(ch,f){
for (i in 1:dim(ch)[1]){
if (sum(ch[i,])==1) next
n2 <- max(which(ch[i,]==1))
ch[i,f[i]:n2] <- NA
}
for (i in 1:dim(ch)[1]){
ch[i,1:f[i]] <- NA
}
return(ch)
}
# Initial values
inits <- function(){list(z = cjs.init.z(y1, f), gamma = rnorm(n_occasions1-1), beta = c(NA, rnorm(1)), p.g = runif(length(unique(group1)), 0, 1))}
# Parameters monitored
parameters <- c("phi.g1", "phi.g2", "p.g", "beta")
# MCMC settings
ni <- 10000
nt <- 3
nb <- 2000
nc <- 3
# Call JAGS from R (BRT 7 min)
cjs.add <- jags(data = jags.data,
inits = inits,
parameters.to.save = parameters,
model.file= "D:/BoxFiles/Box Sync/CodigoR/Chinchilla/R/cjs-add.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: 57460
## Unobserved stochastic nodes: 56824
## Total graph size: 282094
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
##
## Burn-in phase, 2000 iterations x 3 chains
##
##
## Sampling from joint posterior, 8000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.# Summarize posteriors
summary(cjs.add, digits = 3) ## Summary for model 'D:/BoxFiles/Box Sync/CodigoR/Chinchilla/R/cjs-add.jags'
## Saved parameters: phi.g1 phi.g2 p.g beta deviance
## MCMC ran for 179.778 minutes at time 2022-04-30 20:28:33.
##
## For each of 3 chains:
## Adaptation: 100 iterations (sufficient)
## Burn-in: 2000 iterations
## Thin rate: 3 iterations
## Total chain length: 10100 iterations
## Posterior sample size: 2666 draws
##
## Successful convergence based on Rhat values (all < 1.1).
##
## DIC info: (pD = var(deviance)/2)
## pD = 1675 and DIC = 3426.844saveRDS(cjs.add, "D:/BoxFiles/Box Sync/CodigoR/Chinchilla/model/cjs_add1.RData")Check model parameter and convergence
Notice: All Plots are in logit scale.
# nice plots
MCMCtrace(cjs.add,
params = c("p.g[1]", "p.g[2]"),
ISB = FALSE,
exact = TRUE,
pdf = FALSE,
main_den= c("recapture male",
"recapture female"))
MCMCtrace(cjs.add,
params = c("beta[2]"),
ISB = FALSE,
exact = TRUE,
pdf = FALSE,
main_den="Difference male-female survival")
MCMCplot(cjs.add,
params = c("p.g", "beta"),
ci = c(50, 90),
guide_lines = TRUE,
main= "Recap. and dif survival, male[1]")
# Figure of male and female survival
lower.f <- upper.f <- lower.m <- upper.m <- numeric()
for (t in 1:(n_occasions1-1)){
lower.f[t] <- quantile(cjs.add$sims.list$phi.g1[,t], 0.025)
upper.f[t] <- quantile(cjs.add$sims.list$phi.g1[,t], 0.975)
lower.m[t] <- quantile(cjs.add$sims.list$phi.g2[,t], 0.025)
upper.m[t] <- quantile(cjs.add$sims.list$phi.g2[,t], 0.975)
}
plot(x=(1:(n_occasions1-1))-0.1, y = cjs.add$mean$phi.g1, type = "b", pch = 16, ylim = c(-0.01, 1), ylab = "Survival probability", xlab = "Year", bty = "n", cex = 1.5, axes = FALSE, main="Temporal Survival. Black(male)")
axis(1, at = 1:34, labels = rep(NA,34), tcl = -0.25)
axis(1, at = seq(2,34,2), labels = c("1988","1990","1992","1994","1996","1998","2000","2002","2004","2006","2008","2010", "2012", "2014", "2016", "2018", "2020"))
axis(2, at = seq(0, 1, 0.1), labels = c("0",NA, "0.2", NA, "0.4", NA, "0.6", NA, "0.8", NA, "1.0"), las = 1)
segments((1:(n_occasions1-1))-0.1, lower.f, (1:(n_occasions1-1))-0.1, upper.f)
points(x = (1:(n_occasions1-1))+0.1, y = cjs.add$mean$phi.g2, type = "b", pch = 1, lty = 2, cex = 1.5)
segments((1:(n_occasions1-1))+0.1, lower.m, (1:(n_occasions1-1))+0.1, upper.m)
Survival of the two sexes varies in parallel over time, but on the logit scale. Hence, on the probability scale the two curves are not parallel—as the difference becomes smaller the closer the estimates are to 1 or 0.
7.6.2 Fixed Group and Random Time Effects
combine fixed group and random time effects to estimate temporal variability of survival (or recapture) in each group separately. As for the interacting model before, such a model would assume that the temporal variability of each group is independent of that in the other group(s).
# Bundle data
jags.data <- list(y = y1, f = f, nind = dim(y1)[1], n.occasions = dim(y1)[2], z = known.state.cjs(y1), g = length(unique(group1)), group = group1, R = matrix(c(5, 0, 0, 1), ncol = 2), df = 3)
# Initial values
inits <- function(){list(z = cjs.init.z(y1, f), p.g = runif(length(unique(group1)), 0, 1), Omega = matrix(c(1, 0, 0, 1), ncol = 2))}
# Parameters monitored
parameters <- c("eta.phi", "p.g", "Sigma", "mean.phi")
# MCMC settings
ni <- 10000
nt <- 3
nb <- 2000
nc <- 3
# Call JAGS from R (BRT 5 min)
cjs.corr <- jags(jags.data,
inits,
parameters,
"D:/BoxFiles/Box Sync/CodigoR/Chinchilla/R/cjs-temp-corr.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: 57460
## Unobserved stochastic nodes: 56826
## Total graph size: 281452
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase.....
## Adaptive phase.....
## Adaptive phase complete
##
##
## Burn-in phase, 2000 iterations x 3 chains
##
##
## Sampling from joint posterior, 8000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.# Summarize posteriors
summary(cjs.corr, digits = 3) ## Summary for model 'D:/BoxFiles/Box Sync/CodigoR/Chinchilla/R/cjs-temp-corr.jags'
## Saved parameters: eta.phi p.g Sigma mean.phi deviance
## MCMC ran for 98.87 minutes at time 2022-04-30 23:28:21.
##
## For each of 3 chains:
## Adaptation: 300 iterations (sufficient)
## Burn-in: 2000 iterations
## Thin rate: 3 iterations
## Total chain length: 10300 iterations
## Posterior sample size: 2666 draws
##
## Successful convergence based on Rhat values (all < 1.1).
##
## DIC info: (pD = var(deviance)/2)
## pD = 1550.6 and DIC = 3303.144saveRDS(cjs.corr, "D:/BoxFiles/Box Sync/CodigoR/Chinchilla/model/cjs.corr1.RData")Check model parameter and convergence
Notice: All Plots are in logit scale.
MCMCtrace(cjs.corr,
params = c("mean.phi[1]", "mean.phi[2]"),
ISB = FALSE,
exact = TRUE,
pdf = FALSE,
main_den= c("mean group survival male",
"mean group survival female"))
MCMCtrace(cjs.corr,
params = c("p.g[1]", "p.g[2]"),
ISB = FALSE,
exact = TRUE,
pdf = FALSE,
main_den= c("mean group detection male",
"mean group detection female"))
MCMCplot(cjs.corr,
params = c("p.g", "mean.phi"),
ci = c(50, 90),
guide_lines = TRUE,
main= "Recap. and mean survival, male[1]")
MCMCplot(cjs.corr,
params = c("Sigma"),
ci = c(50, 90),
guide_lines = TRUE,
main= "temporal covariances of survival")
The elements Σ12 and Σ21 are the temporal covariances of logit male and logit female survival. This quantity may not be easy to interpret, and we may want to compute the temporal correlation of male and female survival.
corr.coef <- cjs.corr$sims.list$Sigma[,1,2] / sqrt(cjs.corr$sims.list$Sigma[,1,1] * cjs.corr$sims.list$Sigma[,2,2])
df <- data.frame(corr.coef)
media <- data.frame(mean=mean(corr.coef))
p <-ggplot(df, aes(x = corr.coef)) +
geom_histogram() + geom_vline(aes(xintercept = mean(corr.coef)),
color = "#000000", size = 1.25,
linetype = "dashed")
p
mediaprint(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] labeling_0.4.2 knitr_1.36 parallel_4.0.3
## [34] fansi_1.0.2 highr_0.9 Rcpp_1.0.8
## [37] scales_1.1.1 plotrix_3.8-1 unmarked_1.1.1
## [40] jsonlite_1.8.0 farver_2.1.0 rjags_4-10
## [43] digest_0.6.27 stringi_1.7.6 dplyr_1.0.7
## [46] grid_4.0.3 cli_3.0.1 tools_4.0.3
## [49] magrittr_2.0.2 sass_0.4.0 tibble_3.1.6
## [52] crayon_1.5.0 pkgconfig_2.0.3 ellipsis_0.3.2
## [55] MASS_7.3-54 RandomFieldsUtils_0.5.3 RandomFields_3.3.8
## [58] assertthat_0.2.1 rmarkdown_2.11 rstudioapi_0.13
## [61] R6_2.5.1 compiler_4.0.3