Bayesian meta analysis of residual standard deviation using BUGS

Question

I have a small number of studies where the same model was calculated and want to infer the typical standard residual deviation. I have the degrees of freedom, the sum of the squares SSR and the mean of the square of the residuals MSR, but not the residuals themselves. My current plan is to do a bayesian meta analysis using either WinBugs or jags using two steps:

1) Assume the residual standard deviation is really the same over all studies and just try to estimate it more precisely. Let's call the overall residual variance res.var

With this, I already run into practical problems with the BUGS modelling language. I wanted to model the residual sum of squares using a chi square distribution, but have problems scaling the data accordingly. My problem is what I wanted to use is that for each study $$ \frac{SSR}{\text{res.var}} \sim \chi^{2}_{k}. $$ But I have problem putting this into a bugs model. I tried to use

model
{
for (i in 1:N)
{
chisquare[i] <- sumsq[i]/res.var
chisquare[i] ~ dchisqr(df[i])
}
res.var~dunif(0,1)
}

but this gives an multiple definitions of chisquare error. I tried transforming the first definition of chisquare so that sumsq is on the left side, but data can not be put on the left side of a logical node. Any practical advice?

2) Should I get this model to run, I would also like to allow for some variation of the residual standard deviation over studies, due to there being at least some small differences in the study population not included in the meta analysis model but I am not certain what a good hyperdistribution for the variance could look like. Any pointers?

Answer 1

Gauss it with ABC

Gauss it – start with the simplest version that has (most of) the elements of the real question.

ABC – Approximate Bayesian Computation – a direct two stage Monte-Carlo simulation where parameters are first drawn from the prior, then data is drawn using those parameters just drawn from the prior and only those draws that had exactly the same (or approximately) are kept. The distribution of parameters that are kept are a sample from the exact (approximate) posterior.

This can be quickly worked out for the linked k~ Bin(N,p) question (yours will require a generation of data that you only saw summaries of).

k~ Bin(N,p) question

Define a prior for N (preferably with no fixed upper bound but here probably OK)

Reps=10^6

N=sample(1:1000,Reps,replace=TRUE)

Draw possible k from Bin(N,k)

Pk=rbinom(Reps,N,prob=p) #recall p is known and k observed

Just keep draws where pk = k

Use that kept sample as an approximation of the posterior and say get a density estimate. (If you just want likelihood - no prior - then take c * posterior/prior as an approximation of that)

For example say the observed k was 28, this cab be easily run in R and checked against the closed form likelihood which is simple case of Example 4 Integrated Likelihood Methods for Eliminating Nuisance Parameters http://www.stat.duke.edu/~berger/papers/brunero.pdf

R code that can be run

Reps=10^7

N=sample(1:1000,Reps,replace=TRUE)

possible_k=rbinom(Reps,size=N,prob=.1)

posterior_k=N[possible_k==28]

hh=hist(posterior_k)

lik=dbinom(28,1:1000,prob=.1)

lines(1:1000,lik/max(lik) * max(hh$counts))

Now write the bugs or jags code, first for the tiny simple version of the question, note the posteriors are approx. the same, and then scale up to real question.

Maybe do this in private and then destroy any evidence of having to do this warm up exercise. (Recall that “people go Bats* crazy when they see others doing what they do for a living a different way” quote from the Moneyball movie.)

As for posterior probabilities, they are usually just formal and not of relevance to anyone. An exception would be the example by Jim Berger at 2009 International Workshop on Objective Bayes Methodology on HIV Vaccine trial. Really liked this from @Jason Aug 23’11

Bayesian methods are generative, in that they provide a complete "story" for how the data came into existence. Thus, they aren't simply pattern finders, but rather they are able to take into account the full reality of the situation at hand.

Answer 2

I solved the first technical issue on my own - the answer lies in the fact if X is Chi-Distributed with n degrees of freedom, then c*X is gamma distributed with $\alpha = n/2$ and $\beta = 1/2c$. So I just had to to write

sumsq[i] ~ dgamma(0.5*df[i],1/(2*res.var))

and then put a reasonable prior on it - the form of the prior did not matter much since due to the number of degrees of freedom the likelihood dominated the prior.

Answer 3

I have no idea what you are trying to do, but I can answer the technical question. I had exactly the same problem when redefining the node.

If sumsq[i]/res.var are known numbers, just compute the result in R and pass it as data into BUGS. That's all!

Regarding 2), I don't know if it helps in your case, but for uninformative prior for variance one usualy uses this:

tau <- 1/variance # normally, tau is used in dnorm etc.
tau ~ dgamma(0.01, 0.01)