Error NAs were detected in model variables

142 views
Skip to first unread message

Shrinivas Shirke

unread,
Dec 12, 2023, 12:28:33 AM12/12/23
to nimble-users
I am trying to run the code wherein O I X~poisson , X~ negative binomial. Following is the output log,
Checking model calculations [Note] NAs were detected in model variables: x, logProb_x, mu, logProb_O.

running chain 1... warning: problem initializing stochastic node x[1]: logProb is -Inf. warning: problem initializing stochastic node x[2]: logProb is -Inf. warning: problem initializing stochastic node x[3]: logProb is -Inf. warning: problem initializing stochastic node x[4]: logProb is -Inf. warning: problem initializing stochastic node x[6]: logProb is -Inf.
The code is mentioned below for kind reference

# Reference : https://gkonstantinoudis.github.io/nimble/RmarkSLC.html
# ScotLip data downloaded from https://geodacenter.github.io/data-and-lab/scotlip/
ScotLip <- read_sf("ScotLip.shp")
ScotLip
# Negative Binomial distribution with mean lambda and rate parametrization
dmynegbinom <- nimbleFunction(
  run = function(x = double(0),
                 rate = double(0),
                 lambda = double(0),
                 log = integer(0, default = 1)) {
    returnType(double(0))
    # compute binomial coefficient
    lchoose <- lgamma(x+rate) - lgamma(rate) - lgamma(x)
    prob=rate/(rate+lambda)
    # binomial density function
    logProb <- lchoose + x * log(1-prob) + (rate) * log(prob)
    if(log) return(logProb)
    else return(exp(logProb))
  })


rmynegbinom <- nimbleFunction(
  run = function(n= integer(0, default = 1),
                 rate = double(0),
                 lambda = double(0)) {
    returnType(double(0))
     y<-0
      pmf=(rate/(rate+lambda))^rate;
      cdf=pmf;
      j=0;
      u=runif(1);
      while(u>=cdf){
        pmf<-((j+1+rate)/(j+1))*(lambda/(lambda+rate))*pmf
        cdf=cdf+pmf;
        j=j+1
      }
      y=j
     return(y)
  })

assign('dmynegbinom', dmynegbinom, .GlobalEnv)
assign('rmynegbinom', rmynegbinom, .GlobalEnv)

ICARCode <- nimbleCode(
  {
    for (i in 1:N){
      O[i]~dpois(mu[i])
      log(mu[i])<-(x[i]>0)*x[i]*beta+(x[i]==0)*beta
      x[i] ~ dmynegbinom(lambda[i] , r)                              
      log(lambda[i]) <- log(E[i]) + alpha + phi[i]
      SIR[i] <- exp(alpha + phi[i])                    
      resSIR[i] <- exp(phi[i])                          
     
    }
   
 
    phi[1:N] ~ dcar_normal(adj[1:L], weights[1:L], num[1:N], tau.phi, zero_mean = 1)
    # the zero_mean is to impose sum to zero constrains in case you have an intercept in the model
   
    # Priors:
    alpha ~ dflat()
    beta ~ dflat()
    # vague uniform prior
    overallSIR <- exp(alpha)                             # overall SIR across study region
    r~dgamma(1,0.01)
    tau.phi ~ dgamma(1, 0.01)                            # precision of the ICAR component
    sigma2.phi <- 1/tau.phi                              # variance of the ICAR component
   
  }
 
)
N <- dim(ScotLip)[1]

# Format the data for NIMBLE in a list
ScotLipdata = list(
  O = ScotLip$CANCER                 # observed nb of cases
)        

ScotLipConsts <-list(        
  N = N,                             # nb of wards  
  E = ScotLip$CEXP,                  # expected number of cases
 
  # and the elements of the neighboring matrix:
  L = length(nbWB_A$weights),        
  adj = nbWB_A$adj,                            
  num = nbWB_A$num,
  weights = nbWB_A$weights
 
)
inits <- list(
  list(alpha=0.01, beta=0.01, r=1,
       tau.phi=10, phi = rep(0.01,times=N)),  # chain 1
  list(alpha=0.5, beta=0.5, r=2,
       tau.phi=1, phi = rep(-0.01,times=N))   # chain 2
)
params <- c("sigma2.phi", "overallSIR", "resSIR", "SIR", "alpha", "phi","r")
ni <- 50000  # nb iterations
nt <- 10     # thinning interval
nb <- 10000   # nb iterations as burn-in
nc <- 2      # nb chains

ICARCodesamples <- nimbleMCMC(code = ICARCode,
                              data = ScotLipdata,
                              constants = ScotLipConsts,
                              inits = inits,
                              monitors = params,
                              niter = ni,
                              nburnin = nb,
                              thin = nt,
                              nchains = nc,
                              setSeed = 9,
                              progressBar = FALSE,    
                              samplesAsCodaMCMC = TRUE,
                              summary = TRUE,
                              WAIC = TRUE
)

Requesting help in this regard

Chris Paciorek

unread,
Dec 13, 2023, 10:43:44 AM12/13/23
to Shrinivas Shirke, nimble-users
Hi Shrinivas,

We have some suggestions on troubleshooting initialization problems like this in this section of the manual.

Please take a look and try those ideas and then feel free to follow up here if that doesn't help (and if so,
please indicate specifically what you found / where you got stuck).

-Chris

--
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/e61cd7a2-0442-4ad4-ba9f-19348cf552cfn%40googlegroups.com.

Shrinivas Shirke

unread,
Dec 17, 2023, 9:21:57 AM12/17/23
to nimble-users

Hello,
Thank you for the reply.
The code works if x is initialized . It works well when I am using nimbleMCMC.   
But when I am using 
ICAR <- nimbleModel(code = ICARCode, name = 'ICAR', constants = ScotLipConsts,
                    data = ScotLipdata, inits = inits)

it gives error message : Defining model
Error in if (initName %in% names(dL)) { : argument is of length zero
This error does not occur when I am using nimbleMCMC. Hence, I am not able to understand where I am making mistake. Also, I was trying nimbleHMC in place of nimbleMCMC but it states "Distribution dcar_normal does not have support for derivatives. This model cannot be compiled." 
Error: parameterTransform doesn't handle 'dcar_normal' distributions. 
Is there a way out in nimble for handling Intrinsic Autoregressive Models using HMC?

Chris Paciorek

unread,
Dec 21, 2023, 1:00:40 PM12/21/23
to Shrinivas Shirke, nimble-users
As far as the error, it looks like the problem is that the `inits` argument to `nimbleModel` needs a list containing a single set of initial values, while
you are passing a list of lists, with each component list containing a set of initial values.

We're going to fix this so that we handle this kind of input better - either giving a useful error message or just using the first list of the list of lists 
to initialize.

As far as CAR models with nimbleHMC, we don't support that at the moment. We'd like to make it possible, hopefully in the coming few months.

Chris Paciorek

unread,
Feb 4, 2024, 2:57:09 PM2/4/24
to Shrinivas Shirke, nimble-users
Hi Shrinivas,

An update on CAR models with nimbleHMC. With version 1.1.0, which we just released, it's now possible to use nimbleHMC sampling for model nodes following one of NIMBLE's CAR distributions.

-chris


Reply all
Reply to author
Forward
0 new messages