# Code for generating CJS capture histories
# Adapted from code from Marc Kéry and Michael Schaub's excellent book,
# "Bayesian population analysis using WinBUGS: A hierarchical perspective"
# https://www.vogelwarte.ch/de/projekte/publikationen/bpa/
# Define parameter values
n.occasions <- 6 # Number of capture occasions
n.mark <- 60
marked <- rep(n.mark, (n.occasions - 1)) # Number of newly marked individuals each year
# Vectors of values for phi and p for each year
# Each vector needs to have (n.occasions - 1) values and have values between 0 and 1
phi <- c(0.6, 0.4, 0.7, 0.5, 0.45)
p <- c(0.3, 0.29, 0.26, 0.33, 0.31)
# Define matrices with survival and recapture probabilities
PHI <- matrix(phi, ncol = n.occasions-1, nrow = sum(marked), byrow = TRUE)
P <- matrix(p, ncol = n.occasions-1, nrow = sum(marked), byrow = TRUE)
# Define function to simulate a capture-history (CH) matrix
simul.cjs <- function(PHI, P, marked){
n.occasions <- dim(PHI)[2] + 1
CH <- matrix(0, ncol = n.occasions, nrow = sum(marked))
# Define a vector with the occasion of marking
# note: length(marked) is the number of marking occasions
# so this next line indicates which occasion each individual was 1st marked on
mark.occ <- rep(1:length(marked), marked[1:length(marked)])
# Fill the CH matrix
for (i in 1:sum(marked)) {
CH[i, mark.occ[i]] <- 1 # Write a 1 at the release occasion
if (mark.occ[i] == n.occasions) next # ignore animals 1st caught on last occasion
for (t in (mark.occ[i] + 1):n.occasions) {
# Bernoulli trial: does individual survive occasion?
sur <- rbinom(1, 1, PHI[i, t - 1])
if (sur == 0)
break # If dead, move to next individual
# Bernoulli trial: is individual recaptured?
rp <- rbinom(1, 1, P[i, t - 1])
if (rp == 1) CH[i,t] <- 1
} #t
} #i
return(CH)
}
# Simulate capture-histories
CH <- simul.cjs(PHI, P, marked)
CH
# collapse the numbers in each row into character data
ch <- apply(CH, 1, paste, collapse = "")
# summarize how many animals have each observed capture history
table(ch)
# build and output capture history data in text file
# formatted for Program MARK
# adapted from addendum to Chapter 2 of Cooch and White
input.file <- data.frame(eh = ch)
# add the frequency column and semicolon at the end of each row
# as needed for a Program MARK input file, where the frequency
# value is the number of individuals represented the row
# Here, each row provides information for 1 animal
input.file$end <- "1 ;"
# write out the input file
write.table(input.file, file="cjs_data_hw02.inp",
sep = " ",
quote = FALSE,
col.names = FALSE,
row.names = FALSE)