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