Initializing Dynamic N-Mixture Models in NIMBLE

[heather...@gmail.com]

I have a multi-species, multi-site, dynamic n-mixture model that runs in JAGS but is incredibly slow. I'm trying to convert the same model to NIMBLE so that it will run faster. However, I'm really stuck on the initialization part. Using the Nimble initializeinfo() function I'm able to find the nodes that haven't been initialized, which I've been using to track down extra initial values to send to the model. But my model still will not produce results (that aren't all 0). 

I'm fairly new to NIMBLE so I'm hoping my problem is something basic that I've just overlooked.  I also don't fully understand which parameters need to have initial values and which ones do not. For instance, in my model below, I have this line:

log(gamma[i,k,q]) <- gamma.b2[k] + gamma.b3[k]*temp[i,q] +
              gamma.b4[k]*precip[i,q] + gamma.b5[k]*q + gamma.b6[k]*asp[i] +
                gamma.b7[k]*moisture[i]+ gamma.b8[group[k]]

I've provided reasonable initial values for all the gamma.b2 - gamma.b7 nodes, but NIMBLE says gamma is not fully initialized. I've tried running calculate() on gamma, but no luck. What am I missing?

 

My model code and files to run everything are below. I'm hoping someone can spot an obvious error that's not letting the model run properly. 

Thanks in advance! 

-Heather 

nimblebirds <-
  nimbleCode({
  for (i in 1:n.sites) {
     for (k in 1:n.species){
        log(psi[i,k,1]) <- psi.b2[k] +  psi.b3[k]*temp[i,1] + psi.b4[k]*precip[i,1] +
        psi.b5[k]*asp[i] + psi.b6[k]*moisture[i] + psi.b7[group[k]]
        abund[i,k,1] ~dpois(psi[i,k,1])
       
          for (q in 2:n.years){
             abund[i,k,q] <- S[i,k,q] + R[i,k,q]
              S[i,k,q] ~ dbin(phi[i,k,q], abund[i,k,q-1]) #reductions from pop
              R[i,k,q] ~ dpois(abund[i,k,q-1]*gamma[i,k,q]) #additions to pop
              log(gamma[i,k,q]) <- gamma.b2[k] + gamma.b3[k]*temp[i,q] +
              gamma.b4[k]*precip[i,q] + gamma.b5[k]*q + gamma.b6[k]*asp[i] +
                gamma.b7[k]*moisture[i]+ gamma.b8[group[k]]
             
              logit(phi[i,k,q]) <- phi.b2[k] + phi.b3[k]*temp[i,q] +
              phi.b4[k]*precip[i,q] + phi.b5[k]*q + phi.b6[k]*asp[i] +
                phi.b7[k]*moisture[i]+phi.b8[group[k]]
             
             
      } #end q
    } #end k
  } #end i
   
#detection loop

for (q in 1:n.years){
  
      for (k in 1:n.species){
      for (i in 1:n.sites) {
            p.dist.sum[i,k,q] <- sum(pi[i,k,q,1:bin.dist])
 
  obs.pis[i,k,q,1:bin.dist] <- p.dist.c[i,k,q,1:bin.dist]*surv[i,q]
  obs[i,k,q,1:bin.dist] ~ dmulti(obs.pis[i,k,q,1:bin.dist], n[i,k,q]) #what distance was it at when observed
  timeclass.pis[i,k,q,1:4] <- pi.pa.c[i,k,q,1:4]*surv[i,q]
   timeclass[i,k,q,1:4] ~ dmulti(timeclass.pis[i,k,q,1:4], n[i,k,q]) #what time interval first observed
  n[i,k,q] ~ dbin(p.dist.sum[i,k,q]*surv[i,q], abund[i,k,q]) #number detected

  
        sigma[i,k,q] <- exp(p.b0 + p.b3*time[i,q] + p.b4[observer[i,q]])
  for (b in 1:bin.dist){
      pbar[i,k,q,b] <- (sigma[i,k,q]^2 * (1- exp(-bb[b+1]^2/(2*sigma[i,k,q]^2))) -
                      sigma[i,k,q]^2 * (1- exp(-bb[b]^2/(2*sigma[i,k,q]^2)))) *
                      2 * 3.141593/area[b]
      pi[i,k,q,b] <- ps[b]*pbar[i,k,q,b]
     
      p.dist.c[i,k,q,b] <- (pi[i,k,q,b] / p.dist.sum[i,k,q]) #conditional probability
     
  } #end b
    
  
      logit(p.avail[i,k,q]) <-  p.avail.b0 + p.avail.b1[song[k]] #time availability?
      pi.pa.sum[i,k,q] <- sum(pi.pa[i,k,q,1:4])
     
    for (t in 1:4){
        pi.pa[i,k,q,t]  <- (p.avail[i,k,q]*p.dist.sum[i,k,q])* pow(1-(p.avail[i,k,q]*p.dist.sum[i,k,q]), (t-1))
        pi.pa.c[i,k,q,t] <- pi.pa[i,k,q,t]/pi.pa.sum[i,k,q]
        } #end t
      } #end i
  } #end k
  
} #end q
 
for (k in 1:n.species){
  psi.b2[k] ~ dunif(-5, 5)
  psi.b3[k] ~ dunif(-5, 5)
  psi.b4[k] ~ dunif(-5, 5)
  psi.b5[k] ~ dunif(-5, 5)
  psi.b6[k] ~ dunif(-5, 5)
  gamma.b2[k] ~ dnorm(0, 1)
  gamma.b3[k] ~ dunif(-5, 5)
  gamma.b4[k] ~ dunif(-5, 5)
  gamma.b5[k] ~ dunif(-5, 5)
  gamma.b6[k] ~ dunif(-5, 5)
  gamma.b7[k] ~ dunif(-5, 5)
  phi.b2[k] ~ dnorm(0,1)
  phi.b3[k] ~ dunif(-5, 5)
  phi.b4[k] ~ dunif(-5, 5)
  phi.b5[k] ~ dunif(-5, 5)
  phi.b6[k] ~ dunif(-5, 5)
  phi.b7[k] ~ dunif(-5, 5)
 
 } #end k
 

 
for (q in 1:n.years){
  for (k in 1:n.species){
    mean.p[k,q] <- mean(p.dist.sum[1:n.sites,k,q])
    spec_abund[k,q] <- sum(abund[1:n.sites,k,q])
   
  } #end k
 
} #end q

for (i in 1:n.sites){
  for (q in 1:n.years){
    for (k in 1:n.species){
      occ[i,k,q] <- (abund[i,k,q] > 0)
    }
  richness[i,q] <- sum(occ[i,1:n.species,q])
  }
}
 
  ## Stable Priors
    p.avail.b1[1] <- 0
    p.avail.b1[2] ~ dunif(-2, 0)
    p.avail.b1[3] ~ dunif(-2, 0)
    p.avail.b0 ~dunif(-3,3)
    psi.b7[1] ~ dunif(-2, 2)
    psi.b7[2] ~ dunif(-2, 2)
    phi.b8[1] ~ dunif(-2, 2)
    phi.b8[2] ~ dunif(-2, 2)
    gamma.b8[1] ~ dunif(-2, 2)
    gamma.b8[2] ~ dunif(-2, 2)
   
  p.b0 ~ dunif(-3,3)
  p.b3 ~ dunif(-5, 5)

  for (y in 1:observers){
  p.b4[y] ~ dunif(-5, 5)
  }
})

[Wei Zhang]

Hi Heather,

The first point I see is that in the first chunk of code the for loop variable q starts from 2, not 1 so that gamma[, ,1] is a matrix of NA. I guess q-1 might be needed in the code. The second point is you have not provided initial values for stochastic nodes S, R, and abund[, , 1]. Fix these and see how it goes. It is a little slow to run all the code on my machine, so I have not tried it. 

Hope this helps.

Best wishes,

Wei