Network meta-analysis is a technique to combine information from all the randomized comparisons among a set of several treatments. The estimate for the relative effectiveness of two treatments uses both the direct head-to-head comparisons of these two treatments and indirect comparisons where each of the two is compared to a third treatment.

Network meta-analysis was the first such technique to allow for the possibility that the indirect information might be biased, and to assess the degree of agreement between direct and indirect comparisons from the data.

The example presented in the original paper has some typographical errors. This page shows a worked example of a network meta-analysis, from a paper on treatments for glaucoma, by van der Valk and others.

The data are in comma-separated text format. The file analysis.R does the complete analysis. The file glaucoma.rda has preprocessed data in a suitable form for fitting the linear mixed models, with a predictor for each drug, coded +/-1 for the two drugs in a pairwise comparison and zero for drugs not involved. The data have also been reweighted, with variances multiplied by 2 for comparisons within a three-armed trial (since each person appears twice and would otherwise be double-counted).

The models can be fitted with the following code

library(nlme)
load("glaucoma.rda")
trough.lme<-lme(mu ~ 0 + betaxolol + bimatoprost + brimonidine + brinzolamide + 
    dorzolamide + latanoprost + placebo + travoprost,
     ~1|trtpair,data=troughdata,
     weights=varConstPower(form=~sigma, fixed=list(power=1)))
peak.lme<-lme(mu ~ 0 + betaxolol + bimatoprost + brimonidine + brinzolamide + 
    dorzolamide + latanoprost + placebo + travoprost,
    ~1|trtpair,data=peakdata,
    weights=varConstPower(form=~sigma, fixed=list(power=1)))
The output from the models is below. The treatment comparisons are the fixed effect estimates. The incoherence is the intercept random effect standard deviation (effectively zero in these data).
> summary(peak.lme)
Linear mixed-effects model fit by REML
 Data: peakdata 
        AIC       BIC   logLik
  -77.98308 -63.32883 49.99154

Random effects:
 Formula: ~1 | trtpair
         (Intercept)   Residual
StdDev: 1.747687e-06 0.05480966

Variance function:
 Structure: Constant plus power of variance covariate
 Formula: ~sigma 
 Parameter estimates:
    const     power 
0.4123657 1.0000000 
Fixed effects: attr(peakdata, "formula") 
                   Value  Std.Error DF   t-value p-value
betaxolol     0.06763815 0.01331770  9  5.078817  0.0007
bimatoprost  -0.08429832 0.01408424  9 -5.985294  0.0002
brimonidine  -0.01490472 0.01741234  9 -0.855986  0.4142
brinzolamide  0.08437913 0.02505425  9  3.367856  0.0083
dorzolamide   0.04280303 0.01409322  9  3.037135  0.0141
latanoprost  -0.06449534 0.01089720  9 -5.918524  0.0002
placebo       0.18447008 0.01791467  9 10.297150  0.0000
travoprost   -0.05928840 0.01363291  9 -4.348918  0.0019
 Correlation: 
             betxll bmtprs brmndn brnzlm drzlmd ltnprs placeb
bimatoprost  0.071                                           
brimonidine  0.028  0.172                                    
brinzolamide 0.153  0.179  0.062                             
dorzolamide  0.323  0.220  0.086  0.473                      
latanoprost  0.088  0.549  0.314  0.199  0.273               
placebo      0.173  0.328  0.107  0.507  0.534  0.341        
travoprost   0.044  0.449  0.131  0.105  0.136  0.417  0.187 

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.06964198 -0.43427485 -0.03559876  0.66911129  2.07421950 

Number of Observations: 36
Number of Groups: 17 
> summary(trough.lme)
Linear mixed-effects model fit by REML
 Data: troughdata 
        AIC       BIC   logLik
  -81.00163 -67.16256 51.50081

Random effects:
 Formula: ~1 | trtpair
         (Intercept)  Residual
StdDev: 1.626263e-08 0.1495982

Variance function:
 Structure: Constant plus power of variance covariate
 Formula: ~sigma 
 Parameter estimates:
       const        power 
3.641353e-10 1.000000e+00 
Fixed effects: attr(troughdata, "formula") 
                   Value   Std.Error DF   t-value p-value
betaxolol     0.05430220 0.010817085  9  5.020040  0.0007
bimatoprost  -0.06170327 0.011152761  9 -5.532556  0.0004
brimonidine   0.07376523 0.015080105  9  4.891559  0.0009
brinzolamide  0.05888197 0.016221087  9  3.629964  0.0055
dorzolamide   0.05503531 0.010456848  9  5.263088  0.0005
latanoprost  -0.02632711 0.009074042  9 -2.901365  0.0176
placebo       0.14977410 0.012245448  9 12.231002  0.0000
travoprost   -0.02748992 0.011463645  9 -2.398009  0.0400
 Correlation: 
             betxll bmtprs brmndn brnzlm drzlmd ltnprs placeb
bimatoprost  0.088                                           
brimonidine  0.016  0.093                                    
brinzolamide 0.168  0.258  0.043                             
dorzolamide  0.305  0.290  0.053  0.550                      
latanoprost  0.091  0.522  0.178  0.242  0.297               
placebo      0.184  0.436  0.069  0.572  0.604  0.388        
travoprost   0.041  0.324  0.067  0.114  0.135  0.376  0.187 

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.0618580 -0.5756033 -0.1307495  0.4804651  1.9255973 

Number of Observations: 34
Number of Groups: 17 
The file analysis.R gives tider summaries corresponding to Table 1 in the paper
>  nice.summary(peak.lme)

 peak.lme 
Random effects:
 Formula: ~1 | trtpair
         (Intercept)   Residual
StdDev: 1.747687e-06 0.05480966

Variance function:
 Structure: Constant plus power of variance covariate
 Formula: ~sigma 
 Parameter estimates:
    const     power 
0.4123657 1.0000000 
              diff   low    up pvalue
betaxolol     0.07  0.04  0.09  0.000
bimatoprost  -0.08 -0.11 -0.06  0.000
brimonidine  -0.01 -0.05  0.02  0.392
brinzolamide  0.08  0.04  0.13  0.001
dorzolamide   0.04  0.02  0.07  0.002
latanoprost  -0.06 -0.09 -0.04  0.000
placebo       0.18  0.15  0.22  0.000
travoprost   -0.06 -0.09 -0.03  0.000
>  nice.summary(trough.lme)

 trough.lme 
Random effects:
 Formula: ~1 | trtpair
         (Intercept)  Residual
StdDev: 1.626263e-08 0.1495982

Variance function:
 Structure: Constant plus power of variance covariate
 Formula: ~sigma 
 Parameter estimates:
       const        power 
3.641353e-10 1.000000e+00 
              diff   low    up pvalue
betaxolol     0.05  0.03  0.08  0.000
bimatoprost  -0.06 -0.08 -0.04  0.000
brimonidine   0.07  0.04  0.10  0.000
brinzolamide  0.06  0.03  0.09  0.000
dorzolamide   0.06  0.03  0.08  0.000
latanoprost  -0.03 -0.04 -0.01  0.004
placebo       0.15  0.13  0.17  0.000
travoprost   -0.03 -0.05 -0.01  0.016
and forestplots (which the journal didn't want)