Random effects N-mixture models with Nimble

[david ferrer]

Hi everyone,

I’m trying to apply an N‐mixture model using Nimble with camera trap data. We started with the standard N‐mixture model that assumes an ecological process following a Poisson distribution and an observational process modeled as binomial. However, we observed overdispersion in this model. Because of that, we plan to incorporate a random effect into the ecological process (one for each camera location).

The first model produced results identical to those from a frequentist approach, so we are confident in its performance. The problem arises when we add the random effect in the second model. In that case, we encounter two issues: first, some parameters fail to initialize (which appears to be related to NAs in the observational count data), and second, the posterior distributions show poor convergence (as evidenced by the traceplots, Rhat, and effective sample size values). Even when we set initial values for the missing counts, the convergence issues persist.

Interestingly, when we use the stan_pcount function—with a Poisson ecological process and a random effect for each camera—the model runs smoothly (i.e. no initialization issues and good convergence for all parameters). Since stan_pcount has some limitations that prevent us from implementing all desired changes, we want to continue developing the model in Nimble.

Below, I share the code for both models as applied:
• Model 1 (N‐mixture model with Poisson–Binomial)
• Model 2 (N‐mixture model with Poisson (and random effect per camera) and Binomial)

###############

#### Create the model 1 -----------------------------------------------------------
library(nimble)

#Save a list with all the constants that we will employ and other list with the data.
constants <- list(nsurveys = ncol(counts.dry), # number of surveys
                  nsites = nrow(counts.dry), # number of sites
                  ngr = 56) #number of categories of cameras (cameras ID)


data <- list(counts = counts.dry, # Our animals counts (detections)
             d_rivers = cov.site$DIST.RIVER, #Covriate ecological process
             effort = cov.obs_dry #Covariate observational process
)

# Define the model
Poisson01NMIX <- nimbleCode({
 
  # Likelihood
 
  # ECOLOGICAL PROCESS
 
  # For each sampling site...
  for(i in 1:nsites) {
    N[i] ~ dpois(lambda[i])
   
    log(lambda[i]) <- b0 + b_d.river * d_rivers[i]
  }
 
  # OBSERVATIONAL PROCESS
 
  # For each sampling site...
  for(i in 1:nsites) {
    # and sampling survey...
    for(j in 1:nsurveys) {
      counts[i, j] ~ dbin(p[i,j], N[i])
     
      logit(p[i,j]) <- a0 + a_effort * effort[i,j]
    }
  }
 
  # PRIORS:
 
  b_d.river ~ dnorm(mean = 0, sd = 1)
  b0 ~ dnorm(mean = 0, sd = 1)
 
  a0 ~ dnorm(mean = 0, sd = 1)
  a_effort ~ dnorm(mean = 0, sd = 1)
 
}) 

# list of initial values (one for each MCMC chain)
Nst <- apply(counts.dry, 1, max, na.rm = TRUE)+1 #Important to give good initial values for latent N (not start with a lower values than present yet in any count)
inits <-function() {
  list(N = Nst,
       b0 = runif(1, -1, 1),
       b_d.river = runif(1,-1, 1),
       a0 = runif(1, -1, 1),
       a_effort = runif(1, -1, 1))
}

#### Create the model 2 -----------------------------------------------------------
library(nimble)

#Save a list with all the constants that we will employ and other list with the data.
constants <- list(nsurveys = ncol(counts.dry), # number of surveys
                  nsites = nrow(counts.dry), # number of sites
                  ngr = 56, #number of categories of cameras (cameras ID)
                  camID = cov.site$Camera.ID) # number camera ID, employed as random factor (one category per camera)


data <- list(counts = counts.dry, # Our animals counts (detections)
             d_rivers = cov.site$DIST.RIVER, #Covariate ecological process
             effort = cov.obs_dry #Covariate observational process
)

# Define the model
Poisson02NMIX <- nimbleCode({
 
  # Likelihood
 
  # ECOLOGICAL PROCESS
 
  # For each sampling site...
  for(i in 1:nsites) {
    N[i] ~ dpois(lambda[i])
   
    log(lambda[i]) <- b0[camID[i]] + b_d.river * d_rivers[i]
  }
 
  # OBSERVATIONAL PROCESS
 
  # For each sampling site...
  for(i in 1:nsites) {
    # and sampling survey...
    for(j in 1:nsurveys) {
      counts[i, j] ~ dbin(p[i,j], N[i])
     
      logit(p[i,j]) <- a0 + a_effort * effort[i,j]
    }
  }
 
  # PRIORS:
 
  b_d.river ~ dnorm(mean = 0, sd = 1)
  

  for(k in 1:ngr) {
    b0[k] ~ dnorm(mean = mu.b0, sd = sd.b0)
  }

  mu.b0 ~ dnorm(mean = 0, sd = 1)
  sd.b0 ~ dunif(0,10)
 
  a0 ~ dnorm(mean = 0, sd = 1)
  a_effort ~ dnorm(mean = 0, sd = 1)
 
}) 

# list of initial values (one for each MCMC chain)
Nst <- apply(counts.dry, 1, max, na.rm = TRUE)+1 #Important to give good initial values for latent N (not start with a lower values than present yet in any count)
inits <-function() {
  list(N = Nst,
       b0 = runif(56, -1, 1),
       mu.b0 = runif(1, -1, 1),

       sd.b0 = runif(1, 0, 10),

       b_d.river = runif(1,-1, 1),
       a0 = runif(1, -1, 1),
       a_effort = runif(1, -1, 1))
}

###############

For the secon model, when we run 

"####

modelCAM <- nimbleModel(code = PoissonNMIX,

+                         data = data,
+                         constants = constants,
+                         inits = inits(),

+                         calculate = FALSE)

####" 

we obtain the following note: "[Note] This model is not fully initialized. This is not an error. To see which variables are not initialized, use model$initializeInfo(). For more information on model initialization, see help(modelInitialization)."

Checking with "modelCAM$initializeInfo()"  yields:

  "[Note] Missing values (NAs) or non-finite values were found in model variables: lambda, p."

Based on that note, we attempted to set initial values for the NA counts by replacing NA with the maximum value observed at each site: 

"####

# list of initial values (one for each MCMC chain)

Nst <- apply(counts.dry, 1, max, na.rm = TRUE)+1 #Important to give good initial values for latent N (not start with a lower values than present yet in any count)

# Change NA in counts.dry with the maximum value detected on the site
counts_ini <- t(apply(counts.dry, 1, function(r) {
  # r is the camera corresponding vector with al the counts
  r[is.na(r)] <- max(r, na.rm = TRUE)
  return(r)
}))

inits <- function() {
  list(N = Nst,
       counts = counts_ini,
       mu.b0 = runif(1, -1, 1),
       b0 = runif(56, -1, 1),
       sd.b0 = runif(1, 0, 100),
       b_d.river = runif(1, -1, 1),
       a0 = runif(1, -1, 1),
       a_effort = runif(1, -1, 1))
}
####

"

However, we still encounter problems with initialization and convergence once the random effect is added, and we are unsure how to resolve these issues.

Thank you, 

David Ferrer

[Qing Zhao]

Hi David,

Would you mind sharing your data? I may have a look to see if I can help.

Best,

Qing

--
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 visit https://groups.google.com/d/msgid/nimble-users/e78968ff-8afd-47c3-bce0-2cf83126e10an%40googlegroups.com.

[david ferrer]

Hi Quing, 

Of course, here you have the environment required and the script. Here we have three models. Model 1 with the standard N-mixture. Model 2 with the random factor, and incorporating initial values to counts with NAs. And model 3 with random factor and simulating NAs of counts. 

Thank you, 

David

[Qing Zhao]

Hi David & group,

I tried a few things. There are still two parameters that don't quite converge (panel 1 = a0, the intercept of detection, panel 3 = b0, the intercept of abundance), and I'm not sure how to make them converge at this point, so I'd like to see what others do.

I suspected that these two parameters confound with each other, but when I plotted them, there is this negative correlation between them but it's not strong at all.

Code is attached in case anyone is interested in what I did.

Best,

Qing

To view this discussion visit https://groups.google.com/d/msgid/nimble-users/59437065-b7b7-4bca-a83b-cd1480163568n%40googlegroups.com.

[david ferrer]

Hi Qing and group, 

Thank you for your message. The improvements are clear. Regarding the correlation between the intercepts, I will check whether considering other parameters in the abundance and observation processes could be beneficial.

I appreciate the changes made, and I will now take the time to assimilate them and fully understand this procedure. Along with some private comments I have received, I will test the models. As soon as I have results, I will respond in this conversation.

Best regards,
David.

[Perry de Valpine]

Hi David and Qing,

Qing thanks for providing further insight.

I haven't looked at the code, but it may be relevant that sometimes detection and abundance are not really separable in N-mixture estimation. Aspects of this have been described by Kéry and by Knape et al. Ben Goldstein and I included one approach to the issue here, although that was not involving MCMC.

On the software side: nimble generally needs you to provide reasonable initial values. If you don't provide inits, nimble will attempt to simulate from the prior. With something like an N-mixture model, this can give invalid results, for example if the N at a site is less than one of the counts at that site, which can happen because the prior just says N is Poisson (or negative binomial). That is impossible and so you'll get bad results. This won't be an issue in Stan because you must be marginalizing. (Also if you provide missing data in a node nimble, it will be sampled, so it is best not to do this unless you actually need those imputed values as estimates of something unobserved.)

Your report that Stan appeared to mix where nimble did not is of interest. Can you clarify if you believe the two models were fully equivalent? Since Stan can't sample over discrete latent states, I assume you are marginalizing. In doing so, are you manually setting the range of the discrete summation? If so, that could spuriously appear to stabilize the results. nimbleEcology has marginal N-mixture distributions with options either to have the summation range determined dynamically or to set it manually. For use with nimbleHMC, you must set it manually. It may be of interest to compare to that.

HTH!

Perry

To view this discussion visit https://groups.google.com/d/msgid/nimble-users/CAE-HdUEZjnjXy5hYRMQmGHpu4wsJ721hccZPUei9-vHTupNTCg%40mail.gmail.com.

[Jose.J...@uclm.es]

Dear Perry, David...

With appropriate inits (including missing data inits) and less restrictive priors for the fixed effects, David's model converges reasonably well. It's not perfect, but neither does it seem to be the case you pointed out about the difficulty of separating detection and abundance.

Best regards

José