Table of Contents
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.
Example
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) dat
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, se)) 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.
References
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.