Error in model due to Nimble function: "Dimension of x is 0, which does not match its usage in(...)"

[Jaume Adria Badia Boher]

Hi all!

I am a new NIMBLE user with experience in JAGS and working with bayesian ecological models. I am struggling to create a function (nimbleFunction), and any help would be much appreciated.

My function is a part of a spatial component of a population model. The aim of the function is to use some gamma-distributed generated distances (Dist) to find which the closest point to them is (Distmat) based on a departure point (birthterr). 

sortbyterr <- nimbleFunction(

  run = function(Distmat = double(2), birthterr = double(1), Dist = double(1), nterr = double()){

    returnType(double(1))

    n1 <- length(Dist)

    recterr <- numeric(length=n1, init=FALSE)

    

    for(i in 1:n1){

      recterr[i] <- which((abs(Distmat[birthterr[i], 1:nterr] - Dist[i])) == min(abs(Distmat[birthterr[i], 1:nterr] - Dist[i])))

    }

    

    return(recterr)

  })

You can easily replicate the function by creating these objects:

Dist <- rgamma(20, 0.8, 1.1) #Movement distances

Distmat <- matrix(nrow = 20, ncol= 40, rgamma(20*40,  0.8, 1.1)) #Matrix of distances between all departure and arrival points

nterr <- dim(Distmat)[2] #Number of departure/arrival points

birthterr <- sample(1:20, 20, replace=T)  #Departure rows

The function works well outside NIMBLE. However, the model delivers an error, as if the function is not creating the output vector properly:

Error in genVarInfo3() : 

  Dimension of recterr is 0, which does not match its usage in 'phiT[i, t] <- expit(muST[sexvec[i], pcatsurv[recterr[i]]] + epsST[sexvec[i], t, pcatsurv[recterr[i]]])'.

I don't fully understand the reason for this error, but it might be related with the following. In the model, the Dist vector is actually entered as data and is a combination of observed distance values and NAs for those individuals in which movement has not been observed. However, the model is asked to infer distance values in all individuals, and run the function after that:

 for(i in 1:nind){

    Dist[i] ~ dgamma(sh, ra)

  }

  

  sh ~ dgamma(1,1)

  ra ~ dgamma(1,1)

recterr <- sortbyterr(Distmat, birthterr, Dist, nterr)

As I understand, the Dist vector used in the function should no longer contain NAs, as all unknown values should have been simulated at that time. Hence, I do not actually understand the origin of my error and no solution comes to my mind.

Any comment would be more than appreciated! Please let me know if something was not clear or if you need more context info. 

Thanks in advance!

[Perry de Valpine]

Hi Jaume,

Good to see this application and thanks for reaching out for support.

I see a couple of issues going on.

In model code (but not nimbleFunction code), nimble requires square brackets to indicate non-scalars.  So you will need something like:

recterr[1:nterr] <- sortbyterr(Distmat[1:n, 1:nterr], birthterr[1:n], Dist[1:n], nterr)

Sometimes you can use simply recterr[] (etc, with Distmat[,] for a matrix), if either the sizes can be determined by nimble from other model declarations or are provided in the dimensions argument to nimbleModel.

Since you said you have experience with JAGS, let me point you to our guide to converting from JAGS to nimble, in case you haven't seen it yet.  It is linked from our documentation page.

The other issue is that there is a glitch in your nimbleFunction, which is actually a bug in our compilation system.  We encourage trying to compile outside of the model to check a nimbleFunction:

Csortbyterr <- compileNimble(sortbyterr)

You can then use Csortbyterr as a function from R.  In this case, there is an error from compileNimble, and I see it is an issue with how our compiler tries to support which().  I'll file a GitHub issue for us to fix that.

However, this gives me an excuse to comment on efficiency in R vs C++.  In R, users learn to vectorize when possible because for-loops are often slower than vectorized operations such as your `recterr[i] <- which(...)` line.  However, in C++, for loops are fast, and although we (try to) support vectorized operations, they are sometimes faster but sometimes slower than for loops.  In this case, the min() operation must internally iterate (run a for loop) over all elements, and then the which operation must run another for loop.  So the extent 1:nterr is being looped over (at least) twice, once to find the min, then again to find the which.min.  The following should, I suspect, be more efficient, although perhaps not meaningfully so unless nterr is very large.

sortbyterr <- nimbleFunction(
  run = function(Distmat = double(2), birthterr = double(1), Dist = double(1), nterr = double()){
    returnType(double(1))
    n1 <- length(Dist)
    recterr <- numeric(length=n1, init=FALSE)
   
    for(i in 1:n1) {

      min_val <- Inf
      min_ind <- 1
      for(j in 1:nterr) {
        this_val <- abs(Distmat[birthterr[i], j] - Dist[i])
        if(this_val < min_val) {
          min_val <- this_val
          min_ind <- j
        }
      }
      recterr[i] <-  min_ind
    }
   
    return(recterr)
  })

Csortbyterr <- compileNimble(sortbyterr) # compilation seems to work, but please test that the behaves correctly.

This does a single loop to find the min and record its index.

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/051af5e4-758a-487f-9ff4-582ec968767en%40googlegroups.com.

[Jaume Adria Badia Boher]

Hi Perry!

Thank you so much for your feedback. It was really useful and the jags-to-nimble guide is helping me a lot. 

I am translating all the models I'm currently using in my research from JAGS to Nimble, and I realized it is way better to start from the simplest models and add complexity step by step to deal with all errors one by one. Please, if you guys do not mind, I will ask one further question as there is one recurrent error I did not find any information on how to address. Any feedback would be very appreciated.

My aim is to model survival and territorial recruitment in a wild animal population using a multistate model. My base model used constant probabilities for all parameters and ran perfectly. A next step was to add an age-by-sex structure ("ages_sex") to survival and recruitment by means of adding a matrix that indicates age categories for males and females from 1 to 5. When creating the model, it complained about the variables structured by ages_sex not being fully initialized, but the model ran smoothly anyways. Now, I'm stuck when trying to model territorial recapture (pT) as it was in the field, with intervals with no sampling inserted between sampling periods. I am trying to model it as a constant probability, and a further step will be to model it as a random effect. You will find the code below.

Apparently the model does not like this way of parameterizing recapture. Models ran with equivalent structures in JAGS, but in Nimble, when trying to run the model, I have dozens - hundreds of warning messages about logProbs of data nodes from y (observation matrix) being NA or NaN. Perhaps I am missing some coding differences in Nimble, although I tried to follow all guidelines.

Many thans again for your help so far! If you have time to check the code and find to need any more information, do not doubt to ask.

Have a nice day,

Jaume A.

nimp4 <- nimbleCode({

  

  #Priors and constraints

  survT ~ dunif(0,1)

  mean.pNT ~ dunif(0,1)

  mean.rNT ~ dunif(0,1)

  mean.rT ~ dunif(0,1)

  mean.pT ~ dunif(0,1)  

  for(u in 1:10){

    survNT[u] ~ dunif(0,1)

    gam[u] ~ dunif(0,1)

  }

   

  #Modelling territorial recapture: 0s for non-sampling intervals, constant probability for sampling intervals

  pT[1] <- 0

  pT[2] <- mean.pT

  pT[3] <- 0

  pT[4] <- mean.pT

  pT[5] <- 0

  pT[6] <- mean.pT

  pT[7] <- 0

  pT[8] <- mean.pT

  pT[9] <- 0

  pT[10] <- mean.pT

  pT[11] <- 0

  pT[12] <- mean.pT

  pT[13] <- 0

  pT[14] <- mean.pT

  pT[15] <- 0

  pT[16] <- mean.pT

  pT[17] <- 0

  pT[18] <- mean.pT

  pT[19] <- 0

  pT[20] <- mean.pT

  pT[21] <- 0

  pT[22] <- mean.pT

  pT[23] <- 0

  pT[24] <- mean.pT

  

    #Define state-transition and observation matrices

  

  #STATES:

  #1: ALIVE, NON-TERRITORIAL

  #2: ALIVE, TERRITORIAL

  #3: DEAD, NON-TERRITORIAL

  #4: DEAD, TERRITORIAL

  #5: LONG-DEAD

  

  #OBSERVATIONS:

  #1: ALIVE, NT

  #2: ALIVE, T

  #3: DEAD, NT

  #4: DEAD, T

  #5: UNOBSERVED

  

  

  for(i in 1:nind){

    for(t in f[i]:(n.occasions-1)){  #f is the vector of first encounters

      phiNT[i,t] <- survNT[ages_sex[i,t]]

      gamma[i,t] <- gam[ages_sex[i,t]]

      

      #Survival and recruitment are age-dependent    

#Ages_sex is the sex-by-age matrix, with 1:5 for males and 6:10 for females. 

#we changed NAs by 0s between first occasion and first encounter since we already had some issues with NAs in 

#the state matrix initial values (z). The model seemed to accept it.

      

      

      #State transition matrix

      ps[1,1,i,t] <- phiNT[i,t]

      ps[1,2,i,t] <- phiNT[i,t]*gamma[i,t]

      ps[1,3,i,t] <- (1-phiNT[i,t])

      ps[1,4,i,t] <- (1-phiNT[i,t])*gamma[i,t]

      ps[1,5,i,t] <- (1-phiNT[i,t])

      ps[2,1,i,t] <- 0

      ps[2,2,i,t] <- survT

      ps[2,3,i,t] <- 0

      ps[2,4,i,t] <- (1-survT)

      ps[2,5,i,t] <- (1-survT)

      ps[3,1,i,t] <- 0

      ps[3,2,i,t] <- 0

      ps[3,3,i,t] <- 0

      ps[3,4,i,t] <- 0

      ps[3,5,i,t] <- 1

      ps[4,1,i,t] <- 0

      ps[4,2,i,t] <- 0

      ps[4,3,i,t] <- 0

      ps[4,4,i,t] <- 0

      ps[4,5,i,t] <- 1

      ps[5,1,i,t] <- 0

      ps[5,2,i,t] <- 0

      ps[5,3,i,t] <- 0

      ps[5,4,i,t] <- 0

      ps[5,5,i,t] <- 1

      

    } #t

  } #i

  

  

  for(t in f[i]:(n.occasions-1)){

  

  #Observation probabilities

  po[1,1,t] <- mean.pNT

  po[1,2,t] <- 0

  po[1,3,t] <- 0

  po[1,4,t] <- 0

  po[1,5,t] <- (1-mean.pNT)

  po[2,1,t] <- 0

  po[2,2,t] <- pT[t]

  po[2,3,t] <- 0

  po[2,4,t] <- 0

  po[2,5,t] <- (1-pT[t])

  po[3,1,t] <- 0

  po[3,2,t] <- 0

  po[3,3,t] <- mean.rNT

  po[3,4,t] <- 0

  po[3,5,t] <- (1-mean.rNT)

  po[4,1,t] <- 0

  po[4,2,t] <- 0

  po[4,3,t] <- 0

  po[4,4,t] <- mean.rT

  po[4,5,t] <- (1-mean.rT)

  po[5,1,t] <- 0

  po[5,2,t] <- 0

  po[5,3,t] <- 0

  po[5,4,t] <- 0

  po[5,5,t] <- 1

  

  }

  

  # Likelihood

  for (i in 1:nind){

    # Define latent state at first capture

    z[i,f[i]] <- y[i,f[i]]

    for (t in (f[i]+1):n.occasions){

      # State process: draw S(t) given S(t-1)

      z[i,t] ~ dcat(ps[z[i,t-1], 1:5, i, t-1]) 

      # Observation process: draw O(t) given S(t)

      y[i,t] ~ dcat(po[z[i,t], 1:5, t]) 

    } #t

  } #i

})

ages_sex_nimble <- ages_sex[,1:24]

ages_sex_nimble[is.na(ages_sex_nimble)] <- 0

my.data <- list(y=CMRmatrix)  

my.constants <- list(n.occasions=dim(CMRmatrix)[2], nind=dim(CMRmatrix)[1], f=f, ages_sex = ages_sex_nimble)

z2 <- my.inits4(CMRmatrix, f, 5)

z2[is.na(z2)] <- 1

my.inits <- function(){list(survNT=runif(10,0,1),

                            mean.pNT=runif(1,0,1),

                            mean.pT=runif(1,0,1),

                            survT=runif(1,0,1),

                            gam=runif(10,0,1),

                            mean.rNT=runif(1,0,1),

                            mean.rT=runif(1,0,1),

                            z=z2)}

El dia dijous, 1 de juliol de 2021 a les 16:03:23 UTC+2, pdevalpine va escriure:

[Perry de Valpine]

Hi Jaume,

I'm glad that was helpful.

The warnings about variables not being fully initialized are just that, warnings, not errors.  A difference from JAGS is that you do not need all inits and data to be set up when you call nimbleModel, although many find it convenient to do so.  By default, nimbleModel tries calculating the entire model and lets you know if anything is not initialized, in case you thought it was fully initialized.  If you don't want it to do this (e.g. after you have your workflow working smoothly), you can do 'calculate=FALSE' in the nimbleModel call.

If you run MCMC with incomplete initial values, it might work or it might not.  You may get those warnings from a full calculation at the start and then it might run fine.  The MCMC will try to initialize a node that is NA by drawing from its prior.  In some models, that value can make something that depends on it invalid, with a log probability of -Inf, if the value of one node determines the valid range of values for another.  Sometimes the MCMC will successfully mix its way out of invalid initial states, sometimes not.  In summary, if you get a slew of warnings, and then they stop, and then you get reasonable MCMC results, you are probably fine.  Either way, many people find it helpful to achieve full initialization.  The nodes that need to be initialized are stochastic nodes.  A related issue is that even if drawing initial values from priors does not give warnings or errors, those are often not very good initial values.  

From reading your code, I am not quickly seeing what nodes are not fully initialized.  You can use the model object to debug the situation.  For example, if the messages are telling you that survNT is not fully initialized, you can do:

model$survNT # look at the values

model$calculate('survNT') # look at their log probability sum

model$calculate('survNT[3]') # one by one if needed

(assuming "model" is from model <- nimbleModel(<your inputs>)

If the problem is further down and unclear where, you can do something like this:

for(yNode in model$expandNodeNames('y')) {

  logProb <- model$calculate(yNode)

  cat(paste("logProb for", yNode, "=", logProb))

}

If you determine a specific yNode has an invalid logProb, you can check its inputs:

model$po[<indices as needed>]

model$z

and so on.

I would not say nimble is unhappy with writing a model this way, but there is something going on for you to track down.

Also you may be interested in nimbleEcology, which gives marginal distributions for hidden Markov models, which you could use here.  Some of our workshop materials cover this.  A recent short workshop covering this for capture-recapture is here.

If you still can't track it down, I think we'll need a fully reproducible example to help further.  Good luck.

-Perry

  

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/0a6157fb-84ab-465f-a5e5-81dd09a6a0efn%40googlegroups.com.