HMM Capture Mark Recapture with Random effects for cohorts: numerical problem with estimation

[Etienne Rouby]

Hi all,

I am currently coding a Hidden Markov Model on Nimble. My model has two stages :

- PB: Prebreeder
- B: Breeder

One constraint is that the PB stage is 100% sure, from 1 to 5 years included. Then, from age 6, my birds can transit into the breeder stage, with a certain probability psiPB * phiPB. psiPB is the probability of reproducing as a prebreeder. PhiPB is the probability of surviving from one year to the next as a pre-Breeder. My simple model and the one with sexes as nested indexes work well.

Now I am including random effects to estimate differences in psiPB from one cohort to another. I get some weird estimations. And I do not know where it comes from.

First, Here is how I coded my random effect to estimate how the year of birth affects the psiPB parameter on a logit scale:

multistate.cate <- nimbleCode({
 
  # Observation matrix
  omega[1,1] <- 1-pPB
  omega[1,2] <- pPB
  omega[1,3] <- 0
 
  omega[2,1] <- 1-pB
  omega[2,2] <- 0
  omega[2,3] <- pB
 
  omega[3,1] <- 1
  omega[3,2] <- 0
  omega[3,3] <- 0
 
  for(i in 1:n_ind){
    for(t in (first[i]):(n_year-1)) {
     
      phiB[i,t] <- Mu.phiB
      phiPB[i,t] <- Mu.phiPB[age[i,t]]
      logit(psiPB[i,t]) <- logit(Mu.psiPB[age[i,t]]) + epsilon[BIRTH_YEAR[i]]
     
      # Transition matrix
      gamma[1,1,i,t] <- phiPB[i,t] * (1-psiPB[i,t])
      gamma[1,2,i,t] <- phiPB[i,t] * psiPB[i,t]
      gamma[1,3,i,t] <- 1-phiPB[i,t]
     
      gamma[2,1,i,t] <- 0
      gamma[2,2,i,t] <- phiB[i,t]
      gamma[2,3,i,t] <- 1-phiB[i,t]
     
      gamma[3,1,i,t] <- 0
      gamma[3,2,i,t] <- 0
      gamma[3,3,i,t] <- 1
    }
  }
 
  # Priors
  pPB ~ dunif(0,1)
  pB ~ dunif(0,1)
  Mu.phiB ~ dunif(0,1)
 
  for(a in 1:5){
    Mu.phiPB[a] ~ dunif(0,1)
    Mu.psiPB[a] ~ dunif(0,0.01)
  }
 
  for(a in 6:Amax){
    Mu.phiPB[a] ~ dunif(0,1)
    Mu.psiPB[a] ~ dunif(0,1)
  }
 
  for(birth_year in unique(BIRTH_YEAR)) {
    espilon[birth_year] ~ dnorm(0, sigma_eps)
  }
  sigma_eps ~ dunif(0, 1)
 
  # Likelihood
  for(i in 1:n_ind){
    # latent state at first capture
    z[i,first[i]] <- y[i,first[i]] - 1
    for(t in (first[i]+1):n_year){
      # z[i,t] given
      z[i,t] ~ dcat(gamma[z[i,t-1],1:3,i,t-1])
      y[i,t] ~ dcat(omega[z[i,t],1:3])
    }
  }
 
})

And the initial values that I pass to my model :

initial.values <- function(){list(phiPB = matrix(runif(ncol(y)*nrow(y), 0, 1),
                                                 ncol = ncol(y),
                                                 nrow = nrow(y)),
                                  phiB = matrix(runif(ncol(y)*nrow(y), 0, 1),
                                                ncol = ncol(y),
                                                nrow = nrow(y)),
                                  psiPB = matrix(runif(ncol(y)*nrow(y), 0, 1),
                                                 ncol = ncol(y),
                                                 nrow = nrow(y)),
                                  Mu.phiPB = c(runif(seq(1, max(age, na.rm = TRUE)), 0, 1)),
                                  Mu.psiPB = c(runif(seq(1, 5), 0, 0.001), runif(seq(6, max(age, na.rm = TRUE)), 0, 1)),
                                  Mu.phiB = 0.9,
                                  pPB = 0.8,
                                  pB = 0.9,
                                  epsilon = c(rnorm(seq(1, max(birth_year)), 0, runif(1,0,1))),
                                  sigma_eps = runif(1,0,1),
                                  z = zinits)}

Second, in the attached files, this is what I obtained for my parameters epsilon and sigma_eps. I also get the following message: NAs were detected in model variables: epsilon, logProb_espilon. This means that these parameters have problems to initialize i think.

I tried some things to figure it out what the problem is:
- Give more data and less years (so less epsilons to estimate, currently I just have 2)
- Change the prior scale and initial values for epsilon and sigma_eps.

Currently, nothing has worked. I think this is a problem of parametrization. But I need to figure out what to explore next to fix the problem and learn from my mistakes here.

Do you have a clue or something in mind that could help me to fix the problem(s) ?

Best and thank you very much for your reading. Do not hesitate if you need more input from my part.

Etienne

[Etienne Rouby]

UPDATE :

As noticed by a community member, there is a misconception in my code currently. I am defining a prior for epsilon which is :

for(birth_year in unique(BIRTH_YEAR)) {
    espilon[birth_year] ~ dnorm(0, sigma_eps)
  }
  sigma_eps ~ dunif(0, 1)

At the same time, I am initiating values for epsilon. Which makes no sense. So I changed this and now I do not initiate some values for epsilon. I still have a problem, which is the following:
NAs were detected in model variables: epsilon, espilon, logProb_espilon, psiPB, gamma, logProb_z.

I think it is a matter of how I define the epsilon index or the distribution around epsilon. Since I have a birth_year for all individuals in my dataset, there are no NAs in the index. I do not think it is a problem of sigma because when I try to debug, sigma does not seem to be associated with a particular bug. I also tried to run the model the long way by using nimbleModel(), but I got the same information.

Do you have any advice ? If I find something, I keep people updated here.

Best,
Etienne

[Perry de Valpine]

Dear Etienne,

I see a few things that might be going on that I can speculate about.

The first is that for(birth_year in unique(BIRTH_YEAR)) is not supported in model code. Only sequential indices can be used. So you would need something like:

for(i in 1:num_birth_years) {

espilon[unique_birth_year[i]] ~ dnorm(0, sigma_eps)

}

where you would create 

unique_birth_year <- unique(BIRTH_YEAR)

num_birth_years <- length(unique_birth_year)

 in R and provide them in the constants list.

The second is that you might be assuming some variables are visible (in scope) when the initial.values function is called that might not necessarily be in scope. For example, ncol(y) depends on having y in scope, max(birth_year) depends on having birth_year in scope, etc. Maybe they are. I just can't tell. Even worse, the function might be seeing values of those variables using R's lexical scoping that are in (e.g.) the R global environment and could even be different from the ones you have in the data and/or constants lists.

The third is that 

 logit(psiPB[i,t]) <- logit(Mu.psiPB[age[i,t]]) + epsilon[BIRTH_YEAR[i]]

can potentially result in NA if Mu.psiPB[age[i,t]] happens to be too close to 0 or epsilon[BIRTH_YEAR[i]] happens to be invalid.

My guess is that the invalid for() construction is causing multiple of these problems.

HTH

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/c1a71507-8a62-4ef6-a21f-e4e601d375ffn%40googlegroups.com.