Estimating Age at Recruitment

95 views
Skip to first unread message

Michael Roast

unread,
Aug 12, 2021, 8:58:22 AM8/12/21
to nimble-users
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)
})

Michael Roast

unread,
Aug 12, 2021, 8:59:43 AM8/12/21
to nimble-users
These are my current model estimates for time-variant recruitment.
Estimates.png

Michael Roast

unread,
Aug 12, 2021, 9:28:50 AM8/12/21
to nimble-users
To clarify,  I'm trying to estimate the probability an individual will breed for the first time at each given age. So for example I would like the following %s, 0% probability of a non-breeder (chicks and immature lumped) becoming a breeder </= Age 4, 30% at Age 5, 60% at Age 6, 10% at Age 7, 0% at Age >/= 8. By 'recruitment' I mean specifically transitioning into breeding adult status. The distinction between chick and immature isn't really important here other than that birds ringed as chicks are of known-age.

Thanks :)

Michael

Perry de Valpine

unread,
Aug 16, 2021, 12:35:43 PM8/16/21
to Michael Roast, nimble-users
Hi Michael,

Thanks for posting this question and example.

If I follow, your questions involve aspects of both statistical modeling issues (formulating the model and what can reasonably be estimated from the data) and implementation issues (how to set it up in nimble).  For the former, do you know the hmecology (hierarchical modeling in ecology) list?  People often point out the latest in relevant methods on topics such as this.  Our list here is more on the latter (nimble software issues per se).

I take the liberty of copying here from Marc Kéry's email footer from hmecology posts.  He points to three helpful lists for ecological models:

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2020)
---

Have you looked at the nimbleEcology package?  The dHMM and dDHMM give (dynamic) hidden Markov model marginal distributions.  This would allow you to set up a multi-event capture-recapture model in nimble related to how you would do it E-SURGE, I think.  There is a package vignette that I hope is helpful.

The issue of unknown age-within-stage is an issue I've been involved with in other forms in several projects.  I'm not sure if this is where you are going, but it sounds like you would want to include individuals of both known age (first captured as chicks) and unknown age (first captured as immature non-breeders).  You're saying the can set up a model for the former on its own, but to include the latter I suspect the only way to do full bookkeeping would be a hidden Markov model with age-stage combinations as the states.  With care I think one could set up transition matrices from each age-stage category to each other one and observation matrices from each age-stage category to each stage (observable) one.

Let us know how it goes.

-Perry


--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/9e58134a-5898-4e44-866c-b71a8ee57470n%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages