# 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)