Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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?

share|improve this question
    
While certainly not an expert at such things, I can try to help debug if you give some data. –  Matt Albrecht Nov 3 '12 at 1:16
    
Thanks, Matt! Sample data pasted above. –  user1172558 Nov 3 '12 at 2:24

2 Answers 2

up vote 2 down vote accepted

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).

require(rjags)
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)
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)
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)
id  <- rep(1:6, each=8)
df1 <- data.frame(x1, x2, y, id)
x   <- cbind(x1, x2)

# Data, Initial Values, and Parameters
N     <- length(unique(id))
W     <- dim(x)[2]
data  <- list(Ndata = length(y), N=N, M=6, W=W, y=y, x=x,id=id, m=c(-2,2), prec=c(.2,.2), tau.a=1, tau.b=2, tau.a2=1, tau.b2=2, mbeta0=1, precbeta0=.01)
inits <- list(beta0=0, beta=c(1,1), tau=1)
param <- c("beta0", "beta", "sigma", "delta")

# THE MODEL.
modelstring = "
model {
    for( r in 1 : Ndata ) {
        y[r] ~ dlnorm( mu[r] , tau ) #tau[subj[r]]
        mu[r] <- beta0 + delta[ id[r] ] + inprod(beta[] , x[r,]) #b1[ subj[r] ] * x1[r] + b2[ subj[r] ] * x2[r] - for subject varying slopes
    }
        beta0 ~ dnorm(mbeta0, precbeta0)
    for ( s in 1 : N ) {
        delta[s]  ~ dnorm(0, tau2)
#        b1[s]     ~ dnorm(m[l], prec[l]) #for subject varying slopes that do not correlate with the intercept
#        b2[s]     ~ dnorm(m[l], prec[l]) 
#        tau[s]    ~ dgamma( sG , rG ) #for subject varying tau in the likelihood
    }
    for (l in 1:W) {
        beta[l] ~ dnorm(m[l], prec[l])
    } 
    tau ~ dgamma(tau.a, tau.b)   
    sigma <- 1/sqrt(tau)
    tau2 ~ dgamma(tau.a2, tau.b2)
    sigma2 <- 1/sqrt(tau2)
}
" # close quote for modelstring
writeLines(modelstring,con="model.txt")

adaptSteps = 500              # Number of steps to "tune" the samplers.
burnInSteps = 500             # Number of steps to "burn-in" the samplers.
nChains = 3                   # Number of chains to run.
numSavedSteps=50000           # Total number of steps in chains to save.
thinSteps=1                   # Number of steps to "thin" (1=keep every step).
nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.

jagsModel = jags.model( "model.txt" , data=data , inits=inits , 
                    n.chains=nChains , n.adapt=adaptSteps )
cat( "Burning in the MCMC chain...\n" )
update( jagsModel , n.iter=burnInSteps )
cat( "Sampling final MCMC chain...\n" )
codaSamples = coda.samples( jagsModel , variable.names=param , 
                        n.iter=nPerChain , thin=thinSteps )
mcmcChain1 = as.matrix( codaSamples )
apply(mcmcChain1, 2, summary)

# check with lmer
require(lme4)
yl=log(y)
lmerx <- lmer(yl ~ x1 + x2 + (1|id))
lmerx

## User1172558 Model
modelstring = " 
    model {
    for(i in 1:N) {
        for (j in 1:M) {
            y[(i-1)*8+j] ~ dlnorm(mu[(i-1)*8+j], tau)
            mu[(i-1)*8+j] <- beta0 + delta[i] + inprod( beta[], x[(i-1)*8+j , ])
        }
        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)
}
" 
writeLines(modelstring, con="model.txt")
# Create, initialize, and adapt the model:
jagsModel = jags.model( "model.txt" , data=data , inits=inits , 
                    n.chains=nChains , n.adapt=adaptSteps )
update( jagsModel , n.iter=burnInSteps )
codaSamples = coda.samples( jagsModel , variable.names=param , 
                        n.iter=nPerChain , thin=thinSteps )
mcmcChain2 = as.matrix( codaSamples )
apply(mcmcChain2, 2, summary)
lmerx

require(MCMCglmm) #Check with MCMCglmm
mc1 <- MCMCglmm(yl ~ x1 + x2, random=~id, data=df1)
summary(mc1); lmerx; apply(mcmcChain1, 2, summary); apply(mcmcChain2, 2, summary)

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.

share|improve this answer
    
Thanks, Matt! But this was not what I was looking for. The nested loop was necessary given I have longitudinal data. It means, the first 8 rows are for municipality A, the next 8 rows for municipality B, and so on. The way you have it now does not allow for this structure of the data. I also need to use lognormal, not a normal. My y is very skewed to the right. –  user1172558 Nov 3 '12 at 17:29
    
Yes, I forgot about the dlnorm when I was experimenting. I've updated the answer to include this. For the general model, I understood it to indicate different intercepts per municipality but common b1 and b2 slopes across municipality. This was because the 'i' indicated municipality level intercepts. While in the model the 'k' was a proxy for the entirety of the data. If what you want is different slopes per municipality I will update it (they're the sections that have been commented out). –  Matt Albrecht Nov 3 '12 at 18:12
    
Also, while I've taken out the nested loop, I'm pretty sure it's specifying the same model: delta[i] was fitting different intercepts per group and this is what delta[id[r]] does. In lmer format what I've fitted is lmer(y ~ x1 + x2 + (1|municipality)) but with a lognormal distribution. Does this make sense? –  Matt Albrecht Nov 3 '12 at 18:19
    
Matt, can you paste the r code you are using to run your JAGS code? I haven't been able to run your code yet. Thanks –  user1172558 Nov 3 '12 at 19:20
    
@user1172558 Done. Plus extra comments at the end. I'm not sure why we have a relatively large divergence in the posteriors for our different methods. –  Matt Albrecht Nov 4 '12 at 3:37

I found the problem.

JAGS does not accept to have an object receiving multiple values inside different loops. I changed the JAGS code and it worked. See Below:

JAGS CODE:

model {
    for(i in 1:N) {
        for (j in 1:M) {
            y[(i-1)*8+j] ~ dlnorm(mu[(i-1)*8+j], tau)
            mu[(i-1)*8+j] <- beta0 + delta[i] + inprod(x[(i-1)*8+j,], 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)
}

RJAGS CODE:

# Set up Data
y <- br[, 28]
x <- as.matrix(br[, c(29,30,12,22)])

data <- list(N=length(unique(br$munname)), M=8, W=dim(x)[2], 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(delta=rep(0,length(unique(br$munname))), beta0=0, beta=rep(1,dim(x)[2]), tau=1, tau2=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)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.