Hi,
My script for nimbleMCMC is like this, and in each MC loop, I will regenerate the x, y, z, and area information
`code <- nimbleCode({
beta0 ~ dnorm(0, sd = 5)
beta1 ~ dnorm(0, sd = 5)
gamma0 ~ dnorm(0, sd = 5)
gamma1 ~ dnorm(0, sd = 5)
phi ~ dgamma(2, 0.05)
alphaDP ~ dgamma(1, 10)
for (h in 1:(K - 1)) {
v[h] ~ dbeta(1, alphaDP)
}
v[K] <- 1
stick_left[1] <- 1
for (h in 1:(K - 1)) {
omega[h] <- v[h] * stick_left[h]
stick_left[h + 1] <- stick_left[h] * (1 - v[h])
}
omega[K] <- stick_left[K]
for (h in 1:K) {
sigma_b[h] ~ dunif(0.001, 2)
sigma_m[h] ~ dunif(0.001, 2)
rho[h] ~ dunif(-0.95, 0.95)
cond_coef[h] <- rho[h] * sigma_m[h] / sigma_b[h]
cond_sd_m[h] <- sigma_m[h] * sqrt(1 - rho[h] * rho[h])
}
for (a in 1:A) {
cl[a] ~ dcat(omega[1:K])
b_raw[a] ~ dnorm(0, sd = sigma_b[cl[a]])
mean_m_cond[a] <- cond_coef[cl[a]] * b_raw[a]
m_raw[a] ~ dnorm(mean_m_cond[a], sd = cond_sd_m[cl[a]])
b[a] <- b_raw[a]
m[a] <- m_raw[a]
}
for (j in 1:N) {
logit(mu[j]) <- beta0 + beta1 * x[j] + b[area[j]]
logit(p_one[j]) <- gamma0 + gamma1 * z[j] + m[area[j]]
y[j] ~ dOneInflBeta(p_one[j], mu[j], phi)
}
})
constants <- list(
N = N,
A = A,
K = K,
area = area_samp,
x = x_samp,
z = z_samp
)
data <- list(
y = y_samp
)
inits_function <- function() {
list(
beta0 = 0,
beta1 = 0,
gamma0 = 0,
gamma1 = 0,
phi = 40,
alphaDP = 0.5,
v = c(rbeta(K - 1, 1, 2), 1),
sigma_b = runif(K, 0.05, 0.5),
sigma_m = runif(K, 0.05, 0.5),
rho = runif(K, -0.3, 0.3),
cl = sample(1:K, A, replace = TRUE),
b_raw = rnorm(A, 0, 0.1),
m_raw = rnorm(A, 0, 0.1)
)
}
monitors <- c(
"beta0", "beta1",
"gamma0", "gamma1",
"phi",
"alphaDP",
"omega",
"sigma_b", "sigma_m", "rho",
"b", "m",
"cl"
)
fit <- nimbleMCMC(
code = code,
constants = constants,
data = data,
inits = inits_function,
monitors = monitors,
niter = 6000,
nburnin = 2000,
thin = 2,
nchains = 3,
samplesAsCodaMCMC = TRUE)`