Fully Bayesian hierarchical modelling in two stages, with application to meta-analysis

Meta-analysis is often undertaken in two stages, with each study analysed separately in stage 1 and estimates combined across studies in stage 2. The study-specific estimates are assumed to arise from normal distributions with known variances equal to their corresponding estimates. In contrast, a one-stage analysis estimates all parameters simultaneously. A Bayesian one-stage approach offers additional advantages, such as the acknowledgement of uncertainty in all parameters and greater flexibility. However, there are situations when a two-stage strategy is compelling, e.g. when study-specific analyses are complex and/or time consuming. We present a novel method for fitting the full Bayesian model in two stages, hence benefiting from its advantages while retaining the convenience and flexibility of a two-stage approach. Using Markov chain Monte Carlo methods, posteriors for the parameters of interest are derived separately for each study. These are then used as proposal distributions in a computationally efficient second stage. We illustrate these ideas on a small binomial data set; we also analyse motivating data on the growth and rupture of abdominal aortic aneurysms. The two-stage Bayesian approach closely reproduces a one-stage analysis when it can be undertaken, but can also be easily carried out when a one-stage approach is difficult or impossible.


A.1 Generic code for stage two
The second stage of our method is implemented in OpenBUGS using the following code, for a p-dimensional set of parameters of interest.  (R[,], rho) Sigma[1:p, 1:p] <-inverse(Sigma.inv [,]) } The dummy[1:N] variable is merely a 'placeholder' for the missing likelihood from stage one. It simply allows us to create a link between each theta[i,] and the L samples generated for theta [i,] in stage one (samples[i, 1:p, 1:L]). Then we may specify level 2 onwards of the hierarchical model in the standard way. (Note that in BUGS normal distributions, both multivariate (dmnorm) and univariate (dnorm), are parameterised in terms of mean and precision (inverse-variance) rather than mean and variance.) The dproposem(.,.) syntax acts as a flag for BUGS to utilise a specialised updating algorithm, whereby the conditional prior for theta[i,] at its current value, conditional on mu and Sigma.inv, is compared to the conditional prior evaluated at a randomly chosen point in samples[i, 1:p, 1:L], in order to calculate the probability of moving theta[i,] to the randomly chosen point. Note that the values of dummy, samples, c, T, R and rho are specified in the data set. In cases where there is only one parameter of interest (p = 1), the following line is used instead, in combination with assumptions of univariate normality (dnorm), say, as opposed to multivariate normality (dmnorm).

dummy[i]~dpropose(theta[i], samples[i, 1:L])
The dummy[] variable must be 'observed' in order for the new mechanism to work, but its value is irrelevant. We typically specify dummy = c(1,1,...,1) in the data set. The mechanism is somewhat cumbersome but it works with much greater efficiency than more natural alternatives. For example, a more intuitive syntax might be etc. However, in this case the samples[i, 1:p, 1:L] matrix is considered as a multivariate stochastic quantity, which, internally in BUGS, necessitates each element being linked to every other. For large multivariate quantities this carries a substantial computational overhead, particularly when 'building' the data structures required for running the simulation. In the former approach, samples[i, 1:p, 1:L] is simply data and no such linking is required.

A.2 Transferring samples from stage one to stage two
Assuming the theta[1:N, 1:p] matrix has been monitored during stage one, we can access the required set of samples via the coda button on the Sample Monitor Tool. Two windows are created, one containing the samples themselves and another detailing which samples correspond to which elements of the monitored variable. The samples can be reformatted for incorporation into the second stage as follows. Edit the window containing the samples themselves so that a new line is created at the beginning with three numbers on it: (i) the number of studies; (ii) the number of parameters of interest; and (iii) the number of samples for each parameter. For example, we would enter '9 1 10000' for the pre-eclampsia example. Now select Reformat CODA from the Particle menu and a new window will appear containing the full sample for theta[1:N, 1:p] in the correct format for use with dproposem.

A.3 Checking flatness of implied priors
We can check for flatness of the prior distribution for any transformed parameter by simulating from the implied prior as follows. The BUGS code below is appropriate for any scalar θ = f (ϕ), where ϕ is a vector of parameters ϕ j ∼ N(0, 100 2 ), j = 1, ..., p. A sample from the prior p(θ) is obtained by simply storing the values of θ computed in each iteration; these can then be used to construct a density estimate. The final three lines of code impose the constraint a ≤ θ ≤ b, to restrict sampling to a specific region of interest, such as values supported by the likelihood.

A.4 Stage-one BUGS code for pre-eclampsia model
The following BUGS code and data list is used to perform stage one of our pre-eclampsia data analysis. Note that in BUGS the binomial distribution, dbin(), is parameterised such that the 'probability' and 'number of trials' parameters are the first and second arguments, respectively.  x.C = c(14,17,24,18,35,175,20,2,40), n.C = c (136,134,48,40,760,1336,524,103,102), x.T = c(14,21,14,6,12,138,15,6,65), n.T = c(131,385,57,38,1011,1370,506,108,153 Figure B1 shows Kaplan-Meier estimates (with 95% confidence intervals) of the survival function of the Cox-Snell residuals, defined as the estimated cumulative risk function evaluated at the observed event times. These are obtained from each study-specific model separately. If the assumed model fits the data well then we expect the Cox-Snell residuals to have a unit exponential distribution (solid grey line). For reliability of the plots, only studies with 5 or more ruptures are shown.  Figure B1: Kaplan-Meier estimates (black solid lines) with 95% confidence intervals (black dashed lines) of the survival function of Cox-Snell residuals. The solid grey line represents the unit exponential distribution. Only studies with 5 or more ruptures are shown.