• Overview
• Searching The Literature
• Fixed Effects Model
• Random Effects Model
• Evaluating Heterogeneity
• Meta-Regression
• Publication Bias
• Comparing R Packages For Standard Meta-Analysis
• Some Advanced Topics

Overview

The analysis of analyses.

-- Gene V. Glass

Primary, secondary and meta-analysis of research,

Educational Researcher, 1976.

More formally...a meta-analysis is the synthesis of:

• \(K\) compatible effects (\(Y_i\))

• (Preferably, but not necessarily, from randomized controlled trials)

Giving greater weight to studies with:

• Less variance (\(V_i\)), and

• More precision (\(W_i = 1/V_i\))

An "effect" could be almost any aggregate statistic of interest:

• Mean, Mean difference, Mean change

• Risk ratio, Odds ratio, Risk difference

• Incidence rate, Prevalence, Proportion

• Correlation

Source: Nissen SE, Wolski K. N Engl J Med 2007;356:2457-71.

Conducting Meta-Analyses in R

• Official home page: http://www.R-project.org

• Introduction under "What is R?"

• Download base system at: http://cran.r-project.org

• Extend with user-contributed packages -> install.packages

• Find further introductory material with help.start

• R is a free, open-source, & powerful statistical environment

• Run on Windows, Mac OS, and Linux platforms

• Has 20+ meta-analytic packages on CRAN

• Tools for meta-regression, Bayesian meta-analysis, multivariate meta-analyses, etc.

• Easy (in most cases) to customize and extend these tools

• A meta-analysis starts with a systematic review.

• A systematic review is a scientific summary of all available evidence on a specific research question.

• An exhaustive search of the literature will require more than R.

• Note: If available studies are too few or too different a meta-analysis may not be appropriate.

Searching The Literature With R

• Not many packages for helping with early stages of a systematic review.

• But, I created the RISmed package to import metadata from NCBI databases into R.

• Using this package, one can search, store, and easily mine metadata on PubMed articles.

• RISmed tools are not comprehensive enough to complete a systematic review but may be a helpful aid.

1. Create EUtilsSummary object for specified query.

2. Retrieve matching records with EUtilsGet.

Syntax

EUtilsSummary( [query], [db], [search.limits])

• query: String query as given on PubMed site

• db: String name of NCBI database

• search.limits: Additional arguments to restrict search

library(RISmed)  # Load Package

The following code performs a PubMed query of all BMJ articles with "rofecoxib" in the title.

fit <- EUtilsSummary("rofecoxib[ti]+British Medical Journal[jo]", db = "pubmed")

QueryTranslation(fit)  # Extract the translated query

## [1] "rofecoxib[ti] AND (\"Br Med J\"[Journal] OR \"Br Med J (Clin Res Ed)
\"[Journal] OR \"BMJ\"[Journal])"

QueryCount(fit)  # Extract the number of matched records

## [1] 16

Now we can extract the metadata for the queried records.

fetch <- EUtilsGet(fit)
fetch  # Medline Object

## PubMed query:  rofecoxib[ti] AND ("Br Med J"[Journal] OR "
## Br Med J (Clin Res Ed)"[Journal] OR "BMJ"[Journal]) 
## 
## Records:  16

getSlots("Medline")  # Available methods

##                Query                 PMID                 Year 
##          "character"          "character"            "numeric" 
##                Month                  Day               Author 
##            "numeric"            "numeric"               "list" 
##                 ISSN                Title         ArticleTitle 
##          "character"          "character"          "character" 
##          ELocationID         AbstractText          Affiliation 
##          "character"          "character"          "character" 
##             Language      PublicationType            MedlineTA 
##          "character"          "character"          "character" 
##          NlmUniqueID          ISSNLinking                 Hour 
##          "character"          "character"            "numeric" 
##               Minute    PublicationStatus            ArticleId 
##            "numeric"          "character"          "character" 
##               Volume                Issue      ISOAbbreviation 
##          "character"          "character"          "character" 
##           MedlinePgn CopyrightInformation              Country 
##          "character"          "character"          "character" 
##              GrantID              Acronym               Agency 
##          "character"          "character"          "character" 
##       RegistryNumber            RefSource       CollectiveName 
##          "character"          "character"          "character" 
##                 Mesh 
##               "list"

ArticleTitle(fetch)[1:5]

## [1] "Merck pays $1bn penalty in relation to promotion of rofecoxib."      
## [2] "Merck to pay $58m in settlement over rofecoxib advertising."         
## [3] "94% of patients suing Merck over rofecoxib agree to company's offer."
## [4] "Merck to pay $5bn in rofecoxib claims."                              
## [5] "Merck appeals rofecoxib verdict."

Author(fetch)[[1]]

##   LastName       ForeName Initials order
## 1    Tanne Janice Hopkins       JH     1

Year(fetch)

##  [1] 2011 2008 2008 2007 2007 2006 2005 2005 2004 2004 2004 2004 2004 2003
## [15] 2002 2001

Using the "rofecoxib" Medline object,

1. Determine the first year a matching article appeared.

2. What was the title of this article?

3. Do some authors have multiple matching records?

4. If so, which authors?

min(Year(fetch))  # Earliest year

## [1] 2001

ArticleTitle(fetch)[Year(fetch) == 2001]  # Title of earliest record(s)

## [1] "FDA warns Merck over its promotion of rofecoxib."

AuthorList <- Author(fetch)  # Extract list of authors
LastFirst <- sapply(AuthorList, function(x) paste(x$LastName, x$ForeName))
sort(table(unlist(LastFirst)), dec = TRUE)[1:3]  # Tabulate & Sort

## Tanne Janice Hopkins        Charatan Fred      Abenhaim Lucien 
##                    4                    3                    1

Basic Meta-Analysis In R

In no particular order...

• meta (Author: Guido Schwarzer)

• metafor (Author: Wolfgang Viechtbauer)

• rmeta (Author: Thomas Lumley)

1. BCG vaccine trials (from metafor)

2. Amlodipine angina treatment trials (from meta)

• Overview: 13 vaccine trials of Bacillus Calmette–Guérin (BCG) vaccine vs. no vaccine
• Treatment goal: Prevention of tuberculosis
• Primary endpoint: Tuberculosis infection
• Possible explanatory variables:
  • latitude of study region
  • treatment allocation method
  • year published

A version of the BCG dataset is provided by package metafor

library(metafor)  # Load package
data(dat.bcg)  # BCG meta-analytic dataset
str(dat.bcg)  # Describe meta-analysis structure

## 'data.frame':    13 obs. of  9 variables:
##  $ trial : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ author: chr  "Aronson" "Ferguson & Simes" "Rosenthal et al" "Hart & Sutherland" ...
##  $ year  : int  1948 1949 1960 1977 1973 1953 1973 1980 1968 1961 ...
##  $ tpos  : int  4 6 3 62 33 180 8 505 29 17 ...
##  $ tneg  : int  119 300 228 13536 5036 1361 2537 87886 7470 1699 ...
##  $ cpos  : int  11 29 11 248 47 372 10 499 45 65 ...
##  $ cneg  : int  128 274 209 12619 5761 1079 619 87892 7232 1600 ...
##  $ ablat : int  44 55 42 52 13 44 19 13 27 42 ...
##  $ alloc : chr  "random" "random" "random" "random" ...

• Overview: 8 randomized controlled trials (RCTs) of amlodipine vs. placebo
• Treatment goal: Reduce harms of angina (chest pain)
• Primary endpoint: Work capacity (ratio of exercise time after to before intervention)

A version of the amlodipine dataset is provided by package meta

library(meta)  # Load package
data(amlodipine)  # amlodipine meta-analytic dataset
str(amlodipine)  # Describe meta-analysis structure

## 'data.frame':    8 obs. of  7 variables:
##  $ study    : Factor w/ 8 levels "Protocol 154",..: 1 2 3 4 5 6 7 8
##  $ n.amlo   : int  46 30 75 12 32 31 27 46
##  $ mean.amlo: num  0.232 0.281 0.189 0.093 0.162 ...
##  $ var.amlo : num  0.2254 0.1441 0.1981 0.1389 0.0961 ...
##  $ n.plac   : int  48 26 72 12 34 31 27 47
##  $ mean.plac: num  -0.0027 0.027 0.0443 0.2277 0.0056 ...
##  $ var.plac : num  0.0007 0.1139 0.4972 0.0488 0.0955 ...

• An effect size could be almost any summary statistic (e.g. a mean, a difference in proportions, an adjusted odds ratio, etc.)

• Conventional meta-analytic models assume normality of ESs.

• Because of the CLT, this will holds for most ESs given large enough samples.

• To normalize ESs, a log-transform is common.

Effect Size

\[ LOR = \log(\frac{a * d}{ b * c}) \]

Variance

\[ V = 1/a+1/b+1/c+1/d \]

Y <- with(dat.bcg, log(tpos * cneg/(tneg * cpos)))
V <- with(dat.bcg, 1/tpos + 1/cneg + 1/tneg + 1/cpos)
cbind(Y, V)

##              Y        V
##  [1,] -0.93869 0.357125
##  [2,] -1.66619 0.208132
##  [3,] -1.38629 0.433413
##  [4,] -1.45644 0.020314
##  [5,] -0.21914 0.051952
##  [6,] -0.95812 0.009905
##  [7,] -1.63378 0.227010
##  [8,]  0.01202 0.004007
##  [9,] -0.47175 0.056977
## [10,] -1.40121 0.075422
## [11,] -0.34085 0.012525
## [12,]  0.44663 0.534162
## [13,] -0.01734 0.071635

• escalc does the work of calculating ESs.

• Give the necessary data components (i.e. sample size, events in each treatment group, etc.).

• Indicate the ES you want with measure.

• Many, many, types of single group and between-group ESs.

Syntax

ES <- escalc(endpoints, variances, measure, data, ...)

• endpoints: arguments or formula containing endpoint values

• variances: arguments containing endpoint variances

• measure: character value indicating type of ES

• data: data frame containing named variables

ES <- escalc(ai = tpos, bi = tneg, ci = cpos, di = cneg, 
    data = dat.bcg, 
    measure = "OR")

cbind(ES$yi, ES$vi)

##           [,1]     [,2]
##  [1,] -0.93869 0.357125
##  [2,] -1.66619 0.208132
##  [3,] -1.38629 0.433413
##  [4,] -1.45644 0.020314
##  [5,] -0.21914 0.051952
##  [6,] -0.95812 0.009905
##  [7,] -1.63378 0.227010
##  [8,]  0.01202 0.004007
##  [9,] -0.47175 0.056977
## [10,] -1.40121 0.075422
## [11,] -0.34085 0.012525
## [12,]  0.44663 0.534162
## [13,] -0.01734 0.071635

• What if my data is in a "long" format?

• That is, what if I have multiple rows per study, corresponding to difference treatment groups?

• In that case, you may prefer specifying the variables for the ES calculation using a formula.

Syntax

escalc(formula = outcome ~ group | study, data = data, weights = n)

Note: The exact syntax will vary a bit depending on the ES type.

library(reshape2)  # Load package for data reshaping

bcg.long <- melt(dat.bcg[, c("trial", "tpos", "tneg", "cpos", "cneg")], id = "trial")
bcg.long$pos <- ifelse(bcg.long$var == "tpos" | bcg.long$var == "cpos", 1, 0)
bcg.long$group <- ifelse(bcg.long$var == "tpos" | bcg.long$var == "tneg", 1, 0)
head(bcg.long)

##   trial variable value pos group
## 1     1     tpos     4   1     1
## 2     2     tpos     6   1     1
## 3     3     tpos     3   1     1
## 4     4     tpos    62   1     1
## 5     5     tpos    33   1     1
## 6     6     tpos   180   1     1

escalc(factor(pos) ~ factor(group) | factor(trial), 
    weights = value, 
    data = bcg.long, 
    measure = "OR")

##         yi     vi
## 1  -0.9387 0.3571
## 2  -1.6662 0.2081
## 3  -1.3863 0.4334
## 4  -1.4564 0.0203
## 5  -0.2191 0.0520
## 6  -0.9581 0.0099
## 7  -1.6338 0.2270
## 8   0.0120 0.0040
## 9  -0.4717 0.0570
## 10 -1.4012 0.0754
## 11 -0.3408 0.0125
## 12  0.4466 0.5342
## 13 -0.0173 0.0716

• rma is the main modeling function of metafor.

• rma is also a wrapper for escalc, and will compute ESs before modeling (if you like).

• You usually won't work with escalc directly.

• But if you want the ESs (yi) and variances (vi) without modeling, use escalc.

• By default, escalc appends the yi and vi to the dataset.

• To return only yi and vi set append=TRUE.

Summarizing Effects

Normal Assumption \[ Y_i \sim N(\theta, V_i) \]

Summary Effect Size is a Weighted Average \[ \hat{\theta} = \sum_i Y_i W_i / \sum_i W_i,\; \mbox{Var}(\hat{\theta}) = 1/\sum_i W_i \]

Each Study's Contribution \[ \lambda_i = W_i/\sum_i W_i \]

• Fixed

  -> Same mean ES, zero between-study variance

• Random

  -> Different mean ES, between-study variance

• Mixed

  -> Study-level regression for mean ES

• Same mean ES, known variance

\[ Y_i = \theta+e_i, \]

\[ e_i \sim N(0, V_i). \]

• Different mean ES, between-study variance

\[ Y_i = \theta+\theta_i+e_i, \]

\[ \theta_i \sim N(0, \tau^2), \]

\[ e_i \sim N(0, V_i). \]

• Study-level regression for mean ES

\[ Y_i = \boldsymbol{\beta}'\mathbf{x}_i+\theta_i+e_i, \]

\[ \theta_i \sim N(0, \tau^2), \]

\[ e_i \sim N(0, V_i). \]

\(\mathbf{x}_i\) = Study-level covariates

• The FE model is a description of the \(K\) studies.

• The RE model regards the \(K\) studies as a sample of a larger universe of studies.

• The RE model can be used to infer what would likely happen if a new study were performed, the FE model cannot.

• Common practice is to report both fixed and random effects model results.

Summary ES has more uncertainty because of between-study variance.

• All model summaries are made with the rma function.

• rma stands for random effects meta-analysis

• The default method is a REML RE model, but the FE model can also be fit.

REML = Restricted Maximum Likelihood

Syntax

rma(yi, vi, method, ...)

• yi effect size

• vi variances

• method type of model approach

• "FE" = Fixed Effects

• "DL" = DerSimonian-Laird

• "HE" = Hedges estimator

• "ML" = Maximum Likelihood

• "REML" = Restricted ML

Fitting The Fixed Effects Model

result.or <- rma(yi = Y, vi = V, method = "FE")  # Log Odds Ratio

summary(result.or)

## Fixed-Effects Model (k = 13)
## 
##   logLik  deviance       AIC       BIC  
## -76.0290  163.1649  154.0580  154.6229  
## 
## Test for Heterogeneity: 
## Q(df = 12) = 163.1649, p-val < .0001
## 
## Model Results:
## 
## estimate       se     zval     pval    ci.lb    ci.ub          
##  -0.4361   0.0423 -10.3190   <.0001  -0.5190  -0.3533      *** 

args(rma)

## function (yi, vi, sei, weights, ai, bi, ci, di, n1i, n2i, x1i, 
##     x2i, t1i, t2i, m1i, m2i, sd1i, sd2i, xi, mi, ri, ti, sdi, 
##     ni, mods, measure = "GEN", intercept = TRUE, data, slab, 
##     subset, add = 1/2, to = "only0", drop00 = FALSE, vtype = "LS", 
##     method = "REML", weighted = TRUE, knha = FALSE, level = 95, 
##     digits = 4, btt, tau2, verbose = FALSE, control) 
## NULL

Have rma calculate the ESs, if you haven't done it yourself.

result.md <- rma(m1 = mean.amlo, m2 = mean.plac, 
    sd1 = sqrt(var.amlo), sd2 = sqrt(var.plac), 
    n1 = n.amlo, n2 = n.plac, 
    method = "FE", measure = "MD", 
    data = amlodipine)

-> The function rma returns an object of the class rma.

-> This object behaves like a list.

-> You can use the function names to see available elements.

names(result.md)  # Components of rma

names(result.md)

##  [1] "b"         "se"        "zval"      "pval"      "ci.lb"    
##  [6] "ci.ub"     "vb"        "tau2"      "se.tau2"   "k"        
## [11] "k.f"       "k.eff"     "p"         "p.eff"     "parms"    
## [16] "m"         "QE"        "QEp"       "QM"        "QMp"      
## [21] "I2"        "H2"        "int.only"  "int.incl"  "allvipos" 
## [26] "yi"        "vi"        "X"         "yi.f"      "vi.f"     
## [31] "X.f"       "ai.f"      "bi.f"      "ci.f"      "di.f"     
## [36] "x1i.f"     "x2i.f"     "t1i.f"     "t2i.f"     "ni"       
## [41] "ni.f"      "ids"       "not.na"    "slab"      "slab.null"
## [46] "measure"   "method"    "weighted"  "knha"      "robust"   
## [51] "s2w"       "btt"       "intercept" "digits"    "level"    
## [56] "control"   "add"       "to"        "drop00"    "fit.stats"

For the mean difference in the amlodipine trial determine:

1. The summary effect

2. The 95% confidence interval

result.md$b

##           [,1]
## intrcpt 0.1619

c(result.md$ci.lb, result.md$ci.ub)

## [1] 0.0986 0.2252

1. Determine the percentage each study contributed to the overall effect size summary.

2. Which study contributes the most? How much?

3. Use a barplot to show the percentages graphically.

contributions <- 1/result.md$vi/sum(1/result.md$vi) * 100

cbind(contributions)

##      contributions
## [1,]        21.219
## [2,]        11.355
## [3,]        10.923
## [4,]         6.667
## [5,]        17.943
## [6,]        10.848
## [7,]         1.661
## [8,]        19.385

max(contributions)

## [1] 21.22

amlodipine$study[which(contributions == max(contributions))]

## [1] Protocol 154
## 8 Levels: Protocol 154 Protocol 156 Protocol 157 ... Protocol 306

contributions <- 1/result.md$vi/sum(1/result.md$vi) * 100

par(mar = c(5, 10, 5, 5))

barplot(contributions, names = amlodipine$study, 
	xlim = c(0, 50), las = 2, horiz = T, 
	col = "royalblue")

summary(result.md)

## Fixed-Effects Model (k = 8)
## 
##   logLik  deviance       AIC       BIC  
##   4.7834   12.3311   -7.5669   -7.4874  
## 
## Test for Heterogeneity: 
## Q(df = 7) = 12.3311, p-val = 0.0902
## 
## Model Results:
## 
## estimate       se     zval     pval    ci.lb    ci.ub          
##   0.1619   0.0323   5.0134   <.0001   0.0986   0.2252      *** 

coef(result.md)

## intrcpt 
##  0.1619

confint(result.md)  # Heterogeneity measures do not apply for FE model

## 
##        estimate  ci.lb   ci.ub
## tau^2        NA 0.0000  0.1667
## tau          NA 0.0000  0.4082
## I^2(%)       NA 0.0000 95.0658
## H^2          NA 1.0000 20.2667

Fitting The Random Effects Model

• Suppose between-study variance (\(\tau^2\)) is non-zero.

• Methods differ on how they estimate \(\tau^2\).

• Many iterative and non-iterative approaches to estimating \(\tau^2\) have been proposed.

The rma function offers the following estimators:

The rma function offers the following estimators:

No method is universally superior, but Viechtbauer's simulation study (2002) suggests REML has the most recommendable properties.

result.md <- rma(m1 = mean.amlo, m2 = mean.plac, 
    sd1 = sqrt(var.amlo), sd2 = sqrt(var.plac), 
    n1 = n.amlo, n2 = n.plac, 
    method = "REML", measure = "MD", 
    data = amlodipine)

summary(result.md)

## Random-Effects Model (k = 8; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC  
##   3.3094   -6.6188   -2.6188   -2.7270  
## 
## tau^2 (estimated amount of total heterogeneity): 0.0001 (SE = 0.0042)
## tau (square root of estimated tau^2 value):      0.0116
## I^2 (total heterogeneity / total variability):   1.54%
## H^2 (total variability / sampling variability):  1.02
## 
## Test for Heterogeneity: 
## Q(df = 7) = 12.3311, p-val = 0.0902
## 
## Model Results:
## 
## estimate       se     zval     pval    ci.lb    ci.ub          
##   0.1617   0.0326   4.9584   <.0001   0.0978   0.2257      *** 

result.md$tau2

## [1] 0.0001353

result.md$se.tau2

## [1] 0.004239

result.md$tau2 + 1.96 * c(-1, 1) * result.md$se.tau2  #95% CI

## [1] -0.008173  0.008444

estimators <- c("DL", "REML", "HE", "HS", "SJ", "ML", "EB")

taus <- sapply(estimators, function(method) {
    rma(m1 = mean.amlo, m2 = mean.plac, 
    	sd1 = sqrt(var.amlo), sd2 = sqrt(var.plac), 
        n1 = n.amlo, n2 = n.plac, 
        method = method, measure = "MD", 
        data = amlodipine)$tau2
})

plot(y = taus, x = 1:length(taus), 
    type = "h", pch = 19, 
    axes = FALSE, xlab = "Estimators")
    
axis(2, las = 1)

axis(1, at = 1:length(taus), lab = estimators)

DerSimonian-Laird

Method of moments estimator; Most popular approach

\[ \hat{\tau}^2 = \mbox{max} \lbrace 0, \frac{Q - (K-1)}{\sum_i W_i - \sum_i W_i^2/\sum_i W_i } \rbrace \]

\[ Q = \sum_i W_i (Y_i - \bar{Y})^2 \]

$ \bar{Y}_W = \sum W_i Y_i / \sum W_i$

REML

Best properties, in general

\[ \hat{\tau}^2 = \frac{\sum_i \tilde{W}_i^2 \left[ \frac{K}{K-1} (Y_i - \hat{\theta})^2-V_i \right]}{\sum_i \tilde{W}_i^2} \]

K = Number of trials

\(\tilde{W} = (V_i + \hat{\tau}^2)^{-1}\)

\(\hat{\theta}\) = Effect size

Maximum-Likelihood

\[ \hat{\tau}^2 = \frac{\sum_i W_i^2 \left[ (Y_i - \hat{\theta})^2-V_i \right]}{\sum_i W_i^2} \]

Hedges

\[ \hat{\tau}^2 = \frac{\sum_i (Y_i - \bar{Y})^2}{(K-1)}-\frac{ \sum_i V_i }{K} \]

Sidik-Jonkman

\[ \hat{\tau}_0^2 = k^{-1}\sum_i (Y_i - \bar{Y})^2 \]

\[ \hat{\theta} = (\sum_i \tilde{W}_i)^{-1} \sum_i \tilde{W}_i Y_i \]

\[ \hat{\tau}_2 = \hat{\tau}_0^2 /(K-1) \sum_i \tilde{W}_i (Y_i - \hat{\theta})^2 \]

10-minute Break

Evaluating Heterogeneity

\[ Q = \sum_i W_i (Y_i - \hat{\theta})^2 \]

• Q is the weighted deviations about the summary effect size.

• Larger values of Q reflect greater between-study heterogeneity.

• When \(\tau^2 = 0\), \(Q \sim \chi^2 (K-1)\), which leads to a chi-squared test for heterogeneity.

MD <- with(amlodipine, mean.amlo - mean.plac)

W <- 1/with(amlodipine, var.amlo/n.amlo + var.plac/n.plac)

Q <- sum(W * (MD - sum(W * MD)/sum(W))^2)

Q

## [1] 12.33

df <- length(MD) - 1
pchisq(Q, df = df, lower = FALSE)  # HOW LIKELY UNDER NULL?

## [1] 0.09018

curve(1 - pchisq(x, df = df), 0, 20)
abline(v = Q, col = "red", lwd = 2)

result.md$QE

## [1] 12.33

result.md$QEp

## [1] 0.09018

1. Obtain the Q-test for the meta-analysis of log odds ratios with the BCG Vaccine Trials.

2. What does the test suggest about between-study heterogeneity?

result.or <- rma(yi = Y, vi = V, method = "DL")  # DerSimonian-Laird

result.or$QE

## [1] 163.2

result.or$QEp

## [1] 1.189e-28

• The chi-squared approximation is valid when study sample sizes are large.

• Type I error is generally accurate if normal distribution assumption and sample sizes are not too small.

• Q-test has low power (<0.80) when the number of studies and/or sample sizes is small.

Bottom Line: If there are few trials in the meta-analysis (as is usually the case), the Q-test is likely underpowered for detecting true heterogeneity.

• \(\tau^2\)

• Higgins' \(I^2\)

• \(H^2\), \(H\) Index

• Intra-class correlation (ICC)

\[ I^2 = (Q-df)/Q \times 100 \]

Interpretation = Percentage of "unexplained" variance

df = Degrees of Freedom

For random-effects meta-analysis, df = \(K-1\)

Judging the severity of measured heterogeneity is subjective, however Higgins suggests these rules of thumb:

• 0% to 30% \(\to\) Low

• 30% to 60% \(\to\) Moderate

• 50% to 90% \(\to\) Substantial

• 75% to 100% \(\to\) Considerable

(Q - df)/Q * 100

## [1] 43.23

# From rma object
I2 <- with(result.md, (QE - (k - 1))/QE * 100)
I2

## [1] 43.23

Is the ratio of \(Q\) to the Q-test's degrees of freedom,

\[ H^2 = \frac{Q}{df}, \]

\[ 1/H^2 = 1- \frac{I^2}{100}. \]

\(H\) index is the \(\sqrt{H^2}\).

\(H>1\) suggests there is unexplained heterogeneity.

Q/df

## [1] 1.762

1/(1 - I2/100)

## [1] 1.762

After fitting RMA and getting measure of \(\tau^2\), we can compute the intra-class correlation (ICC)

\[ ICC = \frac{\tau^2}{\tau^2 + S^2} \]

\[ S^2 = \frac{\sum W_i (K-1) }{(\sum W_i)^2 - \sum W_i^2} \]

S2 <- sum(W * (length(W) - 1))/(sum(W)^2 - sum(W^2))

result.md$tau2/(result.md$tau2 + S2)

## [1] 0.0154

result.md$I2

## [1] 1.54

• What happened in the previous example?

• We saw that \(I^2 = ICC \times 100\)

• This is because metafor uses the more general definitions of \(I^2\) and \(H^2\), which are based on \(\tau^2\).

• To get the conventional estimates, which do not depend on \(\tau^2\), use method DL.

\[ I^2 = ICC \times 100 \]

\[ H^2 = (\tau^2+\sigma^2)/\sigma^2 \]

where, \(\sigma^2\) is the weighted numerator of the DL \(\tau^2\) estimator

\[ \sigma^2 = \left[ (K-1)(\sum W_i - \sum W_i^2/ \sum W_i )\right]^{-1} \]

result.md$tau2/(result.md$tau2 + S2) * 100

## [1] 1.54

result.md$I2

## [1] 1.54

sigma2 <- (length(Y) - 1) * (sum(W) - sum(W^2)/sum(W))^(-1)

result.md$tau2/sigma2 + 1  #H2

## [1] 1.009

result.md$H2

## [1] 1.016

result.md <- rma(m1 = mean.amlo, m2 = mean.plac, 
    sd1 = sqrt(var.amlo), sd2 = sqrt(var.plac), 
    n1 = n.amlo, n2 = n.plac, 
    method = "DL", measure = "MD", data = amlodipine)

result.md$I2

## [1] 43.23

result.md$H2

## [1] 1.762

• A Q-profile method for an exact confidence interval for \(\tau^2\) is provided with the confint method of rma objects.

• The CI for \(\tau^2\) is used to derive CIs for the remaining heterogeneity indices, which are all monotonic transformations of \(\tau^2\).

confint(result.md)

## 
##        estimate  ci.lb   ci.ub
## tau^2    0.0066 0.0000  0.1667
## tau      0.0812 0.0000  0.4082
## I^2(%)  43.2328 0.0000 95.0658
## H^2      1.7616 1.0000 20.2667

1. Use the confint method to obtain a 95% CI for the ICC of the mean difference DL meta-analysis.

names(confint(result.md))

## [1] "random"   "digits"   "tau2.min"

I2CI <- confint(result.md)$random[3, ]

I2CI/100

## estimate    ci.lb    ci.ub 
##   0.4323   0.0000   0.9507

Visualizing Heterogeneity

Seeing the forest through the trees...

• Is a plot of effect sizes and their precisions

• Is the most common way to report the results of a meta-analysis

• Can help identify patterns across effects

• Can help spot large variation in effects or possible outliers

forest(result.md)  # DEFAULT PLOT

args(forest.rma)

## function (x, annotate = TRUE, addfit = TRUE, addcred = FALSE, 
##     showweight = FALSE, xlim, alim, ylim, at, steps = 5, level = x$level, 
##     digits = 2, refline = 0, xlab, slab, mlab, ilab, ilab.xpos, 
##     ilab.pos, order, transf = FALSE, atransf = FALSE, targs, 
##     rows, efac = 1, pch = 15, psize, col = "darkgray", border = "darkgray", 
##     cex, cex.lab, cex.axis, ...) 
## NULL

Some typical modifications:

• order: Sort by "obs", "fit", "prec", etc.
• slab: Change study labels
• ilab: Add study information
• transf: Apply function to effects
• psize: Symbol sizes

In the following, we modify the study labels and add the (fake) year of publication.

study.names <- paste("Study", letters[1:8])

study.year <- 2000 + sample(0:9, 8, replace = T)

forest(result.md, order = "obs", 
	slab = study.names, 
	ilab = study.year, 
	ilab.xpos = result.md$b - 1, 
	refline = result.md$b)

In the following, we plot the ORs of the BCG trials and order the studies by precision.

forest(result.or, order = "prec", transf = exp, refline = 1)

1. Modify the previous plot by adding the sample size and year of the studies.

dat.bcg$n <- with(dat.bcg, tpos + tneg + cpos + cneg)

forest(result.or, order = "prec", 
	ilab = dat.bcg[, c("n", "year")], 
	ilab.xpos = exp(result.or$b) - c(4, 2), 
	transf = exp, refline = 1)

• You can change the min and max of the drawn region with the argument alim.

• This must include all effects.

• CIs will be clipped if outside the restricted area.

• An arrow will indicate clipped CIs.

forest(result.or, order = "prec", 
    ilab = dat.bcg[, c("n", "year")], 
    ilab.xpos = exp(result.or$b) - c(4, 2), 
    transf = exp, refline = 1, 
    alim = c(0, 4))

Sensitivity Analyses

• A single outlying trial could be the source of substantial heterogeneity.

• To identify suspicious cases, a leave-one-out method can be used whereby we rerun the meta-analysis, iteratively removing studies.

• In the metafor package this is accomplished with the leave1out function.

leave1out(result.md)

##   estimate     se   zval   pval  ci.lb  ci.ub       Q     Qp   tau2
## 1   0.1434 0.0516 2.7759 0.0055 0.0421 0.2446 10.9770 0.0891 0.0080
## 2   0.1449 0.0497 2.9156 0.0036 0.0475 0.2423 11.2868 0.0799 0.0076
## 3   0.1610 0.0519 3.1045 0.0019 0.0593 0.2626 12.2979 0.0556 0.0090
## 4   0.1833 0.0345 5.3173 0.0000 0.1158 0.2509  6.3053 0.3899 0.0004
## 5   0.1595 0.0541 2.9494 0.0032 0.0535 0.2655 12.3252 0.0551 0.0099
## 6   0.1689 0.0505 3.3429 0.0008 0.0699 0.2679 11.7178 0.0686 0.0082
## 7   0.1481 0.0387 3.8295 0.0001 0.0723 0.2239  8.2003 0.2238 0.0028
## 8   0.1623 0.0543 2.9908 0.0028 0.0559 0.2687 12.2425 0.0568 0.0099
##        I2     H2
## 1 45.3404 1.8295
## 2 46.8405 1.8811
## 3 51.2112 2.0496
## 4  4.8426 1.0509
## 5 51.3192 2.0542
## 6 48.7959 1.9530
## 7 26.8316 1.3667
## 8 50.9905 2.0404

1. Which trial contributes the most to the BCG OR meta-analysis?

2. Do any of the trials reduce \(I^2\) to \(<30\)%?

3. Does the removal of any trial change the main conclusion about the efficacy of BCG?

cases <- leave1out(result.or)

which(cases$I2 == min(cases$I2))

## [1] 8

sum(cases$I2 < 30)  # Number with low heterogeneity

## [1] 0

cbind(exp(cases$estimate), cases$pval < 0.05)

##         [,1] [,2]
##  [1,] 0.4784    1
##  [2,] 0.5047    1
##  [3,] 0.4885    1
##  [4,] 0.5173    1
##  [5,] 0.4493    1
##  [6,] 0.4835    1
##  [7,] 0.5025    1
##  [8,] 0.4379    1
##  [9,] 0.4605    1
## [10,] 0.5033    1
## [11,] 0.4510    1
## [12,] 0.4499    1
## [13,] 0.4424    1

Explaining Heterogeneity: Meta-Regression

• Specify study covariates through the mods argument

• The mods argument takes a matrix of \(p\) covariates

result.ormr <- rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, 
    data = dat.bcg, 
    mods = dat.bcg[, "ablat"], 
    measure = "OR", 
    method = "DL")

summary(result.ormr)

## Mixed-Effects Model (k = 13; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0480 (SE = 0.0451)
## tau (square root of estimated tau^2 value):             0.2191
## I^2 (residual heterogeneity / unaccounted variability): 56.17%
## H^2 (unaccounted variability / sampling variability):   2.28
## 
## Test for Residual Heterogeneity: 
## QE(df = 11) = 25.0954, p-val = 0.0088
## 
## Test of Moderators (coefficient(s) 2): 
## QM(df = 1) = 26.1628, p-val < .0001
## 
## Model Results:
## 
##          estimate      se     zval    pval    ci.lb    ci.ub     
## intrcpt    0.3030  0.2109   1.4370  0.1507  -0.1103   0.7163     
## mods      -0.0316  0.0062  -5.1150  <.0001  -0.0437  -0.0195  ***

exp(c(result.or$b, result.ormr$b[1]))

## [1] 0.4736 1.3540

c(result.or$I2, result.ormr$I2)

## [1] 92.65 56.17

• What happened?

• The effect of treatment changed direction.

• Remember: This is a linear not logistic regression.

• As fit, the intercept (treatment log-odds) corresponds to a study conducted in a region with latitude = 0.

1. Determine to what extent the study design (alloc) explains the remaining heterogeneity in the BCG vaccine trials.

2. Center latitude on the median, so that the intercept corresponds to the log-odds effect of BCG at the median latitude.

3. What is the percentage change in \(I^2\) as compared to the RE model?

dat.bcg$random <- ifelse(dat.bcg$alloc == "random", 1, 0)
dat.bcg$cablat <- dat.bcg$ablat - median(dat.bcg$ablat)

result.ormr <- rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, 
	data = dat.bcg, 
	mods = dat.bcg[, c("ablat", "random")], 
	measure = "OR", 
	method = "DL")

summary(result.ormr)

## Mixed-Effects Model (k = 13; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0732 (SE = 0.0677)
## tau (square root of estimated tau^2 value):             0.2706
## I^2 (residual heterogeneity / unaccounted variability): 60.10%
## H^2 (unaccounted variability / sampling variability):   2.51
## 
## Test for Residual Heterogeneity: 
## QE(df = 10) = 25.0624, p-val = 0.0052
## 
## Test of Moderators (coefficient(s) 2,3): 
## QM(df = 2) = 20.0425, p-val < .0001
## 
## Model Results:
## 
##          estimate      se     zval    pval    ci.lb    ci.ub     
## intrcpt    0.3643  0.2596   1.4037  0.1604  -0.1444   0.8731     
## ablat     -0.0307  0.0072  -4.2829  <.0001  -0.0447  -0.0166  ***
## random    -0.2029  0.2124  -0.9551  0.3395  -0.6191   0.2134     

c(result.ormr$I2, result.or$I2)

## [1] 60.10 92.65

(result.or$I2 - result.ormr$I2)/result.or$I2 * 100

## [1] 35.13

Publication Bias

• It is possible that studies showing a significant intervention effect are more often published than studies with null results.

• When a meta-analysis is based only on studies reported in the literature, null studies relegated to the file-drawer could bias the summary intervention effect in the direction of efficacy.

• A funnel plot is a scatter plot of the intervention effect estimates against a measure of study precision.

• Asymmetry (gaps) in the funnel may be indicative of publication bias.

• Some authors argue that judging asymmetry is too subjective to be useful.

• Spurious asymmetry can result from heterogeneity or when ESs are correlated with precision.

metafor method for generating funnel plots from rma objects.

funnel(result.md)

Use addtau2=TRUE to add between-study error.

• Judging asymmetry in the funnel plot can be difficult.

• So you will usually want to consider some additional ways of assessing the threat of publication bias.

• Sensitivity Analyses:

  • Trim-and-Fill
  • Fail Safe N

• The trim and fill method estimates the number of missing NULL studies from the meta-analysis.

• The method trimfill of the metafor package augments the observed data and returns the fitted rma object with the missing studies included.

• These points can be added to the funnel plot.

result.rd <- rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, 
    data = dat.bcg, 
    measure = "RD", 
    method = "DL")  # Risk Differences

trimfill(result.rd)  # Only applicable for FE or RE objects

## Estimated number of missing studies on the right side: 4 
## 
## Random-Effects Model (k = 17; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of total heterogeneity): 0.0000 (SE = 0.0000)
## tau (square root of estimated tau^2 value):      0.0051
## I^2 (total heterogeneity / total variability):   95.83%
## H^2 (total variability / sampling variability):  23.98
## 
## Test for Heterogeneity: 
## Q(df = 16) = 383.6062, p-val < .0001
## 
## Model Results:
## 
## estimate       se     zval     pval    ci.lb    ci.ub          
##  -0.0049   0.0018  -2.7858   0.0053  -0.0084  -0.0015       ** 

funnel(trimfill(result.rd))

• Rosenthal method (sometimes called a ‘file drawer analysis’)

• Is the number of NULL studies that have to be added to reduce the significance of the meta-analysis to \(\alpha\) (usually 0.05)

value <- fsn(y = result.md$yi, v = result.md$vi)

value

## 
## Fail-safe N Calculation Using the Rosenthal Approach 
## 
## Observed Significance Level: <.0001 
## Target Significance Level:   0.05 
## 
## Fail-safe N: 65

value$fsnum

## [1] 65

value$alpha  # Target Significance Level

## [1] 0.05

Other R Packages for Meta-Analysis

• Package metafor is the most comprehensive of currently available R packages for performing meta-analysis, but some may find its design overly complex (think iTunes)

• The package meta has a lot of overlap in provided methods, but it separates modeling functions by endpoint type

• The package rmeta only has DSL random effects modeling and no meta-regression modeling functions, which might be fine for some purposes

• The reliability of all of these packages is very good

library(meta)  # Package meta

Main Functions:

• metabin: Meta-analysis for binary outcome
• metacont: Meta-analysis for continuous outcome
• metareg: Meta-regression
• forest: Forest plot
• funnel: Funnel plot
• trimfill: Trim-and-fill method
• metabias: Test of asymmetry in funnel plot

dat.bcg$tn <- dat.bcg$tpos + dat.bcg$tneg
dat.bcg$cn <- dat.bcg$cpos + dat.bcg$cneg
result.or.meta <- metabin(event.e = tpos, n.e = tn, event.c = cpos, n.c = cn, 
	data = dat.bcg, 
	sm = "OR", 
	method = "Inverse", 
	method.tau = "REML")

names(result.or.meta)

##  [1] "event.e"      "n.e"          "event.c"      "n.c"         
##  [5] "studlab"      "TE"           "seTE"         "w.fixed"     
##  [9] "w.random"     "TE.fixed"     "seTE.fixed"   "lower.fixed" 
## [13] "upper.fixed"  "zval.fixed"   "pval.fixed"   "TE.random"   
## [17] "seTE.random"  "lower.random" "upper.random" "zval.random" 
## [21] "pval.random"  "k"            "Q"            "tau"         
## [25] "se.tau2"      "Q.CMH"        "sm"           "method"      
## [29] "sparse"       "incr"         "allincr"      "addincr"     
## [33] "allstudies"   "MH.exact"     "RR.cochrane"  "incr.e"      
## [37] "incr.c"       "level"        "level.comb"   "comb.fixed"  
## [41] "comb.random"  "hakn"         "df.hakn"      "method.tau"  
## [45] "tau.preset"   "TE.tau"       "method.bias"  "title"       
## [49] "complab"      "outclab"      "label.e"      "label.c"     
## [53] "label.left"   "label.right"  "call"         "warn"        
## [57] "print.byvar"  "print.CMH"    "version"

• w.fixed, w.random: Weight of individual studies

• TE.fixed, TE.random: Estimated overall treatment effect

• lower.fixed, upper.fixed: Lower and upper confidence intervals

• lower.random, upper.random: Lower and upper confidence intervals

• k: Number of studies

• tau: Estimated between-study variance

• Q: Heterogeneity statistic

summary(result.or.meta)

## Number of studies combined: k=13
## 
##                         OR          95%-CI       z  p.value
## Fixed effect model   0.647  [0.595; 0.702] -10.319 < 0.0001
## Random effects model 0.475  [0.330; 0.683]  -4.006 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.3378; H = 3.69 [3.04; 4.47]; I^2 = 92.6% [89.2%; 95%]
## 
## Test of heterogeneity:
##       Q d.f.  p.value
##  163.16   12 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - restricted maximum-likelihood estimator for tau^2

forest(result.or.meta)  # Default like Cochrane forest plot

metabias(result.or.meta, method = "rank")  # Rank-correlation test

## 
##  Rank correlation test of funnel plot asymmetry
## 
## data:  result.or.meta
## z = 0.122, p-value = 0.9029
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
##    ks se.ks 
##  2.00 16.39

No indication of asymmetry for OR analysis.

library(rmeta)  # rmeta package

Key Functions:

• meta.DSL: RE meta-analysis (Binary only)
• meta.MH: FE meta-analysis (Mantel-Haenszel)
• meta.summaries: Fixed/Random given ES and weights
• forestplot: Forest plot
• funnelplot: Funnel plot

• rmeta has the fewest features of the packages we have discussed.

• Focuses on methods for binary data.

• Only option for fixed-effects method is the Mantel-Haenszel method for combining ORs.

• Like Peto's OR, this is a FE model that can be advantageous for handling studies with rare events.

dat.bcg$tn <- with(dat.bcg, tpos + tneg)
dat.bcg$cn <- with(dat.bcg, cpos + cneg)
dat.bcg$tp <- with(dat.bcg, tpos/tn)
dat.bcg$cp <- with(dat.bcg, cpos/cn)

result.mh <- meta.MH(tn, cn, tp, cp, data = dat.bcg)

names(result.mh)

##  [1] "logOR"      "selogOR"    "logMH"      "selogMH"    "MHtest"    
##  [6] "het"        "call"       "names"      "conf.level" "statistic"

• logOR: Log odds ratio

• logMH: Estimated overall log OR

• selogMH: Standard error of overall log OR

• MHtest: Mantel-Haenszel \(\chi^2\)-test that OR=1

• het: Woolf’s chi-square for heterogeneity, df, p-value

summary(result.mh)

## Fixed effects ( Mantel-Haenszel ) meta-analysis
## Call: meta.MH(ntrt = tn, nctrl = cn, ptrt = tp, pctrl = cp, data = dat.bcg)
## ------------------------------------
##         OR (lower  95% upper)
##  [1,] 0.46       0  1.881e+05
##  [2,] 0.20       0  9.548e+05
##  [3,] 0.25       0  5.984e+07
##  [4,] 0.22       0  2.332e+13
##  [5,] 0.92       0  1.368e+14
##  [6,] 0.43       0  4.339e+02
##  [7,] 0.05       0  2.016e+15
##  [8,] 1.01       0  9.527e+15
##  [9,] 0.61       0  1.713e+17
## [10,] 0.25       0  9.272e+08
## [11,] 0.38       0  9.164e+17
## [12,] 1.46       0  4.151e+30
## [13,] 1.04       0  1.035e+30
## ------------------------------------
## Mantel-Haenszel OR =0.36 95% CI ( 0,47.42 )
## Test for heterogeneity: X^2( 12 ) = 0.03 ( p-value 1 )

c(exp(result.or$b), exp(result.mh$logMH))  # Compare with RE model

## [1] 0.4736 0.3635

result.mh$MHtest

## [1] 0.1825 0.6692

Comparing metafor, rmeta, meta

Overall

Calculated Effect Sizes

Synthesis Methods

Computed & Accesible Heterogeneity Measures

Computed & Accesible Confidence Intervals

Graphics

Some Advanced Topics

• No events in one or both treatment groups can create computational problems in standard meta-analytic approaches.

• Most proprietary programs handle this by adding 0.5 to zero counts.

• Mantel-Haenszel only requires correction if zero counts occur for same treatment arm in all studies.

• Peto's odds ratio only requires correction if zero counts in both arms of one or more studies.

• Always preferred but are rarely possible to undertake.

• More often only some studies have IPD and focusing only on these can introduce selection bias.

• If bias is not of concern, synthesize with hierarchical modeling.

• glmm and coxme packages are commonly used for hierarchical modeling in R.

• I have written the ipdmeta package for assessing the power of subgroup effects with IPD meta-analysis.

• Great care is needed in assessing compatibility of effects when considering meta-analysis of non-randomized studies.

• When "treatment" is not randomized, greater heterogeneity between studies is expected.

• The reason for this is that outcomes will be more sensitive to the sample characteristics and adjustment methods a study has used.

• In general, one should use the most fully-adjusted measure of effect.

• We have focused on the summary of effect sizes.

• It is also possible to summarize evidence across trials by combining confidence intervals, the so-called "confidence distribution method".

• Some advantages of combining CDs:

  • Robust to outlying studies
  • Yields exact CI for combining 2 by 2 tables, even with rare events

• Standard FE or RE methodology can be applied to perform meta-analyses of genomic data (e.g., GWAS, microarray, etc.).

• But these data introduce further issues:

  • Missing data
  • Platform discrepancies between studies
  • Multiplicity
  • Efficiency

• There may be multiple endpoints of interest.

• A multivariate meta-analysis combines estimates of multiple outcomes, accounting for their correlation.

• Multivariate meta-regression can also incorporate the effects of study-level predictors.

• A disease may have multiple trial-tested treatments.

• In general, only a subset of treatments will have been considered in any given trial.

• When there is interest in comparing the efficacy among all trials, when not all have been directly compared in available trials, a network (or mixed-treatment) meta-analysis can be performed.

• Proposed network meta-analysis methods usually involve Bayesian approaches and require careful assessment of consistency in treatment comparisons.

• Bayesian meta-analysis focuses on estimating a posterior distribution of effect rather than a summary effect estimate.

• A Bayesian framework can be advantageous for:

  • Accounting for possible bias
  • Making mixed-treatment comparisons
  • Conducting multivariate analyses

Packages For Advanced Meta-Analyses

Remember...

1. Not to begin here!

2. That meta-analytic summaries are all about weighted averages

3. That evaluating bias and heterogeneity are essential steps of meta-analysis

4. You now have a basic knowledge of how to use multiple R packages to perform conventional meta-analyses

Resources For Systematic Reviews

• Cochrane Collaboration Handbook

• Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA Statement)

• Meta-analysis of Observational Studies in Epidemiology (MOOSE Statement)

Further Reading

Overview (1)

• Cochrane Handbook
• CRAN Task View on Meta-Analysis.
• Chen D-G, Peace KE. Applied meta-analysis with R. Boca Raton, Florida: Taylor & Francis Group; 2013.
• Ellis PD. The essential guide to effect sizes : statistical power, meta-analysis, and the interpretation of research results. Cambridge ; New York: Cambridge University Press; 2010.
• Hartung J, Knapp G, Sinha BK. Statistical meta-analysis with applications. Hoboken, N.J.: Wiley; 2008.

Overview (2)

• Hartung J, Knapp G. On tests of the overall treatment effect in meta-analysis with normally distributed responses. Stat Med 2001;20:1771-82.
• Hedges LV, Olkin I. Statistical methods for meta-analysis. Orlando: Academic Press; 1985.
• Hunter JE, Schmidt FL. Methods of meta-analysis : correcting error and bias in research findings. 2nd ed. Thousand Oaks, Calif.: Sage; 2004.

Heterogeneity

• DerSimonian R, Laird N. Meta-analysis in clinical trials. Control Clin Trials 1986;7:177-88.
• Hardy RJ, et al. Detecting and describing heterogeneity in meta-analysis. Stat Med 1998;17:841-56.
• Higgins JP, Thompson SG. Quantifying heterogeneity in a meta-analysis. Stat Med 2002;21:1539-58.
• Knapp G, Biggerstaff BJ, Hartung J. Assessing the amount of heterogeneity in random-effects meta-analysis. Biom J 2006;48:271-85.
• Viechtbauer W. Hypothesis tests for population heterogeneity in meta-analysis. Br J Math Stat Psychol 2007;60:29-60.
• Viechtbauer W. Confidence intervals for the amount of heterogeneity in meta-analysis. Stat Med 2007;26:37-52.

Advanced Topics (1)

• Cai T, Parast L, Ryan L. Meta-analysis for rare events. Stat Med 2010;29:2078-89.
• Dukic V, Gatsonis C. Meta-analysis of diagnostic test accuracy assessment studies with varying number of thresholds. Biometrics 2003;59:936-46.
• Gasparrini A, Armstrong B, Kenward MG. Multivariate meta-analysis for non-linear and other multi-parameter associations. Stat Med 2012;31:3821-39.
• Kovalchik SA. Aggregate-data estimation of an individual patient data linear random effects meta-analysis with a patient covariate-treatment interaction term. Biostatistics 2013;14:273-83.
• Kovalchik SA, Cumberland WG. Using aggregate data to estimate the standard error of a treatment-covariate interaction in an individual patient data meta-analysis. Biom J 2012;54:370-84.

Advanced Topics (1)

• Lin DY, Zeng D. On the relative efficiency of using summary statistics versus individual-level data in meta-analysis. Biometrika 2010;97:321-32.
• Singh K, Xie M, Strawdermann. Combining information from independent sources through confidence distribution. Annals of Statistics 2005;33:159-183.
• Thompson SG, Higgins JP. How should meta-regression analyses be undertaken and interpreted? Stat Med 2002;21:1559-73.
• van Houwelingen HC, Arends LR, Stijnen T. Advanced methods in meta-analysis: multivariate approach and meta-regression. Stat Med 2002;21:589-624.

Applications

• 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:698-702.
• Nissen SE, Wolski K. Effect of rosiglitazone on the risk of myocardial infarction and death from cardiovascular causes. N Engl J Med 2007;356:2457-71.