With this data:
y <- c(1.105808, 1.000000, 5.304166, 33.665875, 139.865451, 109.033703, 176.639245, 1.000000, 28.521134, 44.281121 ,150.478570, 18.465554, 85.096431, 81.907537, 124.631226, 1.000000 , 11.237294 , 20.480519 , 68.642176 , 30.047630, 54.051613 , 134.068889 , 72.215041 , 1.000000 , 31.254480, 6.226026 , 54.340496 , 161.667352 ,1345.948800 , 147.404744, 192.966923 , 1.000000 , 1.150755 , 1.000000 , 99.477430, 72.592598, 107.493014 , 130.201416 , 147.387423 , 1.000000, 27.534944 , 8.492657 , 38.155558, 44.301978 , 81.938633, 75.026848 ,144.926523 , 1.000000 , 1.075801 , 40.001560)
x1 <- c(7.878865, 9.117159 , 9.998539 , 10.300563 , 12.683197 , 12.185060, 101.346385, 168.814861, 4.769977, 4.803769, 5.990192, 6.469412, 7.557664, 9.781595, 264.102447 ,702.321019 , 5.663501 , 6.319843, 6.643405 , 7.147517 , 8.154099 , 8.811370 , 64.089236, 205.728163, 7.218225 , 6.905615 , 9.341990 , 8.554343 , 15.873037 , 15.731554, 227.589294 ,398.435765 , 6.217681 , 10.498929 , 11.088663 , 10.312797, 11.483123 ,9.276521 ,157.311069 ,391.279665 , 3.985544 , 4.389385, 4.663445 , 4.934453 , 6.622184 , 7.770833 , 21.700911 , 33.470576, 5.405404 , 6.531107)
x2 <- c(0.4225000, 0.6619411 ,0.5401000, 0.5138000, 0.5109000, 0.6325681, 0.6425919, 0.6466943, 0.4421000, 0.5870430, 0.5254000, 0.5525000, 0.5392000, 0.6330954, 0.6457942, 0.6582531, 0.4039000, 0.6154006, 0.4620000, 0.4875000, 0.5439000, 0.6494275, 0.6423681, 0.6450814, 0.5110000, 0.6060935, 0.5050000, 0.5200000, 0.5535000, 0.6383363, 0.6684222 ,0.6465682, 0.3963000, 0.5914320, 0.4819000, 0.5361000, 0.5886000, 0.6140674 ,0.6240171, 0.6150484, 0.4797000, 0.6211242, 0.5705000, 0.5709000, 0.6144000 ,0.6412593, 0.6611542, 0.6364444, 0.3599000, 0.6375195)
I'm trying to fit a lognormal random effects model in JAGS. Below my JAGS code:
# Lognormal Model
# N municipalities
# M years
# W Betas
model {
for(i in 1:N) {
for (j in 1:M) {
k <- (i-1)*8 + j
y[k] ~ dlnorm(mu[k], tau)
mu[k] <- beta0 + delta[i] + inprod(x[k,], beta[])
delta[i] ~ dnorm(0, tau2)
}
}
# Prior for betas
beta0 ~ dnorm(mbeta0, precbeta0)
for (l in 1:W) {
beta[l] ~ dnorm(m[l], prec[l])
}
# Prior for precision of Y
tau ~ dgamma(tau.a, tau.b)
sigma <- 1/sqrt(tau)
# Prior for precision of Delta
tau2 ~ dgamma(tau.a2, tau.b2)
sigma2 <- 1/sqrt(tau2)
}
Here my Rjags code:
br <- read.csv("Data/test.csv", header=T, sep=',')
br <- na.omit(br)
br$y <- br$y + 1
# Set up Data
y <- br[, 28]
x <- br[, c(29,30,12,22)]
# Data, Initial Values, and Parameters
N <- length(unique(br$ID))
W <- dim(x)[2]
data <- list(N=N, M=8, W=W, y=y, x=x, m=c(-2,2,4,3), prec=c(.20,.20,.20,.20), tau.a=1, tau.b=2, tau.a2=1, tau.b2=2, mbeta0=1, precbeta0=.01)
inits <- rep(list(list(beta0=0, beta=c(1,1,1,1), tau=1)),5)
param <- c("beta0", "beta", "sigma", "delta")
# Model
sim <- jags.model(file="Code/lognormal2.bug", data=data, inits=inits, n.chains=5, n.adapt=1000)
But I get the following ERROR:
Error in jags.model(file = "Code/lognormal2.bug", data = data, inits = inits, :
RUNTIME ERROR:
Compilation error on line 11.
Missing values in subset expression of y
Any help?
Best Answer
I took out the nested loop to make it clearer for me and stuck the prior for delta in its own loop, giving this for mu: delta[i] -> delta[id[r]]
note the explicit id vector for municipalities that is as long as the data which I've labelled 'id'
I think that's the biggest change I made? Model below (also some extras in case you want to include them).
Edit: Included full script.
My version seems to be more consistent with lmer and MCMCglmm. So I'm not so sure what the difference is between your nested loop for the municipality level intercept and my method of doing it. Did you model your code from somewhere? I would also consider weakening the priors, unless you have a strong reason for them.