# JAGS meta-analysis

I'm using JAGS to run a meta-analysis and I've run into an issue. I have calculated log-response ratios and errors for about 177 studies, where studies fall into one of 8 groups. I'm interested in estimating an overall effect as well as a effects for each group. I've discovered that if I use the model

model{
for(i in 1:N){
RR[i] ~ dnorm(yhat[i], prec_y[i])
yhat[i] <- mu + gamma[type[i]]
}

for(j in 1:J){
gamma[j] ~ dnorm(gamma_mu, gamma_prec)
gamma_true[j] <- gamma[j] - mean(gamma[])
}

mu_true <- mu + mean(gamma[])

mu ~ dnorm(0, 0.001)
mu_prec ~ dgamma(0.01, 0.01)
gamma_mu ~ dnorm(0, 0.001)
gamma_prec ~ dgamma(0.01, 0.01)
}


than the model doesn't fit my data well, because the estimate for Plant-Mych is just way off:

sample image http://www.natelemoine.com/not_working.png

Alternatively, if I use the model:

model{
for(i in 1:N){
RR[i] ~ dnorm(yhat[i], prec_y[i])
yhat[i] ~ dnorm(gamma[type[i]], gamma_prec[type[i]])
}

for(j in 1:J){
gamma[j] ~ dnorm(overall_mu, mu_prec)
gamma_prec[j] <- pow(gamma_sigma[j], -2)
gamma_sigma[j] ~ dunif(0, 100)
}

overall_mu ~ dnorm(0, 0.001)
mu_prec ~ dgamma(0.01, 0.01)

}


then it fits much better:

sample image http://www.natelemoine.com/working.png

Is the fundamental difference that, in the first bad model, the groups are fixed effects whereas in the second model the groups are random? If so, why does this make such a big difference? Is there a way to code the first model so that it fits better? I appreciate the help!

In order to clarify the difference I'm going to simplify and rearrange

The first model becomes

model{

mu ~ dnorm(0, 0.001)
gamma_mu ~ dnorm(0, 0.001)

gamma_prec ~ dgamma(0.01, 0.01)
for(j in 1:J){
gamma[j] ~ dnorm(gamma_mu, gamma_prec)
}

for(i in 1:N){
yhat[i] <- mu + gamma[type[i]]
RR[i] ~ dnorm(yhat[i], prec_y[i])
}
}


and the second model becomes

model{

overall_mu ~ dnorm(0, 0.001)
mu_prec ~ dgamma(0.01, 0.01)

for(j in 1:J){
gamma_sigma[j] ~ dunif(0, 100)
gamma[j] ~ dnorm(overall_mu, mu_prec)
}

for(i in 1:N){
yhat[i] ~ dnorm(gamma[type[i]], gamma_sigma[type[i]]^-2)
RR[i] ~ dnorm(yhat[i], prec_y[i])
}
}


To answer your question both models have the effect of group J on the mean level modelled as a random effect. The first model suffers from multicollinearity because it specifies an intercept mu but then models the random effects as varying about a second non-independent intercept gamma_mu. The second model seems to suffer from a similar problem in the sense that there are two variation terms mu_prec and gamma_sigma being used to model the within group variation. Furthermore both models are incomplete because there is no prior for prec_y

What I think you really want is something like one of the the following.

The first has the differences between the group means as a random effect and the standard deviation for each group as a fixed effect.

model{

mu ~ dnorm(0, 5^-2)

mu_sd ~ dunif(0, 5)
for(j in 1:J){
gamma_sigma[j] ~ dunif(0, 5)
gamma[j] ~ dnorm(0, mu_sd^-2)
}

for(i in 1:N){
yhat[i] ~ mu + gamma[type[i]]
RR[i] ~ dnorm(yhat[i], gamma_sigma[i]^-2)
}
}


While the second has the differences between the group means as a random effect and the standard deviation for each group as a random effect.

model{

mu ~ dnorm(0, 5^-2)

mu_sd ~ dunif(0, 5)
mu_gamma_sigma ~ dnorm(0, 5^-2)
sigma_sigma ~ dunif(0, 5)
for(j in 1:J){
gamma_sigma[j] ~ dlnorm(mu_gamma_sigma, sigma_sigma^-2)
gamma[j] ~ dnorm(0, mu_sd^-2)
}

for(i in 1:N){
yhat[i] ~ mu + gamma[type[i]]
RR[i] ~ dnorm(yhat[i], gamma_sigma[i]^-2)
}
}


I'd appreciate it if you would rerun both new models and confirm they don't have any errors and post the results below.

• I just got it that you are passing prec_y as data so my criticism of the second model is invalid and my suggestions too simple. As you identify your second model in the question is correct. However your suggested answer suffers from the multicollinearity I identified - to eliminate replace gamma_mu with 0 so that each of the gammas is drawn from a normal distribution with a mean of zero, i.e., gamma[j] ~ dnorm(0, gamma_prec). As a result mu will be the overall mean and mu + gamma[j] will be the mean for each group. – joethorley Apr 20 '14 at 2:55

First, thanks joethorley for the answer. But I figured it out and the answer is fairly simple. The first incorrect model can be made equivalent to the second, correct model by:

model{
for(i in 1:N){
RR[i] ~ dnorm(yhat[i], prec_y[i])
yhat[i] ~ dnorm(theta[i], sd_y[type[i]])
theta[i] <- mu + gamma[type[i]]
}

for(j in 1:J){
gamma[j] ~ dnorm(gamma_mu, gamma_prec)
gamma_true[j] <- gamma[j] - mean(gamma[])
sd_y[j] ~ dunif(0, 10)
}

mu_true <- mu + mean(gamma[])

mu ~ dnorm(0, 0.001)
mu_prec ~ dgamma(0.01, 0.01)
gamma_mu ~ dnorm(0, 0.001)
gamma_prec ~ dgamma(0.01, 0.01)
}


In the first, incorrect model, the mean of each study was assumed to vary randomly around the 'true' study mean. Then, the true study mean was a fixed, deterministic outcome of the overall mean plus the effect of type.

In the second model, both the original formula and this one, the reported mean of each study still varies randomly around the true mean. Now, the 'true' means of each study are replicates that vary randomly around the mean of each type (which is either estimated on its own, as in the first correct model or a function of the overall mean + an effect of study type, as here). The problem is that I wasn't allowing for the 'true' response to vary around the mean of each type, which I should have done (because the studies are replicates and therefore expected to vary around some true mean of the type).

The results look great: http://www.natelemoine.com/really_working.png