Hierarchical Model Using Normal Approximation to the Likelihood

Normal approximation to binomial likelihood is a classical method that is commonly used in meta-analysis. Let be the odds ratio for the th study:

You can estimate the treatment effect , , through the means of an approximate normal distribution,

where is the study-specific effect and is the known variance of . For , because all studies have large sample sizes (more than 50 persons in each group in nearly all of the studies), you can use the approximate sampling variance of :

If the odds ratios are exchangeable between the studies, then the treatment effect in each trial can be considered to be a random quantity drawn from some population distribution. In a Bayesian framework, this means that you can place a common prior distribution on the exchangeable random-effects parameters ,

where is the population average of the treatment effect across all studies and is the between-study variation. Suppose the following noninformative priors are placed on the hyperparameters and :

The following DATA step creates the response variable, , and its approximate sampling variance, for .

data multistudies;
   set multistudies;
   logy=log(trt/(trtN-trt)/(ctrl/(ctrlN-ctrl)));
   sigma2=1/trt+1/(trtN-trt)+1/ctrl+1/(ctrlN-ctrl);
run;

The following MCMC procedure statements fit the described random-effects model:

proc mcmc data=multistudies outpost=nlout seed=276 nmc=50000 thin=5
   monitor=(OR Pooled);
   array OR[22];
   parms  mu tau2;
   prior mu ~ normal(0, sd=3);
   prior tau2~ igamma(0.01,s=0.01);
   random theta ~n(mu, var=tau2) subject=studyid;
   OR[studyid]=exp(theta);
   Pooled=exp(mu);
   model logy ~ n(theta, var=sigma2);
run;

The PROC MCMC statement specifies the input data set as Multistudies and the output data set as Nlout, sets a seed for the random number generator, and requests a large simulation size with thinning. The MONITOR= option in the PROC MCMC statement outputs the analyses or saves the posterior samples for the specified variables. The ARRAY statement allocates an array variable OR of size 22, which is the number of random-effects parameters. The PARMS statement declares and as model parameters, and the PRIOR statements specify their prior distributions.

The RANDOM statement declares as random-effects parameters and specifies the normal prior distribution with hyperparameters and . The SUBJECT= option in the RANDOM statement specifies the group index as the data set variable studyid. The RANDOM statement generates 22 random-effects parameters, one for each unique value in the studyid variable. To monitor the function of 22 random-effects parameters, you must allocate an array and store the transformation of every random-effects parameter. The next programming statement, OR[studyid]=exp(theta) , does precisely this. As the procedure steps through the data set, it calculates the exponential of each parameter and stores the result in the th element of the OR array according to the studyid variable value (indexing). For this statement to work, the subject index variable studyid must have an integer value (1, 2, 3, and so on). The symbol Pooled calculates the exponential of the hyperparameter , which is the pooled treatment effect. The MODEL statement specifies the normal likelihood for the log odds ratios.

Figure 1 reports posterior summary statististics for each of the 22 odds ratios and for the Pooled estimate. The odds ratios provide the multiplicative change in the odds of mortality in the treatment group with respect to the control group. Some variation exists among the treatment effects. The pooled estimate of the odds ratio is 0.7849, indicating that mortality is less likely in the treatment groups.

You can use CATER autocall macro to create a caterpillar plot of the odds ratios. The following statement takes the output from the Nlout data set and generates a caterpillar plot of all study-specific odds ratios and the pooled odds ratio:

%CATER(data=nlout, var=OR: Pooled);

Figure 2 shows the posterior means and the 95 equal-tail credible intervals of odds ratios from the 22 studies and for the pooled odds ratio. It looks like the treatment is most effective for Study 7 and least effective for Study 14. This result illustrates another advantage of the random-effects models: a fixed-effects model focuses on the coefficient effects after averaging study differences, whereas a random-effects model can help identify studies that are very different from others, such as ones with very high or low treatment effects.

Figure 2 Caterpillar Plot of the Odds Ratios for the Normal Approximation Model

Hierarchical Model Using Binomial Likelihood

The standard normal approximation method might not be appropriate because the approximation becomes less precise in the extreme probabilities (for example, 0,1). In the Multistudies data set, the mortality rates in the thirteenth study and nineteenth study are closest to 0 (). Therefore, an alternative approach is to use the exact likelihood approach as opposed the normal approximation.

Emulating the binomial model as outlined in Spiegelhalter, Abrams, and Myles (2004), you can use the following model, where the treatment and control groups have their own binomial likelihood functions, with death probabilities and in the treatment group and control group, respectively:

Let the log odds for the control group be and let the log odds for the treatment group be as follows:

Then is the log odds ratio:

If treatment group risks are exchangeable, then you can assume that are random-effects parameters that are drawn from some common prior distribution:

If the control group risks are also exchangeable, then you can consider as other random-effects parameters with the following common prior:

Suppose the following priors are placed on the model parameters:

The following PROC MCMC statements fit the described model:

proc mcmc data=multistudies outpost=blout seed=126 nmc=50000 thin=5
   monitor=(OR Pooled);
   array OR[22];
   parms  mu_theta mu_phi s_theta s_phi;
   prior mu:~ normal(0,sd=3);
   prior s:~ igamma(0.01,s=0.01);
   random theta ~n(mu_theta, var=s_theta) subject=studyid;
   random phi ~n(mu_phi, var=s_phi) subject=studyid;
   p = logistic(theta + phi);
   q = logistic(phi);
   model trt ~ binomial(trtN, p);
   model ctrl ~ binomial(ctrlN, q);
   OR[studyid]=exp(theta);
   Pooled=exp(mu_theta);
run;

The PROC MCMC statement specifies the input data set as Multistudies and the output data set as Blout, sets a seed for the random number generator, and requests a large simulation size with thinning. The MONITOR= option in the PROC MCMC statement saves the posterior samples for the specified variables. The ARRAY statement allocates an array variable OR of size 22. The PARMS statement declares , and as model parameters, and the PRIOR statements specify their prior distributions. The RANDOM statements declare and as random-effects parameters and specify a normal prior distribution for them. The SUBJECT= option in both RANDOM statements indicates the group index as data set variable studyid. The LOGISTIC functions perform logistic transformation from the log odds of the treatment and control groups to the corresponding binomial mortality probabilities. The MODEL statements specify the binomial likelihood for the treatment and the control groups. The programming statement in the next line, OR[studyid]=exp(theta) , calculates the exponential of each ‘theta’ effect and stores the result in the OR array according to the studyid variable value (indexing). The symbol Pooled calculates the exponential of the hyperparameter , which is the pooled treatment effect.

Figure 3 reports posterior summary statististics for each of the 22 odds ratios.

You can use CATER autocall macro to create a caterpillar plot of the odds ratios. The following statement takes the output from the Nlout data set and generates a caterpillar plot of all study-specific treatment effects and the overall treatment effect:

%CATER(data=blout, var=OR: Pooled);

Figure 4 shows credible intervals of odds ratios from the 22 studies. Comparing Figure 2 and Figure 4, you can see that the posterior means and the credible intervals that are obtained by using the approximate normal likelihood and the binomial likelihood are very similar.

Figure 4 Caterpillar Plot of the Odds Ratios for the Binomial Model

Figure 5 compares the kernel density plots of the odds ratios that are produced by using the approximate normal likelihood and the exact binomial likelihood. The first plot shows the comparison for the thirteenth study, and the second plot shows the comparison for the nineteenth study. Even though the observed mortality rates () are close to 0 for both studies, the kernel density plots that are produced by using the approximate normal likelihood and the exact binomial likelihood are very similar. Hence, this example shows very little sensitivity to the likelihood assumption. In general, this can be explained by the large sample sizes of the studies included in the meta-analysis.

Figure 5 Kernel Density Comparison of Odds Ratios

Spiegelhalter, D. J., Abrams, K. R., and Myles, J. P. (2004), Bayesian Approaches to Clinical Trials and Health-Care Evaluation, Chichester, England: John Wiley & Sons.

Yusuf, S., Petro, R., Lewis, J., Collins, R., and Sleight, P. (1985), “Beta Blockade During and After Myocardial Infarction: An Overview of Randomized Trials,” Progress in Cardiovascular Diseases, 27, 335–371.