The metafor Package

A Meta-Analysis Package for R

User Tools

Site Tools


DerSimonian and Kacker (2007)

The Methods and Data

In this paper, the authors describe a variety of methods for estimating the amount of heterogeneity under a random-effects model. In addition to the well-known DerSimonian-Laird and Cochran estimators (the latter is also known as the Hedges or variance component estimator), the author also describe the Paule-Mandel estimator, a two-step Cochran estimator, and a two-step DerSimonian-Laird estimator.

Below, I show how all of these estimators can be obtained for one of the illustrate datasets used in the paper (from CLASP, 1994). For completeness sake, I'll also throw in the other estimator available in the metafor package, including the empirical Bayes estimator (which can actually be shown to be identical to the Paule-Mandel estimator), the maximum-likelihood and restricted maximum-likelihood estimators, the Hunter-Schmidt estimator, the Sidik-Jonkman estimator, and the improved Sidik-Jonkman estimator that uses the Cochran estimator as a starting value.


The dataset from the CLASP (1994) study can be created with:

n1i <- c(156, 303, 565, 1570, 103, 4659)
n2i <- c( 74, 303, 477, 1565, 105, 4650)
ai  <- c(  5,   5,  12,   69,   9,  313)
ci  <- c(  8,  17,   9,   94,  11,  352)
dat <- escalc(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i)
       yi     vi
1 -1.2976 0.3468
2 -1.2649 0.2657
3  0.1208 0.1984
4 -0.3294 0.0265
5 -0.2007 0.2233
6 -0.1285 0.0065

The values given under yi are the observed log odds ratios and the values under vi are the corresponding sampling variances.

After loading the metafor package with library(metafor), we can now fit a random-effects model to these data using the various heterogeneity estimators with:

res.PM  <- rma(yi, vi, method="PM", data=dat)
res.CA  <- rma(yi, vi, method="HE", data=dat)
res.DL  <- rma(yi, vi, method="DL", data=dat)
res.CA2 <- rma(yi, vi, method="GENQ", weights=1/(vi + res.CA$tau2), data=dat)
res.DL2 <- rma(yi, vi, method="GENQ", weights=1/(vi + res.DL$tau2), data=dat)
res.CA2 <- rma(yi, vi, tau2=res.CA2$tau2, data=dat)
res.DL2 <- rma(yi, vi, tau2=res.DL2$tau2, data=dat)
res.EB   <- rma(yi, vi, method="EB",   data=dat)
res.ML   <- rma(yi, vi, method="ML",   data=dat)
res.REML <- rma(yi, vi, method="REML", data=dat)
res.HS   <- rma(yi, vi, method="HS",   data=dat)
res.SJ   <- rma(yi, vi, method="SJ",   data=dat)
res.SJ2  <- rma(yi, vi, method="SJ",   data=dat, control=list(tau2.init=res.CA$tau2))

Objects res.PM, res.CA, and res.DL include the results when using the Paule-Mandel, Cochran, and DerSimonian-Laird estimators (the Cochran estimator is referred to as the Hedges estimator in the metafor package).

The next two models (res.CA2 and res.DL2) use the generalized Q-statistic estimator of $\tau^2$ (referred to as the general method-of-moments estimator in the paper), where the weights are set equal to the usual (random-effects model) inverse-variance weights, once using the Cochran and once using the DerSimonian-Laird estimator values. This will provide us with the two-step Cochran and two-step DerSimonian-Laird estimator values. However, note that the weights specified are used not only to estimate $\tau^2$, but also to estimate $\mu$, the average effect. If we now want to use the usual inverse-variance weights but using the two-step estimator values, we have to fit each model one more time, but this time with specified values for $\tau^2$ as obtained from the two-step procedure.

Next, we obtain the results using the empirical Bayes, maximum-likelihood and restricted maximum-likelihood estimators, the Hunter-Schmidt estimator, and the Sidik-Jonkman estimator.

Finally, to obtain the improved Sidik-Jonkman estimator, we must use the control argument to set the initial estimate of $\tau^2$ to that of the Cochran estimator, which then gets updated using the method described by Sidik and Jonkman.

Next, we can combine the model objects into a list with:

res.all <- list(res.PM, res.CA, res.DL, res.CA2, res.DL2, res.EB, res.ML, res.REML, res.HS, res.SJ, res.SJ2)

Finally, we can extract the estimated values of $\tau$ (i.e., the square-root of the $\tau^2$ estimates), the estimated values of $\mu$, and the corresponding standard errors with:

results <- rbind(
tau = sapply(res.all, function(x) sqrt(x$tau2)),
mu  = sapply(res.all, coef),
se  = sapply(res.all, function(x) sqrt(vcov(x))))
colnames(results) <- c("PM", "CA", "DL", "CA2", "DL2", "EB", "ML", "REML", "HS", "SJ", "SJ2")
round(t(results), 4)

which yields:

        tau      mu     se
PM   0.3681 -0.3811 0.2060
CA   0.4410 -0.4035 0.2327
DL   0.2323 -0.3240 0.1540
CA2  0.3831 -0.3861 0.2115
DL2  0.3254 -0.3655 0.1901
EB   0.3681 -0.3811 0.2060
ML   0.0023 -0.1974 0.0694
REML 0.1843 -0.2980 0.1343
HS   0.1330 -0.2666 0.1125
SJ   0.4572 -0.4079 0.2386
SJ2  0.4084 -0.3941 0.2208

The values for the Paule-Mandel, Cochran, DerSimonian-Laird, two-step Cochran, and two-step DerSimonian-Laird estimators are given in Tables 3 and 4 in the paper and are exactly the same as given above. Note also that using the empirical Bayes estimator yields the exact same results as when using the Paule-Mandel estimator.


DerSimonian, R., & Kacker, R. (2007). Random-effects model for meta-analysis of clinical trials: An update. Contemporary Clinical Trials, 28(2), 105–114.

CLASP (Collaborative Low-dose Aspirin Study in Pregnancy) Collaborative Group (1994). CLASP: A randomised trial of low-dose aspirin for the prevention and treatment of pre-eclampsia among 9364 pregnant women. Lancet, 343(8898), 619–629.

analyses/dersimonian2007.txt · Last modified: 2014/12/15 14:26 by Wolfgang Viechtbauer