nimStep function error

[Mattia Falaschi]

Good morning everyone, 

I am trying to run a ricker model to estimate abundance and detection probability of newts.

The model was written for Jags, but since I have some issues (and because nimble is much faster!) I am trying to run it in nimble.

This model contains a step function (nimStep in nimble). When I try to define the model, I get this error: 

"defining model...

Error in processBUGScode(recurseCode, nextContextID, lineNumber = lineNumber,  : 

  Error: * not allowed in BUGS code in +(1 - nimStep(N[i, t - 1] - 0.01)) * exp(rho + S[i] + T[t])"

Do you know if there is some issue with the nimStep function or if I am doing something wrong?

If you need more to replicate the error, let me know.

Thank you!

Mattia

[Mattia Falaschi]

Here is the code of the model

[Daniel Turek]

I think it's just the way you have the line breaks separating single declarations in your code file, specifically, line 13:

      L[i,t] <- nimStep(N[i,t-1]-0.01) * exp(log(N[i,t-1]) + rho - eta*N[i,t-1] + S[i] + T[t])

then line 14 (which is intended to be a continuation of line 13):
      +

and then, after some comments, the intended continuation of line 13, which is line 17:

        (1 - nimStep(N[i,t-1]-0.01)) * exp(rho + S[i] + T[t])

Rather, put the "+" line at the end of line 13, and I think it should work.

On Mon, Feb 15, 2021 at 6:31 AM Mattia Falaschi <matt...@gmail.com> wrote:

Here is the code of the model

--
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/907bfea3-bef7-4f4c-92b6-4ef690d958e5n%40googlegroups.com.

[Mattia Falaschi]

Thank you very much, I changed the script as you suggested and now it works!

Well, new errors appeared, but they seem to be easier to solve.

Thank you again.

[Mattia Falaschi]

Hello again, 

I am having more troubles with the same model.

I hope you can help me in solving these issues.

Here is the model. I also attach the data and the script to run it.

model.ricker.nb <- nimbleCode({

  for (i in 1:nSite){

    N[i,1] ~ dnegbin(p.nb[i], alpha)

    p.nb[i]<-alpha/(alpha+lambda[i])

    lambda[i] <- exp(b0)

    

    for (t in 2:Tf[i]){

      N[i,t] ~ dpois(L[i,t])

      L[i,t] <- step(N[i,t-1]-0.01) * exp(log(N[i,t-1]) + rho - eta*N[i,t-1] + S[i] + T[t]) +

               (1 - step(N[i,t-1]-0.01)) * exp(rho + S[i] + T[t])

    }

    

    for (t in T0[i]:Tf[i]){

      y[i,t] ~ dbin(p,N[i,t])

    }

  }

  

  # Priors

  for(i in 1:nSite){

    S[i] ~ dnorm(0,tau.S)

  }

  for(t in 2:nYear){

    T[t] ~ dnorm(0,tau.T)

  }

  

  tau.S <- pow(sigma.S,-2)

  sigma.S ~ T(dt(0,1,1),0,)

  

  tau.T <- pow(sigma.T,-2)

  sigma.T ~ T(dt(0,1,1),0,)

  

  b0 ~ dnorm(0,0.01)

 

  p ~ dunif(0,1)

  

  eta ~ dnorm(0,0.01)

  rho ~ dnorm(0,0.01)

  alpha ~ dunif(0,50)

})

When running the "nimbleModel" function, I get some warnings stating

"In dbinom(model$y[getNodeFunctionIndexedInfo(INDEXEDNODEINFO_,  ... : NaN produced)"

Then (I do not know if it is related) I get many times this message when running the model:

"Warning: slice sampler reached maximum number of contractions for 'N[8, 2]'. Current parameter value is 1."

The problem is that I do not understand the meaning of these warnings.

I have many NA in the dataset, but I should have avoided them by adding indexing T0 and Tf into the likelihood.

Thank you again for your patience!

[Daniel Turek]

Mattia, thanks for your persistence and use of NIMBLE.  I tried running your model, and a few comments:

Noting that each lambda[i] is defined the same as exp(b0), so also each p.nb[i] term is exactly the same.  I'm not sure if you meant it this way, or perhaps you're working up to a more complex model?  But if not, the code could be simplified.

I took a look at the errors, and I think it just comes down to choices of initial values.  I made a few changes:

The initial value for the p parameter, you had "rnorm(1,0,1)", but as p is a probability, this was often giving negative values for p, which led to an invalid model state. I changed the initial value for p to "p = runif(1, 0, 1)"

The initial value for eta, you had "rnorm(1, 0, 1)", and while valid, in many cases this was "too large" of a value, and the term from your model:

exp(log(N[i,t-1]) + rho - eta*N[i,t-1] + S[i] + T[t])

when eta was say -1, and a value of N is say 10, this would lead to L being on the order exp(10), which was way too large, and was causing a density evaluation for N (since N[i,t] ~ dpois(L[i,t])) of -Inf.  So, we need a more "neutral" initial value for eta.  I changed the initial value simply to 0, and we can let the MCMC explore from there.  If you truly want a "random" initial value for eta, then I would recommend something much closer to 0, such as "rnorm(1, 0, 0.001)"

The same problem was happening in the initial values for both S and T.  Their initial values, say for S:

S = rnorm(nrow(counts),0,pow(rtruncnorm(1,a=0,,mean=0,sd=1),-2)),

The problem is, when the truncnorm() function returns a small value, near to 0, say 0.1, then pow(0.1, -2) = 100 is the standard deviation for drawing the initial values for S (and same happening for T), and therefore the initial values of S (and T) are often on the order of 100, or 200, then these are also appear in the term:

L[i,t] <- step(N[i,t-1]-0.01) * exp(log(N[i,t-1]) + rho - eta*N[i,t-1] + S[i] + T[t]) + ...

and the exp(S), or exp(T) again gives a huge number, and would cause the density evaluation (the same problem as before) for N to be -Inf.  So, I changed the initial values for S and T to have smaller standard deviations, to give smaller numbers for S and T:

S = rnorm(nrow(counts),0,1)

T = rnorm(ncol(counts),0,1)

All told, the new initial values list looks like:

inits <-list(N = (counts+1), b0 = rnorm(1,0,1),
             p = runif(1, 0, 1), ##rnorm(1,0,1),
             alpha = runif(1,0,50),
             eta = 0, ##rnorm(1,0,1),
             rho = rnorm(1,0,1),
             S = rnorm(nrow(counts),0,1), ##pow(rtruncnorm(1,a=0,,mean=0,sd=1),-2)),
             T = rnorm(ncol(counts),0,1), ##pow(rtruncnorm(1,a=0,,mean=0,sd=1),-2)),
             sigma.S = rtruncnorm(1,a=0,,mean=0,sd=1), sigma.T = rtruncnorm(1,a=0,,mean=0,sd=1),
             tau.S = pow(rtruncnorm(1,a=0,,mean=0,sd=1),-2), tau.T=pow(rtruncnorm(1,a=0,,mean=0,sd=1),-2))

And this seemed to work for me, and not produce any warnings or errors.

To diagnose these problems, you can look at the variables in the model, such as:

Rmodel$T
Rmodel$S
Rmodel$p

And equally useful, you can look at the density evaluations (log-probabilities) of stochastic model variables as:

Rmodel$logProb_y
Rmodel$logProb_N

So, that way you can figure out where the -Infs, or NaNs are coming from.

Hope this helps.  Keep at it.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/a11aa538-8187-44e5-846c-39429e8682f3n%40googlegroups.com.

[Mattia Falaschi]

Thank you very much Daniel for your enlightening reply!

I always have issues with inits, but I am glad that the model itself is fine.

I modified the inits as you suggested and the model runs fine, with no more errors/warnings, even after adding some covariates.

Regarding b0: yes, I am working on a more complex model to estimate the effect of some covariates on the abundance of newts, this is the reason why p.nb is defined for each site.

Again, thank you very much for the help.

Mattia

[Mattia Falaschi]

Here I am again with another question regarding the model above.

First of all, thank you Daniel for helping me until this point and thank you for your patience.

After further modifying the model, it seems that I need your help again.

I added a parameter of colonization probability (gamma), which I assumed constant and will not be related to any covariate.

My colleague ran the model in JAGS and obtained estimates of the gamma parameter (low value, <0.1). However, when I run this in nimble, it seems that gamma returns the prior (uniform between 0 and 1).

Here is the updated script of the model. 

I attach again the data and the script to run it.

Thank you in advance for your help!

Mattia

model:

model.ricker.nb <- nimbleCode({

  for (i in 1:nSite){

    N[i,1] ~ dnegbin(p.nb[i], alpha)

    p.nb[i]<-alpha/(alpha+lambda[i])

    lambda[i] <- exp(b0)

    

    for (t in 2:Tf[i]){

      N[i,t] ~ dpois(L[i,t])

      log(L[i,t]) <- step(N[i,t-1]-0.01) * (log(N[i,t-1]) + rho - eta*N[i,t-1] + S[i] + T[t]) +

        

        (1 - step(N[i,t-1]-0.01)) * log(-log(1-gamma))

      }

    

    for (t in T0[i]:Tf[i]){

      y[i,t] ~ dbin(p,N[i,t])

    }

  }

  

  # Priors

  for(i in 1:nSite){

    S[i] ~ dnorm(0,tau.S)

  }

  for(t in 2:nYear){

    T[t] ~ dnorm(0,tau.T)

  }

  

  tau.S <- pow(sigma.S,-2)

  sigma.S ~ T(dt(0,1,1),0,)

  

  tau.T <- pow(sigma.T,-2)

  sigma.T ~ T(dt(0,1,1),0,)

  

  b0 ~ dnorm(0,0.01)

  p ~ dunif(0,1)

  gamma ~ dunif(0,1)

  

  eta ~ dnorm(0,0.01)

  rho ~ dnorm(0,0.01)

  alpha ~ dunif(0,50)

})

[Daniel Turek]

Mattia, this was subtle to track down.  The root cause comes down to this behaviour in R (which NIMBLE replicates exactly):

> 0 * -Inf
[1] NaN

How does this relate to your problem?  That's a good question.

Briefly, when a value for any N[i,t-1] is proposed to transition to 0, NIMBLE's MCMC calculates the model log-densities associated with that transition.  With N[i,t-1] = 0, the corresponding log(L[i,t]) is calculated using:

log(L[i,t]) <- step(N[i,t-1]-0.01) * (log(N[i,t-1]) + .....) + ....

The step(...) term evaluates to 0, and the (log(N[i,t-1] + ...) evaluates to -Inf, and therefore, regardless of being multiplied by 0, the entire log(L[i,t]) ends up at NaN, and also therefore L[i,t] is NaN.  Finally, the resulting log-density in the model of N[i,t], and note that it doesn't matter what N[i,t] contains, comes from

N[i,t] ~ dpois(L[i,t])

since L[i,t] is NaN, the log density of N[i,t] is NaN, and the proposal of N[i,t-1] to transition to 0 is rejected.  Therefore, the N values *never* transition to 0, and the gamma term is *never* used (or more accurately, it's always multiplied by zero), and therefore posterior(gamma) = prior(gamma), as you observed.

It sounds like JAGS does some numerical calculations in a different way - which aren't entirely consistent with R.  In this case, unfortunately, that's not how NIMBLE operates.

There are probably many different ways to fix this, but the quickest I came up with involves a quick user-defined distribution.  I hope this addition of a custom-written distribution doesn't end up seeming intimidating, but rather perhaps will showcase the simplicity and functionality of NIMBLE.

I think the simple changes to your code, below, will get you going with your model that includes the gamma colonization parameter:

calculateL <- nimbleFunction(
    run = function(N = double(), rho = double(), eta = double(), S = double(), T = double(), gamma = double()) {
        returnType(double())
        if(N == 0) return(log(-log(1-gamma)))
        if(N > 0) return(log(N) + rho - eta*N + S + T)
    }
)

model.ricker.nb <- nimbleCode({
    for (i in 1:nSite){
        N[i,1] ~ dnegbin(p.nb[i], alpha)
        p.nb[i]<-alpha/(alpha+lambda[i])
        lambda[i] <- exp(b0)
        for (t in 2:Tf[i]){
            N[i,t] ~ dpois(L[i,t])

            log(L[i,t]) <- calculateL(N[i,t-1], rho, eta, S[i], T[t], gamma)    ## this line changed
       ## etc......

If you have time, then let us know if this works for you, and also if it agrees with JAGS.

Cheers,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/81ea3f48-35f6-44dc-833d-d555c870a366n%40googlegroups.com.

[Mattia Falaschi]

Again, thank you very much Daniel!

I tried a quick run with 80 000 iteration, and while parameter is still not at convergence, something changed and it seems at least to modify the prior!

I have two additional points:

1) Since gamma is a colonization probability, I need to calculate it when N[i,t-1] == 0 and not when N == 0.

It seems that this is considered in this line

log(L[i,t]) <- calculateL(N[i,t-1], rho, eta, S[i], T[t], gamma)

So that in the calculateL function the N value is taken from t-1 so that the if conditions are already verified at t-1.

Is this correct, or do I need to modify something within the calculateL function?

2) I need to add some covariates for N[i,1] and N[i,t]. 

Assuming that I have three covariates and that coefficients are b1, b2, and b3, if I understood correctly what the custom function is doing, the script should be modified like this:

 calculateL <- nimbleFunction(

    run = function(N = double(), rho = double(), eta = double(), S = double(), T = double(), gamma = double(),

                   b1 = double(), b2 = double(),  b3 = double()) {

        returnType(double())
        if(N == 0) return(log(-log(1-gamma)))

        if(N > 0) return(log(N) + rho - eta*N + S + T + b1*var1[i,t] + b2*var2[i,t] + b3*var3[i,t]) ## Covariates are considered only when N[i,t-1] > 0 

    }
)

model.ricker.nb <- nimbleCode({
    for (i in 1:nSite){
        N[i,1] ~ dnegbin(p.nb[i], alpha)
        p.nb[i]<-alpha/(alpha+lambda[i])

        lambda[i] <- exp(b0 + b1*var1[i,1] + b2*var2[i,1] + b3*var3[i,1]) ## Same as older model

        for (t in 2:Tf[i]){
            N[i,t] ~ dpois(L[i,t])

            log(L[i,t]) <- calculateL(N[i,t-1], rho, eta, S[i], T[t], gamma, b1, b2, b3)    ## Added coefficients here
       ## etc......

Mattia

[Daniel Turek]

You're welcome Mattia.  

1) Yes, the value of N used in the calculateL() function, when it's doing the calculations for log(L[i,t]), the function will use the values of N = N[i,t-1].  I believe that's the behavior that your original code used, and that you want.

2) With the covariates, there are just a few differences from the code you sent.  Mainly, the variables var1, var2, and var3 will also have to be passed into the calculateL function as arguments.  See changes below, in blue:

 calculateL <- nimbleFunction(
    run = function(N = double(), rho = double(), eta = double(), S = double(), T = double(), gamma = double(),

                   b1 = double(), b2 = double(), b3 = double(), var1 = double(), var2 = double(), var3 = double()) {

        returnType(double())
        if(N == 0) return(log(-log(1-gamma)))

        if(N > 0) return(log(N) + rho - eta*N + S + T + b1*var1 + b2*var2 + b3*var3) ## no indexing on var1, var2, var3
    }
)

model.ricker.nb <- nimbleCode({

    ## add priors for b1, b2, b3

    for (i in 1:nSite){
        N[i,1] ~ dnegbin(p.nb[i], alpha)
        p.nb[i]<-alpha/(alpha+lambda[i])

        lambda[i] <- exp(b0 + b1*var1[i,1] + b2*var2[i,1] + b3*var3[i,1]) ## this line is ok

        for (t in 2:Tf[i]){
            N[i,t] ~ dpois(L[i,t])

            log(L[i,t]) <- calculateL(N[i,t-1], rho, eta, S[i], T[t], gamma, b1, b2, b3, var1[i,t], var2[i,t], var3[i,t])   ## changes here
       ## etc......

##

## and, in the constants list used for building the model,

## now also include var1, var2, and var3

##

This should work fine, hopefully if there aren't any typos.

Keep at it,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/0fcce069-2651-40d4-97b6-e80f520e7272n%40googlegroups.com.

[Mattia Falaschi]

1) Yes, this is exactly what I wanted

2) The model works and estimated coefficients agree with the ones of jags. Thank you very much for your help, it is really appreciated!

Now I can also understand the differences in iterations needed to reach convergence between jags and nimble: the NaN issue seems to affect this, since with jags I needed 100-200k iterations, while with nimble I needed 1-2 millions. With the fix you suggested, now the number of iterations needed is the same as jags.

All the best, 

Mattia