by A recent study (Langan et al., 2019) suggests that the two-step DerSimonian and Laird (DL2; Dersimonian & Knacker, 2007) estimator for tau-squared displayed the best properties for random-effects models of meta-analyses for continuous data.
You have the option of calculate the traditional Wald-type confidence intervals and p-values or use the Hartung-Knapp/Sidik-Jonkman (HKSJ) method to calculate the confidence intervals and p-values. Langan et al. (2009) recommend using the HKSJ methods as the default unless your meta-analysis contains fewer than five studies of unequal sample sizes or you do not have heterogeneity in effect sizes.
Finally, you can calculate the I-squared statistic and its confidence intervals using the standard method (see Higgins et al., 2003) or the Q-profile method.
Dependencies
The function requires the dplyr and rlang packages.
Arguments
The function takes the following arguments:
yi: Effect size vi: Variance of effect size data: data.frame or tibble lab: Study name ci: Confidence interval (default = .95) Q.profile: Use Q.profile method for I-squared (default = FALSE) pi: Calculate prediction interval (default = FALSE) hjsj: Use HKSJ adjustment (default = FALSE)
Function
# Meta-Analysis Function --------------------------------------------------
## Summary Effect Size and its Standard Error using DL method (Eqs. 1-3)
mw_est <- function(yi,
vi,
data,
lab,
method = c("DL", "DL2"),
ci = .95,
Q.profile = FALSE,
pi = FALSE,
hksj = FALSE) {
## General Method-of-Moments Estimate for Tau-Squared (Eq. 6)
.tau2MM <- function(yi, vi, ai) {
yw <- sum(ai * yi) / sum(ai)
tau2 <- (sum(ai * (yi - yw)^2) - (sum(ai * vi) - sum(ai^2 * vi) / sum(ai))) / (sum(ai) - sum(ai^2) / sum(ai))
# Return only non-negative values
max(0, tau2)
}
# Keep complete cases
data <- data %>%
select(c({{lab}}, {{ yi }}, {{ vi }})) %>%
filter_all(all_vars(!is.na(.)))
# Extract variables
yi <- pull(data, {{ yi }})
vi <- pull(data, {{ vi }})
# Calculate DL tau2
ai <- 1 / vi
tau2 <- .tau2MM(yi = yi, vi = vi, ai = ai)
# Calculate DL2 if requested
if (method == "DL2") {
ai <- 1 / (vi + tau2)
tau2 <- .tau2MM(yi = yi, vi = vi, ai = ai)
}
## Calculate inverse variance weights
wi <- (1 / (vi + tau2)) # Eq. 2
wy <- wi * yi
mw <- sum(wy) / sum(wi) # Eq. 1
# Get number of studies
k <- length(yi)
# Degrees of freedom
df <- k - 1
### Heterogeneity
## Cochran's Q-statistic
# Use DL inverse-variance weight (ai = 1 / vi)
ai <- 1 / vi
# Calculate Q-statistic
Q <- sum(ai * ((yi - (sum(ai * yi) / sum(ai)))^2))
# Calculate p-value
Q.p <- pchisq(Q, df = df, lower.tail = FALSE)
z.cdf <- qnorm(1 - ((1 - ci) / 2))
# Method used to calculate heterogeneity and confidence intervals
if (Q.profile) {
## Q-profile method to calculate heterogeneity and confidence
## intervals (Eq. 14)
## Calculate I-squared and H
# a. Use DL inverse-variance weight (ai = 1 / vi)
ai <- 1 / vi
# b. Calculate V
V <- (k - 1) * sum(ai) / (sum(ai)^2 - sum(ai^2))
# c. Calculate I-squared
I.sq <- 100 * tau2 / (tau2 + V)
# d. Calculate H-squared
H <- (tau2 + V) / V
# Calculate confidence intervals if k > 2
if (k > 2) {
## Generalised Q-statistics bounds
Q.lb <- qchisq((1 - ci) / 2, df, lower = TRUE)
Q.ub <- qchisq((1 - ci) / 2, df, lower = FALSE)
## Lower confidence intervals
# Initialise values
tau2.tmp <- 0
F.tau2 <- 1
while (F.tau2 > 0) {
# i. calculate weight for tau2
W <- 1 / (vi + tau2.tmp)
yW <- sum(yi * W) / sum(W)
# ii. Calculate delta tau2
F.tau2 <- sum(W * (yi - yW)^2) - Q.ub
F.tau2 <- max(F.tau2, 0)
delta.tau2 <- F.tau2 / sum((W^2) * (yi - yW)^2)
# iii. Next iterative step
tau2.tmp <- tau2.tmp + delta.tau2
}
tau2.tmp <- max(0, tau2.tmp)
I.lower <- 100 * tau2.tmp / (tau2.tmp + V)
H.lower <- (tau2.tmp + V) / V
## Upper confidence intervals
# Initialise values
tau2.tmp <- 0
F.tau2 <- 1
while (F.tau2 > 0) {
# i. calculate weight for tau2
W <- 1 / (vi + tau2.tmp)
yW <- sum(yi * W) / sum(W)
# ii. Calculate delta tau2
F.tau2 <- sum(W * (yi - yW)^2) - Q.lb
F.tau2 <- max(F.tau2, 0)
delta.tau2 <- F.tau2 / sum((W^2) * (yi - yW)^2)
# iii. Next iterative step
tau2.tmp <- tau2.tmp + delta.tau2
}
tau2.tmp <- max(0, tau2.tmp)
I.upper <- 100 * tau2.tmp / (tau2.tmp + V)
H.upper <- (tau2.tmp + V) / V
} else {
I.lower <- NA
I.upper <- NA
H.lower <- NA
H.upper <- NA
pi.res <- NULL
}
} else {
## I-squared
I.sq <- max(0, 100 * ((Q - df) / Q))
## H
H <- sqrt(Q / df)
# Calculate the SE of ln(H) aka B
if (k > 2) {
if (Q > k) {
B <- 0.5 * ((log(Q) - log(df)) / (sqrt(2 * Q) - sqrt(2 * df - 1)))
} else {
B <- sqrt(1 / ((2 * df - 1) * (1 - (1 / (3 * (df - 1)^2)))))
}
# H confidence intervals
H.lower <- max(1, exp(log(H) - z.cdf * B))
H.upper <- exp(log(H) + z.cdf * B)
# I-Squared
I.L <- exp(0.5 * log(Q / df) - (z.cdf * B))
I.U <- exp(0.5 * log(Q / df) + (z.cdf * B))
I.lower <- max(0, 100 * ((I.L^2 - 1) / I.L^2))
I.upper <- max(0, 100 * ((I.U^2 - 1) / I.U^2))
} else {
I.lower <- NA
I.upper <- NA
H.lower <- NA
H.upper <- NA
pi.res <- NULL
}
}
# Save all heterogeneity statistics in a list
het.test <- list(
Q = Q, Q.p = Q.p,
tau2 = tau2,
I.sq = I.sq, I.lb = I.lower, I.ub = I.upper,
H = H, H.lb = H.lower, H.ub = H.upper
)
# HKSJ adjustment
if (hksj & k >= 3) {
# Calculate HKSJ-adjusted standard error
sew <- sqrt(sum(wi * (yi - mw)^2) / (df * sum(wi)))
# Calculate t-value
t <- mw / sew
# Calculate p-value based on t-statistics and k - 1 df
t.p <- 2 * pt(-abs(t), df = k - 1)
# Output test for summary effect size
mw.test <- list(t = t, p = t.p)
# Confidence Intervals
mw.lower <- mw - sew * qt(1 - ((1 - ci) / 2), df = df)
mw.upper <- mw + sew * qt(1 - ((1 - ci) / 2), df = df)
# Wald-type statistics
} else {
# Calculate standard error
sew <- 1 / sqrt(sum(wi)) # Eq. 3
# Calculate z-statistic and p-value
z <- mw / sew
z.p <- 2 * pnorm(-abs(z))
# Confidence Intervals
mw.lower <- mw - z.cdf * sew
mw.upper <- mw + z.cdf * sew
# Output test for summary effect size
mw.test <- list(z = z, p = z.p)
}
# Prediction Interval (based on t-distribution)
if (all(k > 2, pi)) {
if (hksj) {
pi.lower <- mw - sqrt(tau2 + sew^2) * qt(1 - ((1 - ci) / 2), df = df)
pi.upper <- mw + sqrt(tau2 + sew^2) * qt(1 - ((1 - ci) / 2), df = df)
} else {
pi.lower <- mw - sqrt(tau2 + sew^2) * z.cdf
pi.upper <- mw + sqrt(tau2 + sew^2) * z.cdf
}
# Output prediction interval
pi.res <- list(pi.lb = pi.lower, pi.ub = pi.upper)
} else {
pi.res <- NULL
}
# Print data set
cat(
"n", "Data used for meta-analysis",
"n",
"n",
sep = ""
)
print(data)
cat(
"n",
sep = ""
)
# Output results as a list
return(
c(
method = method,
est = mw, se.est = sew,
mw.test, ci.lb = mw.lower, ci.ub = mw.upper,
het.test, pi.res
)
)
}
Example
We are going to recreate the analyses in DerSimonian and Knacker (2007).
First we recreate the CLASP dataset and then calculate the log odds ratio and its variance:
# CLASP Dataset
dat_clasp <- tibble(
study = c("Baum", "Anand", "CLASP", "Bradley", "Llovet", "Teo"),
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_clasp <- metafor::escalc(
measure = "OR",
ai = ai, n1i = n1i,
ci = ci,
n2i = n2i, data = dat_clasp
)
Your data should look like the following:
study n1i n2i ai ci yi vi 1 Baum 156 74 5 8 -1.2976 0.3468 2 Anand 303 303 5 17 -1.2649 0.2657 3 CLASP 565 477 12 9 0.1208 0.1984 4 Bradley 1570 1565 69 94 -0.3294 0.0265 5 Llovet 103 105 9 11 -0.2007 0.2233 6 Teo 4659 4650 313 352 -0.1285 0.0065
Next we run the function. I have requested the DL2 estimator and Q.profile method to recreate the results of DerSimonian and Knacker (2007), who reported a tau-squared of 0.11, a summary effect size (log odds ratio) of -0.37 (SE = 0.19).
mw_est( yi = yi, vi = vi, data = dat_clasp, lab = study, method = "DL2", pi = FALSE, Q.profile = TRUE, hksj = FALSE )
You will get the following output, which is consistent with the results of DerSimonian and Knacker (2007).
Data used for meta-analysis
study yi vi
1 Baum -1.2976287 0.34677403
2 Anand -1.2648771 0.26567574
3 CLASP 0.1207924 0.19838951
4 Bradley -0.3293774 0.02647708
5 Llovet -0.2006707 0.22329680
6 Teo -0.1285341 0.00649856
$method
[1] "DL2"
$est
[1] -0.3654628
$se.est
[1] 0.1900638
$z
[1] -1.922843
$p
[1] 0.05449979
$ci.lb
[1] -0.737981
$ci.ub
[1] 0.007055387
$Q
[1] 9.677867
$Q.p
[1] 0.08489455
$tau2
[1] 0.1058598
$I.sq
[1] 64.72927
$I.lb
[1] 0
$I.ub
[1] 97.18495
$H
[1] 2.835212
$H.lb
[1] 1
$H.ub
[1] 35.5233
A future post will present a Shiny application that presents a nicely-formatted table and forest plot using the above function.
References
Cornell, J. E., Mulrow, C. D., Localio, R., Stack, C. B., Meibohm, A. R., Guallar, E., & Goodman, S. N. (2014). Random-effects meta-analysis of inconsistent effects: A time for change. Annals of Internal Medicine, 160(4), 267–270. https://doi.org/10.7326/M13-2886
DerSimonian, R., & Kacker, R. (2007). Random-effects model for meta-analysis of clinical trials: An update. Contemporary Clinical Trials, 28(2), 105–114. https://doi.org/10.1016/j.cct.2006.04.004
Higgins, J. P. T., Thompson, S. G., Deeks, J. J., & Altman, D. G. (2003). Measuring inconsistency in meta-analyses. BMJ, 327(7414), 557–560. https://doi.org/10.1136/bmj.327.7414.557
Langan, D., Higgins, J. P. T., Jackson, D., Bowden, J., Veroniki, A. A., Kontopantelis, E., … Simmonds, M. (2019). A comparison of heterogeneity variance estimators in simulated random-effects meta-analyses. Research Synthesis Methods, 10(1), 83–98. https://doi.org/10.1002/jrsm.1316
