I have a capture-recapture dataset for a long-lived seabird species with delayed recruitment and I would like to estimate the probability of recruitment at different ages. I have read papers that apparently are able to do this in a multi-event model framework (e.g. Genovart et al. 2016 J. Appl. Ecol.), but these were implemented in E-SURGE software, and I am struggling to work out how to implement and code similar models using Nimble.
My data set contains ~1800 individuals that were first captured as chicks (~45%; known-age), immature non-breeders (~10%; unknown-age) and breeding adults (~45%; unknown-age). The dataset spans 19 timepoints (years) and has variable but generally increasing survey effort. In this species, chicks fledge, and are typically not easily observed for several years while immature non-breeder, before successfully recruiting to become breeding adults. Occasionally immature non-breeders are detected, with some uncertainty around breeding status. Chicks captured in most recent (4-5) years have not yet returned as adult breeders.
Up to now, I have a model that seems to work (below), that estimates time-variant transition probability from immature non-breeder (NB) to breeding adult (BA) (i.e. recruitment), by lumping chicks and non-breeders into a single NB category. Recruitment varies from 0.01-0.13 each year, while in a time-invariant model recruitment probability is estimated as 0.04. (This seems infeasibly low?)
I have tried to modify this model in different ways to estimate the probability of recruitment at different ages. Mainly, I have tried to use the subset of birds banded as chicks that are of known age, and tried to code the time-variant estimation of state transition so estimates are relative to the first capture of each known-age individual, such that the time-variant transition estimates equate to probability of transition at each age. No other models beyond this model appear to run successfully or produce anything meaningful.
I'm still new to not only Nimble, but also this type of capture-recapture data analysis, so there is a good chance I'm missing something fundamental too.
If anyone can offer any advice on how I can estimate probability of recruitment at different ages, or age-dependent state transition with Nimble it would be greatly appreciated! Thanks in advance.
Michael
# States:
# 1. NB Alive Immature Non-breeder (Chicks, and individuals yet to breed)
# 2. BA Alive Breeding Age Adult (including sabbatical non-breeders)
# 3. D Dead
# Observations:
# 1. ND Non-detected
# 2. NB Alive Non-breeder (Individuals not yet observed to breed)
# 3. BA Alive Breeding Age Adult (Successful, failed, or on sabbatical year)
#Year of first capture per individual
first <- apply(y, 1, function(x) min(which(x!=1)))
#Model
hmm.phi.p.psit <- nimbleCode({
# Priors
phiNB ~ dunif(0,1) # prior for Chick/Immature Non-breeder survival φNB
phiBA ~ dunif(0,1) # prior for Breeding Adult survival φBA
pNB ~ dunif(0,1) # prior for probability of detecting immature NB
pBA ~ dunif(0,1) # prior for probability of detecting breeding adult
betaNB ~ dunif(0,1) # prior for probability of correct assignment to NB
betaBA ~ dunif(0,1) # prior for probability of correct assignment to BA
pi ~ dunif(0,1) # prior for probability of NB at initial capture
# Vector of initial state probabilities δ
delta[1] <- pi # Pr(Initial capture = NB)
delta[2] <- 1-pi # Pr(Initial capture = BA)
delta[3] <- 0 # Pr(Initial capture = D)
for(t in 1:(T-1)){
psiNBBA[t] ~ dunif(0,1) # prior for probability of NB transitioning to BA (recruitment)
# Transition matrix Γ
gamma[1,1,t] <- phiNB * (1-psiNBBA[t]) # transition Pr(NB t --> NB t+1)
gamma[1,2,t] <- phiNB * psiNBBA[t] # transition Pr(NB t --> BA t+1)
gamma[1,3,t] <- 1 - phiNB # transition Pr(NB t --> D t+1)
gamma[2,1,t] <- 0 # transition Pr(BA t --> NB t+1)
gamma[2,2,t] <- phiBA # transition Pr(BA t --> BA t+1)
gamma[2,3,t] <- 1-phiBA # transition Pr(BA t --> D t+1)
gamma[3,1,t] <- 0 # transition Pr(D t --> NB t+1)
gamma[3,2,t] <- 0 # transition Pr(D t --> BA t+1)
gamma[3,3,t] <- 1 # transition Pr(D t --> D t+1)
}
# Observation matrix Ω
omega[1,1] <- 1-pNB # observation Pr(NB t --> not-detected t)
omega[1,2] <- pNB * betaNB # observation Pr(NB t --> detected NB t)
omega[1,3] <- pNB * (1-betaNB) # observation Pr(NB t --> detected BA t)
omega[2,1] <- 1-pBA # observation Pr(BA t --> not-detected t)
omega[2,2] <- pBA * (1-betaBA) # observation Pr(BA t --> detected NB t)
omega[2,3] <- pBA * betaBA # observation Pr(BA t --> detected BA t)
omega[3,1] <- 1 # observation Pr(D t --> not-detected t)
omega[3,2] <- 0 # observation Pr(D t --> detected NB t)
omega[3,3] <- 0 # observation Pr(D t --> detected BA t)
# Observation matrix Ω at first capture
omega.init[1,1] <- 0 # observation Pr(NB t --> not-detected t)
omega.init[1,2] <- betaNB # observation Pr(NB t --> detected NB t)
omega.init[1,3] <- 1-betaNB # observation Pr(NB t --> detected BA t)
omega.init[2,1] <- 0 # observation Pr(BA t --> not-detected t)
omega.init[2,2] <- 1-betaBA # observation Pr(BA t --> detected NB t)
omega.init[2,3] <- betaBA # observation Pr(BA t --> detected BA t)
omega.init[3,1] <- 1 # observation Pr(D t --> not-detected t)
omega.init[3,2] <- 0 # observation Pr(D t --> detected NB t)
omega.init[3,3] <- 0 # observation Pr(D t --> detected BA t)
# Likelihood
for (i in 1:N){
# latent state at first capture
z[i,first[i]] ~ dcat(delta[1:3]) # y[i,first[i]]-1 # Draws an initial state
y[i,first[i]] ~ dcat(omega.init[z[i,first[i]],1:3])
for (j in (first[i]+1):T){
z[i,j] ~ dcat(gamma[z[i,j-1],1:3,j-1]) # z(t) given z(t-1)
y[i,j] ~ dcat(omega[z[i,j],1:3]) # y(t) given z(t)
}
}
})
# Defines constants in the model
my.constants <- list(N = nrow(y), # N individual
T = ncol(y), # N of time points
first = first) # occasions of first capture
# Defines data in the model
# 1 is non-detected, higher numbers are detected states
my.data <- list(y=y)
# Defines initial values for states and probabilities
# Initial latent states are 1 (NB), until 2 (BA), then 2 thereafter
zinits <- y-1
for (i in 1:nrow(zinits)) {
for (j in 1:ncol(zinits)) {
if(j < first[i]){zinits[i,j] <- 1}
if(zinits[i,j] == 2){zinits[i,j:ncol(zinits)] <- 2}
if(zinits[i,j] == 0){zinits[i,j] <- 1}
}}
zinits <- as.matrix(zinits)
# Defines initial values for the model
initial.values <- function() list(phiNB = runif(1,0,1),
phiBA = runif(1,0,1),
pNB = runif(1,0,1),
pBA = runif(1,0,1),
betaNB = runif(1,0,1),
betaBA = runif(1,0,1),
psiNBBA = runif(1,0,1),
pi = runif(1,0,1),
z = zinits)
# Specfies the parameters of interest to record in output
parameters.to.save <- c("phiNB","phiBA", "psiNBBA","pNB","pBA","pi","betaNB","betaBA")
#MCMC model details
n.iter <- 20000
n.burnin <- 2000
n.chains <- 2
system.time({ # Wrapper to report model run-time.
mcmc.phi.p.psit <- nimbleMCMC(code = hmm.phi.p.psit, #model code defined matrices
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)
})