Dynamic Indexing issues translating from JAGs to Nimble

[Dave Pearce]

Hi all,

I am trying to translate a JAGs model to NIMBLE that models the number of vocalizations from acoustic recorders to estimate a species abundance. The reason for moving over to NIMBLE is that my abundance model is underdispersed and I don't have a way of handling it in JAGs, but could through a custom distribution in NIMBLE. However, I am running into some issue moving it over. The model uses a lot of dynamic indexing which I can deal with in part by using custom functions , but part of the model is a ztNB hurdle model and indexes through sites for surveys that have at least one vocalization.

Like so (JAGs code)

 # Site
  for (s in 1:S) {

    # Surveys with Vocalizations
    for (j in 1:J.r[s]) {
   
    # Zero Truncated Negative Binomial
    v[s, A.times[s, j]] ~ dpois((delta[s, A.times[s, j]] * N[s] + omega) * phi[s, A.times[s, j]] * y[s, A.times[s, j]]) T(1, )

    # PPC calls  
    v.pred[s, j] ~ dpois((delta[s, A.times[s, j]] * N[s] + omega) * phi[s, A.times[s, j]] * y[s, A.times[s, j]]) T(1, )
    mu.v[s, j] <- ((delta[s, A.times[s, j]] * N[s] + omega) * phi[s, A.times[s, j]]) / (1 - exp(-1 * ((delta[s, A.times[s, j]] * N[s] + omega) * phi[s, A.times[s, j]])))
    resid.v[s, j] <- pow(pow(v[s, A.times[s, j]], 0.5) - pow(mu.v[s, j], 0.5), 2)
    resid.v.pred[s, j] <- pow(pow(v.pred[s, j], 0.5) - pow(mu.v[s, j], 0.5), 2)
   
    } # End J.r
  } # End S

The issue that I'm dealing with is figuring out a way to index through J.r,  where J.r contains the number of surveys at each acoustic site that contains at least 1 detected vocalization.
> J.r
 [1] 13 10  9 13 14  4 13 14  8 14 13 13 14 11 14  9 13  9 11 13 12  3  8
[24]  2  6  5  2

I know there are some other issues that I may be dealing with, but I'm not sure how to handle a work around for indexing through a vector in NIMBLE.

Suggestions and guidance appreciated.

Thanks!
-Dave

[Perry de Valpine]

Hi Dave,

Thanks for posting a question. Can you please clarify a bit more what you need and are stuck on? Dynamic indexing (i.e. an index that is itself a model node) should work fine in nimble. I'm not clear what aspect of indexing related to J.r is making you stuck. Thanks.

Perry

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 visit https://groups.google.com/d/msgid/nimble-users/caa1531f-a20c-49b1-90d8-57524a671d4bn%40googlegroups.com.

[Dave Pearce]

Hi Perry,

Thank you for responding and apologies for the confusion. It seems that having J.r in the constants fixed that issue. But, now im dealing with some initialization issues that I'm not sure how to handle. Below is the error and the code for the model along with the initial values I've already stated. 

Thank you and apologies, I'm very new to NIMBLE.

-Dave

> model <- nimbleModel(
+   code = acoustic_model,
+   constants = constants_list,  
+   data = Wolfe14.data,
+   inits = inits
+ )
Defining model
Building model
Setting data and initial values
Running calculate on model
  [Note] Any error reports that follow may simply reflect missing values in model variables.
Checking model sizes and dimensions
  [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).
> model$initializeInfo()
  [Note] Missing values (NAs) or non-finite values were found in model variables: lifted_oPdelta_oBs_comma_A_times_oBs_comma_j_cB_cB_times_N_oBs_cB_plus_omega_cP_times_phi_oBs_comma_A_times_oBs_comma_j_cB_cB_times_y_oBs_comma_A_times_oBs_comma_j_cB_cB_L29, K, lifted_v_oBsites_a_oBs_cB_comma_val_times_oBs_comma_j_cB_cB_minus_K_oBs_comma_j_cB_L33.
  [Note] This is not an error, but some or all variables may need to be initialized for certain algorithms to operate properly.
  [Note] For more information on model initialization, see help(modelInitialization).



# Initial Values
inits <- list(
  # Abundance
  N = rep(1, S),        
  beta0 = rnorm(1, 0, 1),
  beta1 = rnorm(1, 0, 1),
  beta2 = rnorm(1, 0, 1),
 
  # Detection
  alpha0 = 0,            
  mu_alpha = 0.5,  
  alpha1 = runif(1, 0, 1),
  alpha2 = rnorm(1, 0, 1),
  alpha3 = rnorm(1, 0, 1),
 
  # Vocalization
  omega = runif(1, 0, 1),  
  gamma0 = log(10),  
 
  # Survey random effect
  Jraneff = rep(0, J_A),
  mu_j = 1,
  tau_j = 1,
 
  # Overdispersion
  mu_phi = 1,  
  tau_phi = 1,
  phi = matrix(1, nrow = S, ncol = J_A)

)


 

 
# -----------------------------
# Model Statement
# -----------------------------
acoustic_model <- nimbleCode({
 
  # ----------------------
  # Abundance Priors
  # ----------------------
 
  # Intercept
  beta0 ~ dnorm(0, 1)
 
  # Covariate effect
  beta1 ~ dnorm(0, 1)  
  beta2 ~ dnorm(0, 1)  
 
  # ------------------------
  # Detection Priors
  # ------------------------
 
  # Intercept
  alpha0 <- logit(mu_alpha)  
  mu_alpha ~ dunif(0, 1)
 
  # True individuals
  alpha1 ~ dunif(0, 1000)  
 
  # Covariate effect
  alpha2 ~ dnorm(0, 1)
  alpha3 ~ dnorm(0, 1)  
 
  # ------------------------
  # Call Rate Priors
  # ------------------------
 
  # False positive rate
  omega  ~ dunif(0, 1000)  
 
  # Intercept
  gamma0 ~ dunif(log(1), log(60))
 
  # Survey random effect  
  mu_j ~ dgamma(0.01, 0.01)
  tau_j ~ dgamma(0.01, 0.01)
  for (j in 1:n_days) {
    Jraneff[j] ~ dnorm(0, tau_j)
  }
 
 
  # Overdispersion
  mu_phi ~ dgamma(0.01, 0.01)    
  tau_phi ~ dgamma(0.01, 0.01)  

  for (s in 1:S) {

    for (j in 1:J_A) {
      phi[s, j] ~ dgamma(mu_phi, tau_phi)
    }
  }
 
  # ------------------------
  # Likelihood and Process Model
  # ------------------------

 
  # Site
  for (s in 1:S) {
   

    # ---------------------------------
    # Abundance Submodel  
    # ---------------------------------
   
    # Poisson
    log(lambda[s]) <- beta0 + beta1 * X_abund[s, 7] + beta2 * X_abund[s, 12]
    N[s] ~ dpois(lambda[s])
   
    # Survey
    for (j in 1:J_A) {
     
    # ---------------------------------
    # Detection Submodel  
    # ---------------------------------
     
    # True Positives + Vegetation Density + Wind
    logit(p_a[s, j]) <- alpha0 + alpha1 * N[s] + alpha2 * X_abund[s, 21] + alpha3 * X_det[j, 2]
     
    # ---------------------------------
    # Call rate Submodel  
    # ---------------------------------
     
    # Intercept + Survey Random Effect
    log(delta[s, j]) <- gamma0 + Jraneff[j]
     
    # ---------------------------------
    # Observations
    # ---------------------------------
    y[s, j] ~ dbin(p_a[s, j], 1)
     
    # ---------------------------------
    # True Positives
    # ---------------------------------
    tp[s, j] <- delta[s, j] * N[s] / (delta[s, j] * N[s] + omega)
     
   
    } # End J
   
    # ---------------------------------
    # Vocalizations  
    # ---------------------------------

    # Surveys with Vocalizations
    for (j in 1:J_r[s]) {

      # Zero Truncated Negative Binomial

      v[s, A_times[s, j]] ~ T(dpois((delta[s, A_times[s, j]] * N[s] + omega) * phi[s, A_times[s, j]] * y[s, A_times[s, j]]), 1, )

      }# End J_R
   
  } # End S
 
  # ------------------------
  # Manual Validation
  # ------------------------

  for (s in 1:S_val) {

    for (j in 1:J_val[s]) {
     
      K[s, j] ~ dbin(tp[sites_a[s], j], v[sites_a[s], val_times[s, j]])
      k[s, val_times[s, j]] ~ dhyper_nimble(K[s, j], v[sites_a[s], val_times[s, j]] - K[s, j], n[s, val_times[s, j]])
    } # End J
  } # End S


 
  # -------------------------------------------
  # Derive Parameters
  # -------------------------------------------
 
  # Abundance
  N_tot <- sum(N[1:S])
})
# ---------------------------- End Model ----------------------------

 

[Perry de Valpine]

Hi Dave,

Glad we can help get you started. First, as the message states, it is not necessarily an error to not initialize your model fully.

The nodes listed as not initialized are some lifted nodes (see our User Manual at r-nimble.org if you need background on lifted nodes -- basically nodes nimble has inserted to implement your model) and K.

It looks like K is a posterior predictive node (from your code comments), so it shouldn't be a problem that it's not initialized. 

The first lifted node name is a text conversion of "(delta[s, A_times[s, j]] * N[s] + omega) * phi[s, A_times[s, j]] * y[s, A_times[s, j]])". The second is a text conversion of "v[sites_a[s], val_times[s, j]] - K[s, j], n[s, val_times[s, j]]". The second one depends on K[s, j], so it makes sense that if K[s, j] is not initialized neither is something that depends on it.

I'm not following in enough detail to see what is up with the first one, but you can start inspecting the model as a way to find out. 

E.g.

model$calculate("v") #calculate all of it

model$calculate("v[1, ]") #calcualte part of it, etc., to isolate where it is using values that are not initialized

model$getParam("v[1, 1]", "mean") # inspect directly its mean parameter value

model$N # and/or look at variables that are used in the calculation

model$phi #etc.

HTH

Perry

To view this discussion visit https://groups.google.com/d/msgid/nimble-users/c037483d-5b42-4d9f-b5d8-7ca31fa12e6en%40googlegroups.com.