ORIGINAL ARTICLE
Use of Generalized Linear Mixed Models for Network Meta-analysis Yu-Kang Tu, PhD
In the past decade, a new statistical method—network meta-analysis—has been developed to address limitations in traditional pairwise meta-analysis. Network metaanalysis incorporates all available evidence into a general statistical framework for comparisons of multiple treatments. Bayesian network meta-analysis, as proposed by Lu and Ades, also known as ‘‘mixed treatments comparisons,’’ provides a flexible modeling framework to take into account complexity in the data structure. This article shows how to implement the Lu and Ades model in the frequentist generalized linear mixed model. Two examples are provided to demonstrate how centering the covariates for random effects estimation within each trial can yield
correct estimation of random effects. Moreover, under the correct specification for random effects estimation, the dummy coding and contrast basic parameter coding schemes will yield the same results. It is straightforward to incorporate covariates, such as moderators and confounders, into the generalized linear mixed model to conduct meta-regression for multiple treatment comparisons. Moreover, this approach may be extended easily to other types of outcome variables, such as continuous, counts, and multinomial. Key words: randomized controlled trials; network meta-analysis; mixed treatments comparisons; generalized linear mixed models. (Med Decis Making 2014;34:911–918)
M
comparison. The evidence base, therefore, varies across comparisons. One potential consequence is that the results from multiple pairwise meta-analyses are not always consistent, and some comparisons may not have been conducted yet.11–13 Network meta-analysis incorporates all available evidence into a general statistical framework for comparisons of all available treatments to resolve the limitations of traditional pairwise meta-analysis. The idea of indirect comparisons was first proposed in the 1990s.14,15 Lumley2 coined the term ‘‘network meta-analysis’’ in 2002 and proposed a linear mixed-model approach to comparisons of multiple treatments within the same statistical model. Lu and Ades6 developed a sophisticated Bayesian hierarchical model, providing a flexible statistical framework to take into account the complexity of the data structure within multiarm trials. Their statistical approach, widely known as mixed treatments comparison or Bayesian network meta-analysis, is currently the standard approach to comparisons of multiple treatments.7,16–18 Comparing all available treatments simultaneously is an attractive idea to clinicians and health care providers. The number of articles that report the results of network meta-analysis has increased exponentially in the past few years, as have their citations,18,19 despite the fact that implementing
eta-analysis is an important research tool for synthesis of evidence on the efficacy versus harm of available medical interventions.1 One recent development in meta-analysis methodology is network meta-analysis for comparisons of multiple treatments.2–10 When several treatments are available for a medical condition, traditional metaanalysis undertakes multiple pairwise comparisons, using different sets of clinical trials for each
Received 24 April 2014 from the Institute of Epidemiology & Preventive Medicine, College of Public Health, National Taiwan University, Taipei, Taiwan (Y-KT). The author declares no conflict of interests. This project was partly funded by a grant from the National Science Council in Taiwan (Grant No. NSC 101-2314-B-002-197-MY2). Revision accepted for publication 30 June 2014. Supplementary material for this article is available on the Medical Decision Making Web site at http://mdm.sagepub.com/supplemental. Address correspondence to Yu-Kang Tu, Institute of Epidemiology & Preventive Medicine, College of Public Health, National Taiwan University, 17 Xu-Zhou Road, Taipei, Taiwan; telephone: (886) 2 3366 8039; fax: (886) 2351 1955; e-mail: [email protected]. Ó The Author(s) 2014 Reprints and permission: http://www.sagepub.com/journalsPermissions.nav DOI: 10.1177/0272989X14545789
MEDICAL DECISION MAKING/OCTOBER 2014
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
911
TU
Bayesian network meta-analysis is quite challenging. Several recent articles have attempted to elaborate further on the statistical issues related to implementing the Lu and Ades approach within the generalized linear mixed-modeling framework using nonBayesian statistical software packages,3,7,17,20 with the intention of making network meta-analysis more accessible to clinicians and meta-analysts who are not yet familiar with Bayesian statistics. This article provides a general, non-Bayesian approach to network meta-analysis based on the Lu and Ades model. I first demonstrate how a pairwise meta-analysis may be implemented in a generalized linear mixed model and show how different parameterizations may affect results from random-effects models. I then extend the discussion to multiple treatment comparisons and demonstrate how the Lu and Ades model may be implemented in a generalized linear mixed model with an appropriate parameterization of random effects estimation. Finally, I provide an example to show that 2 proposed generalized linear mixed-modeling approaches to network metaanalysis are different parameterizations of the same model. GENERALIZED LINEAR MIXED MODEL FOR PAIRWISE META-ANALYSIS Fixed-effects Model for Pairwise Meta-analysis A fixed-effects model assumes a common effect behind the observed effects, while a random-effects model assumes that the ‘‘true’’ effect follows a distribution. The basic fixed-effects model for comparing multiple treatments is as follows21:
g yij 5 sj 1 ti ;
ðEquation 1Þ
where yij is the outcome for treatment arm i in study j; sj is the study effect; and ti is the treatment effect. The function g(.) is a suitable link function for transformation of yij. For example, when the outcome is probability of an event for a subject receiving treatment i in study j, the link function is logit. In a traditional pairwise meta-analysis with p studies that compare treatment B to treatment A, the generalized linear regression model for Equation 1 may be expressed as follows22: g y^ij 5 b0 1 bj studyj 1 bp Ti ; i ¼ A; B; j ¼ 1; :::; p 1
ðModel 1Þ
912
where y^ij is the estimated outcome for each arm i in study j, b0 is the intercept, and b1 to bp–1 are regression coefficients for study1 to studyp–1, respectively, while bp is the regression coefficient for dummy variable Ti, which is coded as 1 for treatment i = B and 0 for treatment i = A. Consequently, b0 is the estimated treatment effect for A in study p; and bj, where j = 1 to p – 1, is the difference in treatment effects between either treatment i in study j and the same treatment i in study p. The regression coefficient, bp, is the average difference in treatment effects between B and A across all trials. The Lu and Ades model for fixed-effects metaanalysis comparing 2 treatments may be written as follows4,21: hkj 5
mj k 5 A ; mj 1 d k 5 B
ðEquation 2Þ
where hkj is the estimated outcome for treatment k in study j—that is, gð^ ykj Þ in Model 1. Consequently, d is the difference in treatment effects between A and B. The subscript j for m indicates that the treatment effect for A is estimated separately for each study, so Equation 2 and Equation 1 are equivalent and may be implemented in Model 1. Random-effects Model for Pairwise Meta-analysis Equation 2 assumes that the difference in treatment effects between A and B is a fixed value, d, which is estimated by bp in Model 1. In a randomeffects models, the difference in treatment effects is no longer a fixed value but rather a random variable that follows a normal distribution4,21:
k 5A mj mj 1 dj k 5 B dj ;Nðd; s2 Þ: hkj 5
ðEquation 3Þ
The mean of dj is d and its variance is s2 (i.e., the difference in treatment effects between A and B is a random variable). The generalized linear mixed regression model for Equation 3 may be written as g y^ij 5 b0 1 bj studyj 1 bpj Ti ; i ¼ A; B; j ¼ 1; :::; p 1 bpj 5 bp 1 uj ;
ðModel 2Þ
where y^ij is the estimated treatment effect for each arm i in study j, and b0 is the intercept, b1 to bp–1 are regression coefficients for study1 to studyp–1, respectively, and bpj is the regression coefficient for
MEDICAL DECISION MAKING/OCTOBER 2014
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
GLMM FOR NETWORK META-ANALYSIS
dummy variable Ti in study j. The average difference in treatment effects is bp , and uj is the random effect for bp (i.e., the difference between bp and bpj Þ. However, Model 2 needs some fine-tuning before it will yield a correct estimation for bp and the variance of uj, as the following example illustrates. Colditz and others23 provided data from 13 randomized controlled trials on the efficacy of Bacillus Calmette–Gue´rin (BCG) vaccine in the prevention of tuberculosis. Applying data presented by Borenstein and others24 and using the command metan for the statistical software Stata (College Station, Texas) will allow us to undertake DerSimonian and Laird random-effects meta-analysis.25 The pooled odds ratio (OR) is 0.474 (95% confidence interval [CI] 0.325–0.691). The between-study variance, usually denoted as t2 , is 0.3682. To undertake the random-effects analysis using a generalized linear mixed model, such as Model 2, one may use the Stata command melogit or meqrlogit in version 13 or xtmelogit in version 12. The following Stata code used melogit: melogit r bcg study1-study12, || trial: bcg, noconstant binomial(n) or
Each study has 2 data entries: 1 data entry for the BCG-vaccinated group and 1 data entry for the control group. Variable r is the number of tuberculosis cases; n is the total number of subjects in each group. Variable trial was the identifier for each trial, and study1 to study12 are 12 dummy variables for trial 1 to trial 12, respectively. The option noconstant indicates that the intercept is excluded from the estimation of random effects (i.e., this model only estimated random effects for the slope of variable bcg, which is a dummy variable with the vaccinated group coded as 1 and the control group coded as 0). The difference in treatment effects among the control groups from different trials is estimated as fixed effects by the dummy variables study1 to study12. The variations in the treatment effects of control groups may be estimated as a random variable, but, as Senn22 noted, this is not the standard approach in meta-analysis. The regression coefficient for bcg is –0.764 (OR 0.466; 95% CI 0.332–0.654), which is smaller than that given by metan; s2 , the estimated variance for random effects, is 0.2811. To obtain a similar OR, one must center bcg (i.e., recode the vaccinated group as 0.5; control group, –0.5). The results then show that the OR for centered bcg is 0.475 (95% CI 0.337–0.770) and s2 is 0.295. Note that for the fixed-effects model, centering would not cause a difference in the OR.
Figure 1 illustrates how centering causes the differences in random effects estimation. When the control group is coded 0 and the BCG-vaccinated group is coded 1, random effects, as a random slope for bcg, are applied to the vaccinated group only (Figure 1A); when bcg is centered, random effects are applied to both groups (i.e., to the difference between vaccinated and control groups, Figure 1B).7 Consequently, the random slope for centered bcg can now rotate around the midpoint between the vaccinated and control groups, so the difference between the 2 groups may be estimated as random effects under the constraint that the sum of the random effects within a trial is zero.7 The centered bcg therefore gives the model greater ‘‘freedom’’ to estimate the differences between the 2 groups and should yield a better model fit: The log likelihoods for uncentered and centered models are –90.50 and –88.64, respectively. This approach with centered bcg is similar to the Bayesian meta-analysis model proposed by Smith and others26: di 5 logit pTi logit pCi T mi 5 ðlogit pi 1 logit pCi Þ=2 ; logit pCi 5 mi 12 di T 1 logit pi 5 mi 1 di 2
where pTi and pC i are the event probability for treatment and control groups. Lu and Ades6 extended this model to multiple treatment comparisons. Centering the covariates for random effects estimation is crucial, especially when the heterogeneity across trials is large. For traditional meta-analysis, the 2 groups are coded – 0.5 and 0.5 after centering. For network meta-analysis with multiarm trials, the mean value of a dummy variable is subtracted from its individual value. For instance, in a 3-arm trial, treatments A, B, and C would be coded 1, 0, and 0, respectively, in the dummy variable for treatment A. However, after centering, they would instead be coded 2/3, –1/3, and –1/3, respectively. GENERALIZED LINEAR MIXED MODEL FOR NETWORK META-ANALYSIS Fixed-effects Model for Network Meta-analysis We must now extend our generalized linear mixed model to multiple treatment comparisons. For the fixed-effects model, the basic equations are the same as those for Model 1, with more treatments to compare. The regression model may be written as
ORIGINAL ARTICLE
913
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
TU
Figure 1 Diagram for the impact of centering on random effects estimation. (A) Dummy coding for variable bcg where the control group is coded 0 and the vaccinated group coded 1. (B) Centered bcg where the control group is coded –0.5 and the vaccinated group coded 0.5.
g y^ij 5 b0 1 bj studyj þ ci Ti ; i ¼ 1; :::; k; j ¼ 1; :::; p: ðModel 3Þ
where y^ij is the estimated outcome for each arm i in study j, b0 is the intercept, and b1 to bp are regression coefficients for dummy variables study1 to studyp, respectively. Furthermore, c1 to cK are the regression coefficients for dummy variables T1 to TK, where treatment k is coded 1 and the other arms are coded 0 for Tk. Suppose one uses the first treatment and the first trial as the primary reference groups, so b1 and c1 are fixed at 0. Consequently, b0 is the estimated treatment effect for treatment 1 in study 1; and bj, where j = 2 to p, is the difference in treatment effects between treatment i in study j and treatment i in study 1. The regression coefficient ci is the difference in treatment effects between treatment i and treatment 1 across all trials. The fixed-effects Lu and Ades model may be specified as 8 mjb > > < b 5 A; B; C; . . . if k 5 b hjk 5 m 1 d 5 m 1 d d ; > bk Ak Ab jb > : jb 0 k 5 B; C; D; . . . if k is 0 after B
ðEquation 4Þ
where hjk is the estimated outcome for arm k in study j, and mjb is the baseline treatment in trial j. The difference between the other treatment k and treatment b in the same trial will be estimated by expressing them in
914
terms of effects relative to treatment A, which is the primary reference treatment within the network. For reasons of identification and because of its interpretation as the effect of treatment A compared with itself, dAA is fixed at 0, and Lu and Ades called dAB to dAk the ‘‘basic parameters.’’ The advantage of expressing all treatment comparisons as the relation between basic parameters is that the number of pairwise comparisons to be estimated for a network meta-analysis involving k treatments is reduced to k 2 1 for the fixed effects.27 For instance, suppose one has 3 trials. Trial 1 compares treatments A, B and C, while trial 2 compares A and B, and trial 3 compares B and C. Table 1 shows the difference in coding the covariates for treatments. Variables tA, tB, and tC are dummy variables for the 3 treatments, and because A is the primary reference treatment, only tB and tC are included in the model. Variables tAA tAB, and tAC use the contrast coding for basic parameters.27 Furthermore, because the regression coefficient for tAA is fixed at 0, only tAB and tAC are included in the model. Both coding schemes will yield the same results for the fixedeffects model. If we exclude the intercept and include all 3 dummy variables for the trials (i.e., trial 1, trial 2, and trial 3 in Table 1) in the model, the basic parameter coding has 1 advantage, which is that the dummy variables for each trial have a clear interpretation: They represent the trial-specific baseline treatments for each trial.
MEDICAL DECISION MAKING/OCTOBER 2014
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
GLMM FOR NETWORK META-ANALYSIS
Table 1
Two Coding Schemes for Covariates Included in the Generalized Linear Mixed Model for Treatment Effects Estimation Dummy Coding
Trial
1 1 1 2 2 3 3
Contrast Coding
Study 1
Study 2
Study 3
Treatment
Baseline
tA
tB
tC
tAA
tAB
tAC
1 1 1 0 0 0 0
0 0 0 1 1 0 0
0 0 0 0 0 1 1
A B C A B B C
A A A A A B B
1 0 0 1 0 0 0
0 1 0 0 1 1 0
0 0 1 0 0 0 1
0 –1 –1 0 –1 0 0
0 1 0 0 1 0 –1
0 0 1 0 0 0 1
Sample Data The sample dataset contains the results of 26 studies that directly compared 3 treatments for the prevention of first bleeding in patients with liver cirrhosis28: Treatment A was the control, nonactive treatment; B, sclerotherapy; and C, beta-blocker. (The dataset is available in Table 1 of the online appendix.) Among the 26 studies, 2 studies were 3arm trials, 7 studies compared A with C, and 17 studies compared B with C. The following Stata code undertakes fixed-effects network meta-analysis using the contrast coding scheme for treatments: melogit r tAB tAC study1-study26, noconstant binomial(n) /*code-1*/
The option noconstant removes the intercept from the model, so all the dummy variables, study1 to study26, for estimating the effects of the baseline treatment for the 26 trials can be included. The variable r is the number of events (i.e., bleeding); n is the number of patients in each arm. Variables tAB and tAC use the contrast coding scheme for the differences between B and A and between C and A, respectively, such as tAB and tAC in Table 1. The same results may also be obtained using the following dummy coding scheme: melogit r tB tC study1-study26, noconstant binomial(n) /*code-2*/
Variables tB and tC are dummy variables for treatments B and C, shown as tB and tC in Table 1. Random-effects Model for Network Meta-analysis For the random-effects network meta-analysis, dbk in Equation 4 is replaced by djbk , the trial-specific effect of treatment k relative to the trial-specific baseline treatment b. These trial-specific effects are then
drawn from a normal distribution: djbk ;Nðdbk ; s2 Þ. Thus, dbk is expressed in terms of the basic parameters, dbk 5 dAk dAb , while dAA is fixed at zero.4,7,21 A regression model using the first treatment as the reference treatment may be written as K P g y^ij 5 bj studyj 1 djbk tiA k5B
djAB ;NðdAB ; s2 Þ; . . . ; djAK ;NðdAK ; s2 Þ dAA 5 0
;
ðModel 4Þ
where y^ij is the estimated outcome for each arm i in study j, b1 to bp are regression coefficients for dummy variables study1 to studyp, respectively, and djAB to djAK are the regression coefficients for variables tAB to tAK, which use the contrast coding scheme. When trials with more than 2 arms are included in the network meta-analysis, djbk is no longer independent and the correlations must be accounted for in the estimation. As discussed in the previous section on randomeffects pairwise meta-analysis, random effects estimation is affected by covariate centering—and the basic parameter coding scheme is no exception. To estimate random-effects models where covariates are coded as those shown in Table 1 for tAB and tAC, one may use the following Stata code: melogit r tAB tAC study1-study26, noconstant || id: tAB tAC, noconstant covariance(identity) collinear binomial(n) /*code-3*/
The option noconstant just after the comma removes the intercept from the fixed-effects estimation. The statement after || specifies the random effects estimation in the model. The variable id is the trial identifier ranged from 1 to 26. Variables tAB and tAC specify random slopes. The option noconstant indicates that the intercept is not
ORIGINAL ARTICLE
915
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
TU
involved in the random effects estimation. The covariance structure is identity (i.e., the 2 random effects for slopes are independent with the same variance). The log OR for the A-B comparison is –0.641 (95% CI –1.145 to –0.136) while for the A-C comparison the log OR is –0.741 (95% CI –1.145 to –0.136). The variance of random effects is 0.888. One may then center tAB and tAC but keep the originals in fixed-effects estimation for easy interpretation: melogit r tAB tAC study1-study26, noconstant || id: zAB zAC, noconstant covariance(identity) collinear binomial(n) /*code-4*/
The zAB and zAC are centered tAB and tAC. Results show that the log OR for the A-B comparison is –0.606 (95% CI –1.135 to –0.076) and for the C-A comparison is –0.688 (95% CI –1.456 to 0.080). The variance of random effects is 1.013. However, this analysis so far does not take into account that 2 trials have 3 arms and, therefore, that their random effects are correlated. Lu and Ades suggested setting the correlation between random effects in multiarm trials at 0.5. However, the generalized linear mixed-model commands in Stata do not allow users to set the correlation between random effects at a fixed value. Therefore, one must use a ‘‘trick’’ to implement the fixed correlation: specify a 3-by-3 correlation matrix, M, with its diagonal element being unity and off-diagonal elements being 0.5. One then undertakes singular value decomposition for M (i.e., M = UDUT), as M is a symmetrical matrix. The D matrix is a diagonal matrix with 3 singular1 values of M. Next, create a new diagonal matrix, D2 , with the square root of the 3 singular values at its diagonal. 1 Finally, multiply U and D2 to obtain 3 new variables for random effects estimation xAA, xAB, and xAC: 1 X 5 ZUD2 , where the 3 columns of X are xAA, xAB, and xAC, and the 3 columns of Z are zAA, zAB, and zAC. The new model can be estimated using the follow Stata code: melogit r tAB tAC study1-study26, noconstant || id: xAB xAC, noconstant covariance(identity) collinear binomial(n) /*code-5*/
The log OR for the A-B comparison is –0.581 (95% CI –1.107 to –0.056) and the log OR for the A-C comparison is –0.740 (95% CI –1.490 to 0.010). The variance of the random effects is 1.026. These results are exactly the same as those reported by Jones and others7 in their study on the use of the frequentist approach to network meta-analysis. Their model,
916
however, used the following dummy coding scheme: g y^ij 5 bj studyj 1 gij Ti ; i ¼ B; :::; k; j ¼ 1; :::; p gBj ;NðcB ; s2 Þ; . . . ; gkj ;Nðck ; s2 Þ
:
ðModel 5Þ
Jones and others also noted that the results of their model changed when different treatments were used as the primary reference. They thought that the nonlinearity of the link function led to the different estimates for treatment differences. However, I think the differences are caused by the specification of random effects estimation. Jones and others proposed a ‘‘symmetric’’ constraint for coding the covariates used for random effects estimation (zA, zB, and zC in Table 1 of the online appendix).7 Their approach creates exactly the same 3 covariates for random effects estimation as centering does (i.e., zA and zAA, zB and zAB, and zC and zAC are the same). 1 (the new They then multiply zA, zB, and zC by =2 variables are named xA, xB, and xC, respectively) to create a correlation of 0.5 among the 3 random effects: melogit r tB tC study1-study26, noconstant || id: xA xB xC, noconstant covariance(identity) collinear binomial(n) /*code-6*/
The above Stata code yields the same fixed-effects estimates, and the variance for random effects (i.e., s2 in Model 4) is 1.026. If one substitutes zA, zB, and zC for xA, xB, and xC in the above Stata code, the new model yields the same fixed-effects estimates, but the variance for random effects becomes 0.513, which is about half of 1.026. As zA, zB, and zC estimate random effects for treatments A, B, and C, respectively, the estimated random effects should be independent. The variance of the random effects for the difference between 2 treatments is the sum of the variances of random effects for each treatment, when the 2 treatments are independent. The variances for zA, zB, and zC, which are assumed to be equal, are therefore half of those for xA, xB, and xC. Nevertheless, estimating the variance of random effects for each treatment is more straightforward than estimating the dependent random effects for the difference in treatment effects between 2 treatments. Complete Stata codes and output for all the analyses and SAS (SAS Institute, Cary, North Carolina) codes for network meta-analysis are available in the online appendix.
MEDICAL DECISION MAKING/OCTOBER 2014
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
GLMM FOR NETWORK META-ANALYSIS
DISCUSSION Although the Bayesian approach has many attractive features, it is useful to implement the Lu and Ades model under the frequentist statistical framework, allowing researchers to gain a better understanding of the complexity and intricacy of network meta-analysis. It is also straightforward to incorporate covariates, such as moderators and confounders, into the generalized linear mixed model to conduct meta-regression for multiple treatment comparisons. Moreover, this approach may be extended easily to other types of outcome variables, such as continuous, counts, and multinomial. The statistical models discussed in this article implement treatment arm–based network metaanalysis. Each row in the dataset represents a treatment arm in a trial. Some proposed network metaanalyses use a trial-based analysis strategy, with each data entry representing a treatment difference.2,12 For multiarm trials, more than 1 entry will be needed; the correlation between the random effects must be dealt with carefully. However, because calculating treatment differences usually relies on closed formulas for noncontinuous outcomes, these formulas may not work well when the outcome event is rare. In contrast, our treatment arm–based approach uses the exact likelihood for statistical modeling and, therefore, is more accurate. For instance, the Stata command mvmeta for the multivariate meta-analysis has been used for network meta-analysis.12,29,30 First, a reference treatment A is selected, and data augmentation is undertaken for the trial without treatment A. Trials containing 0 cells have 1 individual with 0.5 success added to each arm.29 Then, pairwise differences in treatment effects between other treatments and A within the same trial are calculated for network meta-analysis. For multiarm trials, the within-trial correlation between pairs of treatment comparisons can be set up. Using the sample dataset that compared 3 treatment options for prevention of first bleeding in patients with liver cirrhosis, results from mvmeta show that the log ORs for the A-B and A-C comparisons are –0.584 (95% CI –1.110 to –0.056) and –0.716 (95% CI –1.476 to 0.043), respectively. The variance of random effects is 0.996. Those small differences in results are likely due to the use of a closed formula for the variance of log OR and continuity correction for zero cells. Recently, Biondi-Zoccai and D’Ascenzo31 discussed differences and similarities between frequentist and Bayesian network meta-analysis as well as the pros and cons of various statistical packages. In Table
3 of the online appendix, I summarize the pros and cons of frequentist approaches to network metaanalysis by making the distinction between armand trial-based approaches. ACKNOWLEDGMENTS The author thanks the 2 anonymous reviewers for their helpful comments.
REFERENCES 1. Sutton AJ, Higgins JPT. Recent developments in meta-analysis. Stat Med. 2008;27(5):625–50. 2. Lumley T. Network meta-analysis for indirect treatment comparisons. Stat Med. 2002;21(16):2313–24. 3. Lu G, Welton NJ, Higgins JPT, White IR, Ades AE. Linear inference for mixed treatment comparison meta-analysis: a two-stage approach. Res Synth Methods. 2011;2(1):43–60. 4. Dias S, Sutton AJ, Ades AE, Welton NJ. Evidence synthesis for decision making 2: a generalized linear modeling framework for pairwise and network meta-analysis of randomized controlled trials. Med Decis Making. 2013;33(5):607–17. 5. Jansen JP, Fleurence R, Devine B, et al. Interpreting indirect treatment comparisons and network meta-analysis for health-care decision making: report of the ISPOR task force on indirect treatment comparisons good research practices: part 1. Value Health. 2011;14(4):417–28. 6. Lu G, Ades AE. Combination of direct and indirect evidence in mixed treatment comparisons. Stat Med. 2004;23(20):3105–24. 7. Jones B, Roger J, Lane PW, et al. Statistical approaches for conducting network meta-analysis in drug development. Pharm Stat. 2011;10(6):523–31. 8. Song F, Altman DG, Glenny A-M, Deeks JJ. Validity of indirect comparison for estimating efficacy of competing interventions: empirical evidence from published meta-analyses. BMJ. 2003; 326(7387):472. 9. Song F, Loke YK, Walsh T, Glenny A-M, Eastwood AJ, Altman DG. Methodological problems in the use of indirect comparisons for evaluating healthcare interventions: survey of published systematic reviews. BMJ. 2009;338(7700):b1147. 10. Glenny A-M, Altman D, Song F, et al. Indirect Comparisons of Competing Interventions. Health Technology Assessment. 2005; 9(26):1–148. 11. Dias S, Welton NJ, Caldwell DM, Ades AE. Checking consistency in mixed treatment comparison meta-analysis. Stat Med. 2010;29(7–8):932–44. 12. Higgins JPT, Jackson D, Barrett JK, Lu G, Ades AE, White IR. Consistency and inconsistency in network meta-analysis: concepts and models for multi-arm studies. Res Synth Methods. 2012;3(2): 98–110. 13. Dias S, Welton NJ, Sutton AJ, Caldwell DM, Lu G, Ades AE. Evidence synthesis for decision making 4: inconsistency in networks of evidence based on randomized controlled trials. Med Decis Making. 2013;33(5):641–56.
ORIGINAL ARTICLE
917
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
TU
14. Bucher HC, Guyatt GH, Griffith LE, Walter SD. The results of direct and indirect treatment comparisons in meta-analysis of randomized controlled trials. J Clin Epidemiol. 1997;50(6):683–91. 15. Hasselblad V. Meta-analysis of multitreatment studies. Med Decis Making. 1998;18(1):37–43. 16. Donegan S, Williamson P, D’Alessandro U, Tudur Smith C. Assessing key assumptions of network meta-analysis: a review of methods. Res Synth Methods. 2013;4(4):291–323. 17. Senn S, Gavini F, Magrez D, Scheen A. Issues in performing a network meta-analysis. Stat Methods Med Res. 2013;22(2): 169–89. 18. Cipriani A, Higgins JPT, Geddes JR, Salanti G. Conceptual and technical challenges in network meta-analysis. Ann Intern Med. 2013;159(2):130–7. 19. Lee AW. Review of mixed treatment comparisons in published systematic reviews shows marked increase since 2009. J Clin Epidemiol. 2014;67(2):138–43. 20. Hong H, Carlin BP, Shamliyan TA, et al. Comparing Bayesian and frequentist approaches for multiple outcome mixed treatment comparisons. Med Decis Making. 2013;33(5):702–14. 21. Hoaglin DC, Hawkins N, Jansen JP, et al. Conducting indirecttreatment-comparison and network-meta-analysis studies: report of the ISPOR task force on indirect treatment comparisons good research practices: part 2. Value Health. 2011;14(4):429–37. 22. Senn S. The many modes of meta. Drug Inf J. 2000;34(2): 535–49. 23. Colditz GA, Brewer TF, Berkey CS, et al. Efficacy of BCG vaccine in the prevention of tuberculosis: meta-analysis of the published literature. JAMA. 1994;271(9):698–702.
918
24. Borenstein M, Hedges LV, Higgins JPT, Rothstein HR. Meta-regression. In: Borenstein M, Hedges LV, Higgins JPT, Rothstein HR, eds. Introduction to Meta-Analysis. New York: John Wiley & Sons; 2009. p 187–203. 25. Harris R, Bradburn M, Deeks J, Harbord R, Altman D, Sterne J. metan: fixed- and random-effects meta-analysis. Stata J. 2008;8(1): 3–28. 26. Smith TC, Spiegelhalter DJ, Thomas A. Bayesian approaches to random-effects meta-analysis: a comparative study. Stat Med. 1995;14(24):2685–99. 27. Piepho HP, Williams ER, Madden LV. The use of two-way linear mixed models in multitreatment meta-analysis. Biometrics. 2012;68(4):1269–77. 28. Pagliaro L, D’Amico G, So¨rensen TIA, et al. Prevention of first bleeding in cirrhosis: a meta-analysis of randomized trials of nonsurgical treatment. Ann Intern Med. 1992;117(1):59–70. 29. White IR. Multivariate random-effects meta-regression: updates to mvmeta. Stata J. 2011;11(2):255–70. 30. White IR, Barrett JK, Jackson D, Higgins JPT. Consistency and inconsistency in network meta-analysis: model estimation using multivariate meta-regression. Res Synth Methods. 2012;3(2): 111–25. 31. Biondi-Zoccai G, D’Ascenzo F. Choosing between frequentist and Bayesian frameworks and the corresponding statistical package. In: Biondi-Zoccai G, ed. Network Meta-Analysis: Evidence Synthesis with Mixed Treatment Comparison. Hauppauge, NY: Nova Science Publishers; 2014. p 101–13.
MEDICAL DECISION MAKING/OCTOBER 2014
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
Use of Generalized Linear Mixed Models for Network Meta-analysis Yu-Kang Tu, PhD
In the past decade, a new statistical method—network meta-analysis—has been developed to address limitations in traditional pairwise meta-analysis. Network metaanalysis incorporates all available evidence into a general statistical framework for comparisons of multiple treatments. Bayesian network meta-analysis, as proposed by Lu and Ades, also known as ‘‘mixed treatments comparisons,’’ provides a flexible modeling framework to take into account complexity in the data structure. This article shows how to implement the Lu and Ades model in the frequentist generalized linear mixed model. Two examples are provided to demonstrate how centering the covariates for random effects estimation within each trial can yield
correct estimation of random effects. Moreover, under the correct specification for random effects estimation, the dummy coding and contrast basic parameter coding schemes will yield the same results. It is straightforward to incorporate covariates, such as moderators and confounders, into the generalized linear mixed model to conduct meta-regression for multiple treatment comparisons. Moreover, this approach may be extended easily to other types of outcome variables, such as continuous, counts, and multinomial. Key words: randomized controlled trials; network meta-analysis; mixed treatments comparisons; generalized linear mixed models. (Med Decis Making 2014;34:911–918)
M
comparison. The evidence base, therefore, varies across comparisons. One potential consequence is that the results from multiple pairwise meta-analyses are not always consistent, and some comparisons may not have been conducted yet.11–13 Network meta-analysis incorporates all available evidence into a general statistical framework for comparisons of all available treatments to resolve the limitations of traditional pairwise meta-analysis. The idea of indirect comparisons was first proposed in the 1990s.14,15 Lumley2 coined the term ‘‘network meta-analysis’’ in 2002 and proposed a linear mixed-model approach to comparisons of multiple treatments within the same statistical model. Lu and Ades6 developed a sophisticated Bayesian hierarchical model, providing a flexible statistical framework to take into account the complexity of the data structure within multiarm trials. Their statistical approach, widely known as mixed treatments comparison or Bayesian network meta-analysis, is currently the standard approach to comparisons of multiple treatments.7,16–18 Comparing all available treatments simultaneously is an attractive idea to clinicians and health care providers. The number of articles that report the results of network meta-analysis has increased exponentially in the past few years, as have their citations,18,19 despite the fact that implementing
eta-analysis is an important research tool for synthesis of evidence on the efficacy versus harm of available medical interventions.1 One recent development in meta-analysis methodology is network meta-analysis for comparisons of multiple treatments.2–10 When several treatments are available for a medical condition, traditional metaanalysis undertakes multiple pairwise comparisons, using different sets of clinical trials for each
Received 24 April 2014 from the Institute of Epidemiology & Preventive Medicine, College of Public Health, National Taiwan University, Taipei, Taiwan (Y-KT). The author declares no conflict of interests. This project was partly funded by a grant from the National Science Council in Taiwan (Grant No. NSC 101-2314-B-002-197-MY2). Revision accepted for publication 30 June 2014. Supplementary material for this article is available on the Medical Decision Making Web site at http://mdm.sagepub.com/supplemental. Address correspondence to Yu-Kang Tu, Institute of Epidemiology & Preventive Medicine, College of Public Health, National Taiwan University, 17 Xu-Zhou Road, Taipei, Taiwan; telephone: (886) 2 3366 8039; fax: (886) 2351 1955; e-mail: [email protected]. Ó The Author(s) 2014 Reprints and permission: http://www.sagepub.com/journalsPermissions.nav DOI: 10.1177/0272989X14545789
MEDICAL DECISION MAKING/OCTOBER 2014
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
911
TU
Bayesian network meta-analysis is quite challenging. Several recent articles have attempted to elaborate further on the statistical issues related to implementing the Lu and Ades approach within the generalized linear mixed-modeling framework using nonBayesian statistical software packages,3,7,17,20 with the intention of making network meta-analysis more accessible to clinicians and meta-analysts who are not yet familiar with Bayesian statistics. This article provides a general, non-Bayesian approach to network meta-analysis based on the Lu and Ades model. I first demonstrate how a pairwise meta-analysis may be implemented in a generalized linear mixed model and show how different parameterizations may affect results from random-effects models. I then extend the discussion to multiple treatment comparisons and demonstrate how the Lu and Ades model may be implemented in a generalized linear mixed model with an appropriate parameterization of random effects estimation. Finally, I provide an example to show that 2 proposed generalized linear mixed-modeling approaches to network metaanalysis are different parameterizations of the same model. GENERALIZED LINEAR MIXED MODEL FOR PAIRWISE META-ANALYSIS Fixed-effects Model for Pairwise Meta-analysis A fixed-effects model assumes a common effect behind the observed effects, while a random-effects model assumes that the ‘‘true’’ effect follows a distribution. The basic fixed-effects model for comparing multiple treatments is as follows21:
g yij 5 sj 1 ti ;
ðEquation 1Þ
where yij is the outcome for treatment arm i in study j; sj is the study effect; and ti is the treatment effect. The function g(.) is a suitable link function for transformation of yij. For example, when the outcome is probability of an event for a subject receiving treatment i in study j, the link function is logit. In a traditional pairwise meta-analysis with p studies that compare treatment B to treatment A, the generalized linear regression model for Equation 1 may be expressed as follows22: g y^ij 5 b0 1 bj studyj 1 bp Ti ; i ¼ A; B; j ¼ 1; :::; p 1
ðModel 1Þ
912
where y^ij is the estimated outcome for each arm i in study j, b0 is the intercept, and b1 to bp–1 are regression coefficients for study1 to studyp–1, respectively, while bp is the regression coefficient for dummy variable Ti, which is coded as 1 for treatment i = B and 0 for treatment i = A. Consequently, b0 is the estimated treatment effect for A in study p; and bj, where j = 1 to p – 1, is the difference in treatment effects between either treatment i in study j and the same treatment i in study p. The regression coefficient, bp, is the average difference in treatment effects between B and A across all trials. The Lu and Ades model for fixed-effects metaanalysis comparing 2 treatments may be written as follows4,21: hkj 5
mj k 5 A ; mj 1 d k 5 B
ðEquation 2Þ
where hkj is the estimated outcome for treatment k in study j—that is, gð^ ykj Þ in Model 1. Consequently, d is the difference in treatment effects between A and B. The subscript j for m indicates that the treatment effect for A is estimated separately for each study, so Equation 2 and Equation 1 are equivalent and may be implemented in Model 1. Random-effects Model for Pairwise Meta-analysis Equation 2 assumes that the difference in treatment effects between A and B is a fixed value, d, which is estimated by bp in Model 1. In a randomeffects models, the difference in treatment effects is no longer a fixed value but rather a random variable that follows a normal distribution4,21:
k 5A mj mj 1 dj k 5 B dj ;Nðd; s2 Þ: hkj 5
ðEquation 3Þ
The mean of dj is d and its variance is s2 (i.e., the difference in treatment effects between A and B is a random variable). The generalized linear mixed regression model for Equation 3 may be written as g y^ij 5 b0 1 bj studyj 1 bpj Ti ; i ¼ A; B; j ¼ 1; :::; p 1 bpj 5 bp 1 uj ;
ðModel 2Þ
where y^ij is the estimated treatment effect for each arm i in study j, and b0 is the intercept, b1 to bp–1 are regression coefficients for study1 to studyp–1, respectively, and bpj is the regression coefficient for
MEDICAL DECISION MAKING/OCTOBER 2014
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
GLMM FOR NETWORK META-ANALYSIS
dummy variable Ti in study j. The average difference in treatment effects is bp , and uj is the random effect for bp (i.e., the difference between bp and bpj Þ. However, Model 2 needs some fine-tuning before it will yield a correct estimation for bp and the variance of uj, as the following example illustrates. Colditz and others23 provided data from 13 randomized controlled trials on the efficacy of Bacillus Calmette–Gue´rin (BCG) vaccine in the prevention of tuberculosis. Applying data presented by Borenstein and others24 and using the command metan for the statistical software Stata (College Station, Texas) will allow us to undertake DerSimonian and Laird random-effects meta-analysis.25 The pooled odds ratio (OR) is 0.474 (95% confidence interval [CI] 0.325–0.691). The between-study variance, usually denoted as t2 , is 0.3682. To undertake the random-effects analysis using a generalized linear mixed model, such as Model 2, one may use the Stata command melogit or meqrlogit in version 13 or xtmelogit in version 12. The following Stata code used melogit: melogit r bcg study1-study12, || trial: bcg, noconstant binomial(n) or
Each study has 2 data entries: 1 data entry for the BCG-vaccinated group and 1 data entry for the control group. Variable r is the number of tuberculosis cases; n is the total number of subjects in each group. Variable trial was the identifier for each trial, and study1 to study12 are 12 dummy variables for trial 1 to trial 12, respectively. The option noconstant indicates that the intercept is excluded from the estimation of random effects (i.e., this model only estimated random effects for the slope of variable bcg, which is a dummy variable with the vaccinated group coded as 1 and the control group coded as 0). The difference in treatment effects among the control groups from different trials is estimated as fixed effects by the dummy variables study1 to study12. The variations in the treatment effects of control groups may be estimated as a random variable, but, as Senn22 noted, this is not the standard approach in meta-analysis. The regression coefficient for bcg is –0.764 (OR 0.466; 95% CI 0.332–0.654), which is smaller than that given by metan; s2 , the estimated variance for random effects, is 0.2811. To obtain a similar OR, one must center bcg (i.e., recode the vaccinated group as 0.5; control group, –0.5). The results then show that the OR for centered bcg is 0.475 (95% CI 0.337–0.770) and s2 is 0.295. Note that for the fixed-effects model, centering would not cause a difference in the OR.
Figure 1 illustrates how centering causes the differences in random effects estimation. When the control group is coded 0 and the BCG-vaccinated group is coded 1, random effects, as a random slope for bcg, are applied to the vaccinated group only (Figure 1A); when bcg is centered, random effects are applied to both groups (i.e., to the difference between vaccinated and control groups, Figure 1B).7 Consequently, the random slope for centered bcg can now rotate around the midpoint between the vaccinated and control groups, so the difference between the 2 groups may be estimated as random effects under the constraint that the sum of the random effects within a trial is zero.7 The centered bcg therefore gives the model greater ‘‘freedom’’ to estimate the differences between the 2 groups and should yield a better model fit: The log likelihoods for uncentered and centered models are –90.50 and –88.64, respectively. This approach with centered bcg is similar to the Bayesian meta-analysis model proposed by Smith and others26: di 5 logit pTi logit pCi T mi 5 ðlogit pi 1 logit pCi Þ=2 ; logit pCi 5 mi 12 di T 1 logit pi 5 mi 1 di 2
where pTi and pC i are the event probability for treatment and control groups. Lu and Ades6 extended this model to multiple treatment comparisons. Centering the covariates for random effects estimation is crucial, especially when the heterogeneity across trials is large. For traditional meta-analysis, the 2 groups are coded – 0.5 and 0.5 after centering. For network meta-analysis with multiarm trials, the mean value of a dummy variable is subtracted from its individual value. For instance, in a 3-arm trial, treatments A, B, and C would be coded 1, 0, and 0, respectively, in the dummy variable for treatment A. However, after centering, they would instead be coded 2/3, –1/3, and –1/3, respectively. GENERALIZED LINEAR MIXED MODEL FOR NETWORK META-ANALYSIS Fixed-effects Model for Network Meta-analysis We must now extend our generalized linear mixed model to multiple treatment comparisons. For the fixed-effects model, the basic equations are the same as those for Model 1, with more treatments to compare. The regression model may be written as
ORIGINAL ARTICLE
913
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
TU
Figure 1 Diagram for the impact of centering on random effects estimation. (A) Dummy coding for variable bcg where the control group is coded 0 and the vaccinated group coded 1. (B) Centered bcg where the control group is coded –0.5 and the vaccinated group coded 0.5.
g y^ij 5 b0 1 bj studyj þ ci Ti ; i ¼ 1; :::; k; j ¼ 1; :::; p: ðModel 3Þ
where y^ij is the estimated outcome for each arm i in study j, b0 is the intercept, and b1 to bp are regression coefficients for dummy variables study1 to studyp, respectively. Furthermore, c1 to cK are the regression coefficients for dummy variables T1 to TK, where treatment k is coded 1 and the other arms are coded 0 for Tk. Suppose one uses the first treatment and the first trial as the primary reference groups, so b1 and c1 are fixed at 0. Consequently, b0 is the estimated treatment effect for treatment 1 in study 1; and bj, where j = 2 to p, is the difference in treatment effects between treatment i in study j and treatment i in study 1. The regression coefficient ci is the difference in treatment effects between treatment i and treatment 1 across all trials. The fixed-effects Lu and Ades model may be specified as 8 mjb > > < b 5 A; B; C; . . . if k 5 b hjk 5 m 1 d 5 m 1 d d ; > bk Ak Ab jb > : jb 0 k 5 B; C; D; . . . if k is 0 after B
ðEquation 4Þ
where hjk is the estimated outcome for arm k in study j, and mjb is the baseline treatment in trial j. The difference between the other treatment k and treatment b in the same trial will be estimated by expressing them in
914
terms of effects relative to treatment A, which is the primary reference treatment within the network. For reasons of identification and because of its interpretation as the effect of treatment A compared with itself, dAA is fixed at 0, and Lu and Ades called dAB to dAk the ‘‘basic parameters.’’ The advantage of expressing all treatment comparisons as the relation between basic parameters is that the number of pairwise comparisons to be estimated for a network meta-analysis involving k treatments is reduced to k 2 1 for the fixed effects.27 For instance, suppose one has 3 trials. Trial 1 compares treatments A, B and C, while trial 2 compares A and B, and trial 3 compares B and C. Table 1 shows the difference in coding the covariates for treatments. Variables tA, tB, and tC are dummy variables for the 3 treatments, and because A is the primary reference treatment, only tB and tC are included in the model. Variables tAA tAB, and tAC use the contrast coding for basic parameters.27 Furthermore, because the regression coefficient for tAA is fixed at 0, only tAB and tAC are included in the model. Both coding schemes will yield the same results for the fixedeffects model. If we exclude the intercept and include all 3 dummy variables for the trials (i.e., trial 1, trial 2, and trial 3 in Table 1) in the model, the basic parameter coding has 1 advantage, which is that the dummy variables for each trial have a clear interpretation: They represent the trial-specific baseline treatments for each trial.
MEDICAL DECISION MAKING/OCTOBER 2014
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
GLMM FOR NETWORK META-ANALYSIS
Table 1
Two Coding Schemes for Covariates Included in the Generalized Linear Mixed Model for Treatment Effects Estimation Dummy Coding
Trial
1 1 1 2 2 3 3
Contrast Coding
Study 1
Study 2
Study 3
Treatment
Baseline
tA
tB
tC
tAA
tAB
tAC
1 1 1 0 0 0 0
0 0 0 1 1 0 0
0 0 0 0 0 1 1
A B C A B B C
A A A A A B B
1 0 0 1 0 0 0
0 1 0 0 1 1 0
0 0 1 0 0 0 1
0 –1 –1 0 –1 0 0
0 1 0 0 1 0 –1
0 0 1 0 0 0 1
Sample Data The sample dataset contains the results of 26 studies that directly compared 3 treatments for the prevention of first bleeding in patients with liver cirrhosis28: Treatment A was the control, nonactive treatment; B, sclerotherapy; and C, beta-blocker. (The dataset is available in Table 1 of the online appendix.) Among the 26 studies, 2 studies were 3arm trials, 7 studies compared A with C, and 17 studies compared B with C. The following Stata code undertakes fixed-effects network meta-analysis using the contrast coding scheme for treatments: melogit r tAB tAC study1-study26, noconstant binomial(n) /*code-1*/
The option noconstant removes the intercept from the model, so all the dummy variables, study1 to study26, for estimating the effects of the baseline treatment for the 26 trials can be included. The variable r is the number of events (i.e., bleeding); n is the number of patients in each arm. Variables tAB and tAC use the contrast coding scheme for the differences between B and A and between C and A, respectively, such as tAB and tAC in Table 1. The same results may also be obtained using the following dummy coding scheme: melogit r tB tC study1-study26, noconstant binomial(n) /*code-2*/
Variables tB and tC are dummy variables for treatments B and C, shown as tB and tC in Table 1. Random-effects Model for Network Meta-analysis For the random-effects network meta-analysis, dbk in Equation 4 is replaced by djbk , the trial-specific effect of treatment k relative to the trial-specific baseline treatment b. These trial-specific effects are then
drawn from a normal distribution: djbk ;Nðdbk ; s2 Þ. Thus, dbk is expressed in terms of the basic parameters, dbk 5 dAk dAb , while dAA is fixed at zero.4,7,21 A regression model using the first treatment as the reference treatment may be written as K P g y^ij 5 bj studyj 1 djbk tiA k5B
djAB ;NðdAB ; s2 Þ; . . . ; djAK ;NðdAK ; s2 Þ dAA 5 0
;
ðModel 4Þ
where y^ij is the estimated outcome for each arm i in study j, b1 to bp are regression coefficients for dummy variables study1 to studyp, respectively, and djAB to djAK are the regression coefficients for variables tAB to tAK, which use the contrast coding scheme. When trials with more than 2 arms are included in the network meta-analysis, djbk is no longer independent and the correlations must be accounted for in the estimation. As discussed in the previous section on randomeffects pairwise meta-analysis, random effects estimation is affected by covariate centering—and the basic parameter coding scheme is no exception. To estimate random-effects models where covariates are coded as those shown in Table 1 for tAB and tAC, one may use the following Stata code: melogit r tAB tAC study1-study26, noconstant || id: tAB tAC, noconstant covariance(identity) collinear binomial(n) /*code-3*/
The option noconstant just after the comma removes the intercept from the fixed-effects estimation. The statement after || specifies the random effects estimation in the model. The variable id is the trial identifier ranged from 1 to 26. Variables tAB and tAC specify random slopes. The option noconstant indicates that the intercept is not
ORIGINAL ARTICLE
915
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
TU
involved in the random effects estimation. The covariance structure is identity (i.e., the 2 random effects for slopes are independent with the same variance). The log OR for the A-B comparison is –0.641 (95% CI –1.145 to –0.136) while for the A-C comparison the log OR is –0.741 (95% CI –1.145 to –0.136). The variance of random effects is 0.888. One may then center tAB and tAC but keep the originals in fixed-effects estimation for easy interpretation: melogit r tAB tAC study1-study26, noconstant || id: zAB zAC, noconstant covariance(identity) collinear binomial(n) /*code-4*/
The zAB and zAC are centered tAB and tAC. Results show that the log OR for the A-B comparison is –0.606 (95% CI –1.135 to –0.076) and for the C-A comparison is –0.688 (95% CI –1.456 to 0.080). The variance of random effects is 1.013. However, this analysis so far does not take into account that 2 trials have 3 arms and, therefore, that their random effects are correlated. Lu and Ades suggested setting the correlation between random effects in multiarm trials at 0.5. However, the generalized linear mixed-model commands in Stata do not allow users to set the correlation between random effects at a fixed value. Therefore, one must use a ‘‘trick’’ to implement the fixed correlation: specify a 3-by-3 correlation matrix, M, with its diagonal element being unity and off-diagonal elements being 0.5. One then undertakes singular value decomposition for M (i.e., M = UDUT), as M is a symmetrical matrix. The D matrix is a diagonal matrix with 3 singular1 values of M. Next, create a new diagonal matrix, D2 , with the square root of the 3 singular values at its diagonal. 1 Finally, multiply U and D2 to obtain 3 new variables for random effects estimation xAA, xAB, and xAC: 1 X 5 ZUD2 , where the 3 columns of X are xAA, xAB, and xAC, and the 3 columns of Z are zAA, zAB, and zAC. The new model can be estimated using the follow Stata code: melogit r tAB tAC study1-study26, noconstant || id: xAB xAC, noconstant covariance(identity) collinear binomial(n) /*code-5*/
The log OR for the A-B comparison is –0.581 (95% CI –1.107 to –0.056) and the log OR for the A-C comparison is –0.740 (95% CI –1.490 to 0.010). The variance of the random effects is 1.026. These results are exactly the same as those reported by Jones and others7 in their study on the use of the frequentist approach to network meta-analysis. Their model,
916
however, used the following dummy coding scheme: g y^ij 5 bj studyj 1 gij Ti ; i ¼ B; :::; k; j ¼ 1; :::; p gBj ;NðcB ; s2 Þ; . . . ; gkj ;Nðck ; s2 Þ
:
ðModel 5Þ
Jones and others also noted that the results of their model changed when different treatments were used as the primary reference. They thought that the nonlinearity of the link function led to the different estimates for treatment differences. However, I think the differences are caused by the specification of random effects estimation. Jones and others proposed a ‘‘symmetric’’ constraint for coding the covariates used for random effects estimation (zA, zB, and zC in Table 1 of the online appendix).7 Their approach creates exactly the same 3 covariates for random effects estimation as centering does (i.e., zA and zAA, zB and zAB, and zC and zAC are the same). 1 (the new They then multiply zA, zB, and zC by =2 variables are named xA, xB, and xC, respectively) to create a correlation of 0.5 among the 3 random effects: melogit r tB tC study1-study26, noconstant || id: xA xB xC, noconstant covariance(identity) collinear binomial(n) /*code-6*/
The above Stata code yields the same fixed-effects estimates, and the variance for random effects (i.e., s2 in Model 4) is 1.026. If one substitutes zA, zB, and zC for xA, xB, and xC in the above Stata code, the new model yields the same fixed-effects estimates, but the variance for random effects becomes 0.513, which is about half of 1.026. As zA, zB, and zC estimate random effects for treatments A, B, and C, respectively, the estimated random effects should be independent. The variance of the random effects for the difference between 2 treatments is the sum of the variances of random effects for each treatment, when the 2 treatments are independent. The variances for zA, zB, and zC, which are assumed to be equal, are therefore half of those for xA, xB, and xC. Nevertheless, estimating the variance of random effects for each treatment is more straightforward than estimating the dependent random effects for the difference in treatment effects between 2 treatments. Complete Stata codes and output for all the analyses and SAS (SAS Institute, Cary, North Carolina) codes for network meta-analysis are available in the online appendix.
MEDICAL DECISION MAKING/OCTOBER 2014
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
GLMM FOR NETWORK META-ANALYSIS
DISCUSSION Although the Bayesian approach has many attractive features, it is useful to implement the Lu and Ades model under the frequentist statistical framework, allowing researchers to gain a better understanding of the complexity and intricacy of network meta-analysis. It is also straightforward to incorporate covariates, such as moderators and confounders, into the generalized linear mixed model to conduct meta-regression for multiple treatment comparisons. Moreover, this approach may be extended easily to other types of outcome variables, such as continuous, counts, and multinomial. The statistical models discussed in this article implement treatment arm–based network metaanalysis. Each row in the dataset represents a treatment arm in a trial. Some proposed network metaanalyses use a trial-based analysis strategy, with each data entry representing a treatment difference.2,12 For multiarm trials, more than 1 entry will be needed; the correlation between the random effects must be dealt with carefully. However, because calculating treatment differences usually relies on closed formulas for noncontinuous outcomes, these formulas may not work well when the outcome event is rare. In contrast, our treatment arm–based approach uses the exact likelihood for statistical modeling and, therefore, is more accurate. For instance, the Stata command mvmeta for the multivariate meta-analysis has been used for network meta-analysis.12,29,30 First, a reference treatment A is selected, and data augmentation is undertaken for the trial without treatment A. Trials containing 0 cells have 1 individual with 0.5 success added to each arm.29 Then, pairwise differences in treatment effects between other treatments and A within the same trial are calculated for network meta-analysis. For multiarm trials, the within-trial correlation between pairs of treatment comparisons can be set up. Using the sample dataset that compared 3 treatment options for prevention of first bleeding in patients with liver cirrhosis, results from mvmeta show that the log ORs for the A-B and A-C comparisons are –0.584 (95% CI –1.110 to –0.056) and –0.716 (95% CI –1.476 to 0.043), respectively. The variance of random effects is 0.996. Those small differences in results are likely due to the use of a closed formula for the variance of log OR and continuity correction for zero cells. Recently, Biondi-Zoccai and D’Ascenzo31 discussed differences and similarities between frequentist and Bayesian network meta-analysis as well as the pros and cons of various statistical packages. In Table
3 of the online appendix, I summarize the pros and cons of frequentist approaches to network metaanalysis by making the distinction between armand trial-based approaches. ACKNOWLEDGMENTS The author thanks the 2 anonymous reviewers for their helpful comments.
REFERENCES 1. Sutton AJ, Higgins JPT. Recent developments in meta-analysis. Stat Med. 2008;27(5):625–50. 2. Lumley T. Network meta-analysis for indirect treatment comparisons. Stat Med. 2002;21(16):2313–24. 3. Lu G, Welton NJ, Higgins JPT, White IR, Ades AE. Linear inference for mixed treatment comparison meta-analysis: a two-stage approach. Res Synth Methods. 2011;2(1):43–60. 4. Dias S, Sutton AJ, Ades AE, Welton NJ. Evidence synthesis for decision making 2: a generalized linear modeling framework for pairwise and network meta-analysis of randomized controlled trials. Med Decis Making. 2013;33(5):607–17. 5. Jansen JP, Fleurence R, Devine B, et al. Interpreting indirect treatment comparisons and network meta-analysis for health-care decision making: report of the ISPOR task force on indirect treatment comparisons good research practices: part 1. Value Health. 2011;14(4):417–28. 6. Lu G, Ades AE. Combination of direct and indirect evidence in mixed treatment comparisons. Stat Med. 2004;23(20):3105–24. 7. Jones B, Roger J, Lane PW, et al. Statistical approaches for conducting network meta-analysis in drug development. Pharm Stat. 2011;10(6):523–31. 8. Song F, Altman DG, Glenny A-M, Deeks JJ. Validity of indirect comparison for estimating efficacy of competing interventions: empirical evidence from published meta-analyses. BMJ. 2003; 326(7387):472. 9. Song F, Loke YK, Walsh T, Glenny A-M, Eastwood AJ, Altman DG. Methodological problems in the use of indirect comparisons for evaluating healthcare interventions: survey of published systematic reviews. BMJ. 2009;338(7700):b1147. 10. Glenny A-M, Altman D, Song F, et al. Indirect Comparisons of Competing Interventions. Health Technology Assessment. 2005; 9(26):1–148. 11. Dias S, Welton NJ, Caldwell DM, Ades AE. Checking consistency in mixed treatment comparison meta-analysis. Stat Med. 2010;29(7–8):932–44. 12. Higgins JPT, Jackson D, Barrett JK, Lu G, Ades AE, White IR. Consistency and inconsistency in network meta-analysis: concepts and models for multi-arm studies. Res Synth Methods. 2012;3(2): 98–110. 13. Dias S, Welton NJ, Sutton AJ, Caldwell DM, Lu G, Ades AE. Evidence synthesis for decision making 4: inconsistency in networks of evidence based on randomized controlled trials. Med Decis Making. 2013;33(5):641–56.
ORIGINAL ARTICLE
917
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
TU
14. Bucher HC, Guyatt GH, Griffith LE, Walter SD. The results of direct and indirect treatment comparisons in meta-analysis of randomized controlled trials. J Clin Epidemiol. 1997;50(6):683–91. 15. Hasselblad V. Meta-analysis of multitreatment studies. Med Decis Making. 1998;18(1):37–43. 16. Donegan S, Williamson P, D’Alessandro U, Tudur Smith C. Assessing key assumptions of network meta-analysis: a review of methods. Res Synth Methods. 2013;4(4):291–323. 17. Senn S, Gavini F, Magrez D, Scheen A. Issues in performing a network meta-analysis. Stat Methods Med Res. 2013;22(2): 169–89. 18. Cipriani A, Higgins JPT, Geddes JR, Salanti G. Conceptual and technical challenges in network meta-analysis. Ann Intern Med. 2013;159(2):130–7. 19. Lee AW. Review of mixed treatment comparisons in published systematic reviews shows marked increase since 2009. J Clin Epidemiol. 2014;67(2):138–43. 20. Hong H, Carlin BP, Shamliyan TA, et al. Comparing Bayesian and frequentist approaches for multiple outcome mixed treatment comparisons. Med Decis Making. 2013;33(5):702–14. 21. Hoaglin DC, Hawkins N, Jansen JP, et al. Conducting indirecttreatment-comparison and network-meta-analysis studies: report of the ISPOR task force on indirect treatment comparisons good research practices: part 2. Value Health. 2011;14(4):429–37. 22. Senn S. The many modes of meta. Drug Inf J. 2000;34(2): 535–49. 23. Colditz GA, Brewer TF, Berkey CS, et al. Efficacy of BCG vaccine in the prevention of tuberculosis: meta-analysis of the published literature. JAMA. 1994;271(9):698–702.
918
24. Borenstein M, Hedges LV, Higgins JPT, Rothstein HR. Meta-regression. In: Borenstein M, Hedges LV, Higgins JPT, Rothstein HR, eds. Introduction to Meta-Analysis. New York: John Wiley & Sons; 2009. p 187–203. 25. Harris R, Bradburn M, Deeks J, Harbord R, Altman D, Sterne J. metan: fixed- and random-effects meta-analysis. Stata J. 2008;8(1): 3–28. 26. Smith TC, Spiegelhalter DJ, Thomas A. Bayesian approaches to random-effects meta-analysis: a comparative study. Stat Med. 1995;14(24):2685–99. 27. Piepho HP, Williams ER, Madden LV. The use of two-way linear mixed models in multitreatment meta-analysis. Biometrics. 2012;68(4):1269–77. 28. Pagliaro L, D’Amico G, So¨rensen TIA, et al. Prevention of first bleeding in cirrhosis: a meta-analysis of randomized trials of nonsurgical treatment. Ann Intern Med. 1992;117(1):59–70. 29. White IR. Multivariate random-effects meta-regression: updates to mvmeta. Stata J. 2011;11(2):255–70. 30. White IR, Barrett JK, Jackson D, Higgins JPT. Consistency and inconsistency in network meta-analysis: model estimation using multivariate meta-regression. Res Synth Methods. 2012;3(2): 111–25. 31. Biondi-Zoccai G, D’Ascenzo F. Choosing between frequentist and Bayesian frameworks and the corresponding statistical package. In: Biondi-Zoccai G, ed. Network Meta-Analysis: Evidence Synthesis with Mixed Treatment Comparison. Hauppauge, NY: Nova Science Publishers; 2014. p 101–13.
MEDICAL DECISION MAKING/OCTOBER 2014
Downloaded from mdm.sagepub.com at Univ. of Tasmania Library on May 21, 2015
