Supplementary R code for Statistical
Methods for Meta-Analysis
Supplementary R code for Statistical
Methods for Meta-Analysis
Crippa Alessio
Last compiled on April 30, 2024
Description
This page is a supplementary material for the lecture Statistical Methods for Meta-Analysis, in the course Systematic Review and Meta-Analysis.
It contains the R code for exemplifying the theoretical aspects presented in the slides using real worked example. See the homepage of the metafor package for a much more comprehensive collection of analysis examples.
To reproduce the results you need install the metafor
package by typing install.package("metafor"). Load the
package
library(metafor)
library(dosresmeta)
library(tidyverse)Fixed-effect model
Data (Hine et al.)
The example from Hine et al. (1989) contains the results from 6 studies on the mortality risk due to prophylactic use of lidocaine after a heart attack.
dat.hine1989 study source n1i n2i ai ci
1 1 Chopra et al. 39 43 2 1
2 2 Mogensen 44 44 4 4
3 3 Pitt et al. 107 110 6 4
4 4 Darby et al. 103 100 7 5
5 5 Bennett et al. 110 106 7 3
6 6 O'Brien et al. 154 146 11 4Effect size
Among different possibilities, we choose to compute risk differences as preferable effect size.
hine1989_rd <- escalc(measure = "RD", n1i = n1i, n2i = n2i, ai = ai, ci = ci,
data = dat.hine1989)
hine1989_rd
study source n1i n2i ai ci yi vi
1 1 Chopra et al. 39 43 2 1 0.0280 0.0018
2 2 Mogensen 44 44 4 4 0.0000 0.0038
3 3 Pitt et al. 107 110 6 4 0.0197 0.0008
4 4 Darby et al. 103 100 7 5 0.0180 0.0011
5 5 Bennett et al. 110 106 7 3 0.0353 0.0008
6 6 O'Brien et al. 154 146 11 4 0.0440 0.0006 Fixed-effect model
fit_hine1989 <- rma(yi = yi, vi = vi, data = hine1989_rd, method = "FE")
summary(fit_hine1989)
Fixed-Effects Model (k = 6)
logLik deviance AIC BIC AICc
14.2466 0.8597 -26.4932 -26.7014 -25.4932
I^2 (total heterogeneity / total variability): 0.00%
H^2 (total variability / sampling variability): 0.17
Test for Heterogeneity:
Q(df = 5) = 0.8597, p-val = 0.9731
Model Results:
estimate se zval pval ci.lb ci.ub
0.0294 0.0131 2.2531 0.0243 0.0038 0.0551 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Forest plot
forest(fit_hine1989,
ilab = cbind(dat.hine1989$ai, dat.hine1989$n1i, dat.hine1989$ci, dat.hine1989$n2i),
ilab.xpos = c(-.3, -.25, -.2, -.15), cex = .75, header = "Author(s)",
mlab = "Summary (fixed-effect model)", xlim = c(-.5, .4), slab = dat.hine1989$source)
text(c(-.3, -.25, -.2, -.15), 8, c("d", "n", "d", "n"), cex =.75, font=2)
text(c(-.27, -.17), 8.5, c("Lidocaine", "Control"), cex = .75, font = 2)
Random-effects model
Data (Berkey et al.)
The example from Berkey et al. (1995) contains the results with the BCG vaccine dataset (Colditz et al., 1994).
dat.bcg trial author year tpos tneg cpos cneg ablat alloc
1 1 Aronson 1948 4 119 11 128 44 random
2 2 Ferguson & Simes 1949 6 300 29 274 55 random
3 3 Rosenthal et al 1960 3 228 11 209 42 random
4 4 Hart & Sutherland 1977 62 13536 248 12619 52 random
5 5 Frimodt-Moller et al 1973 33 5036 47 5761 13 alternate
6 6 Stein & Aronson 1953 180 1361 372 1079 44 alternate
7 7 Vandiviere et al 1973 8 2537 10 619 19 random
8 8 TPT Madras 1980 505 87886 499 87892 13 random
9 9 Coetzee & Berjak 1968 29 7470 45 7232 27 random
10 10 Rosenthal et al 1961 17 1699 65 1600 42 systematic
11 11 Comstock et al 1974 186 50448 141 27197 18 systematic
12 12 Comstock & Webster 1969 5 2493 3 2338 33 systematic
13 13 Comstock et al 1976 27 16886 29 17825 33 systematicEffect size
We calculate the log risk ratios and corresponding sampling variances
bcg_rr <- escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg,
slab = paste(author, year, sep=", "))Random-effects model
fit_bcg <- rma(yi, vi, data = bcg_rr, method = "REML")
summary(fit_bcg)
Random-Effects Model (k = 13; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
-12.2024 24.4047 28.4047 29.3746 29.7381
tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664)
tau (square root of estimated tau^2 value): 0.5597
I^2 (total heterogeneity / total variability): 92.22%
H^2 (total variability / sampling variability): 12.86
Test for Heterogeneity:
Q(df = 12) = 152.2330, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
-0.7145 0.1798 -3.9744 <.0001 -1.0669 -0.3622 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretation on the natural scale (exponentiating the log results)
predict(fit_bcg, transf = exp)
pred ci.lb ci.ub pi.lb pi.ub
0.4894 0.3441 0.6962 0.1546 1.5490 Forest plot
forest(fit_bcg, atransf = exp, at = log(c(.05, .25, 1, 4)), xlim = c(-16, 6),
ilab = cbind(dat.bcg$tpos, dat.bcg$tneg, dat.bcg$cpos, dat.bcg$cneg),
ilab.xpos = c(-9.5, -8, -6, -4.5), cex = .75, header = "Author(s) and Year",
mlab = "")
text(c(-9.5, -8, -6, -4.5), 15, c("TB+", "TB-", "TB+", "TB-"), cex = .75, font = 2)
text(c(-8.75, -5.25), 16, c("Vaccinated", "Control"), cex = .75, font = 2)
# add text with Q-value, dfs, p-value, and I^2 statistic
text(-16, -1, pos = 4, cex = 0.75,
bquote(paste("RE Model (Q = ",
.(formatC(fit_bcg$QE, digits = 2, format = "f")), ", df = ",
.(fit_bcg$k - fit_bcg$p), ", p = ",
.(formatC(fit_bcg$QEp, digits=2, format="f")), "; ", I^2, " = ",
.(formatC(fit_bcg$I2, digits=1, format="f")), "%)")))
Meta regression
fit_mr_bcg <- rma(yi, vi, mods = ~ ablat, data = bcg_rr, method = "REML")
fit_mr_bcg
Mixed-Effects Model (k = 13; tau^2 estimator: REML)
tau^2 (estimated amount of residual heterogeneity): 0.0764 (SE = 0.0591)
tau (square root of estimated tau^2 value): 0.2763
I^2 (residual heterogeneity / unaccounted variability): 68.39%
H^2 (unaccounted variability / sampling variability): 3.16
R^2 (amount of heterogeneity accounted for): 75.62%
Test for Residual Heterogeneity:
QE(df = 11) = 30.7331, p-val = 0.0012
Test of Moderators (coefficient 2):
QM(df = 1) = 16.3571, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2515 0.2491 1.0095 0.3127 -0.2368 0.7397
ablat -0.0291 0.0072 -4.0444 <.0001 -0.0432 -0.0150 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fit_mr_bcg, fit_bcg)
df AIC BIC AICc logLik LRT pval QE tau^2 R^2
Full 3 22.1746 23.3683 25.6032 -8.0873 30.7331 0.0764
Reduced 2 28.4047 29.3746 29.7381 -12.2024 8.2301 0.0041 152.2330 0.3132 75.6245% # calculate predicted risk ratios for 0 to 60 degrees absolute latitude
preds <- data.frame(predict(fit_mr_bcg, newmods = c(0:60), transf = exp))
# radius of points will be proportional to the inverse standard errors
# hence the area of the points will be proportional to inverse variances
size <- 1 / sqrt(bcg_rr$vi)
size <- size / max(size)
library(ggplot2)
library(scales)
p_bcg_rr <- ggplot(bcg_rr, aes(ablat, exp(yi), text = author)) +
geom_point(aes(size = size)) +
geom_text(aes(label = trial), hjust = 0, nudge_x = 3, size = 3) +
geom_hline(yintercept = 1, linetype = "dotted") +
geom_line(data = preds, aes(x = 0:60, y = pred, text = "line")) +
geom_ribbon(data = preds, aes(x = 0:60, y = pred, ymin = ci.lb, ymax = ci.ub, text = "line"), alpha = .2) +
geom_point(data = bcg_rr, aes(ablat, exp(yi), )) +
scale_y_continuous(trans = "log", breaks = pretty_breaks()) +
labs(x = "Absolute Latitude", y = "Risk Ratio") +
theme_classic()
p_bcg_rr
# make it interactive
library(plotly)
ggplotly(p_bcg_rr)Error in train(..., self = self): unused argument (list("Absolute Latitude", "Risk Ratio", "author", "size", "trial", "yintercept", "ci.lb", "ci.ub"))Heterogeneity
Data (Viechtbauer et al.)
The example from Viechtbauer et al. (2007) contains the results from 9 studies on the effectiveness of diuretics in pregnancy for preventing pre-eclampsia (Collins et al. (1985)).
dat.collins1985b <- dat.collins1985b[, 1:7]
dat.collins1985b id author year pre.nti pre.nci pre.xti pre.xci
1 1 Weseley & Douglas 1962 131 136 14 14
2 2 Flowers et al. 1962 385 134 21 17
3 3 Menzies 1964 57 48 14 24
4 4 Fallis et al. 1964 38 40 6 18
5 5 Cuadros & Tatum 1964 1011 760 12 35
6 6 Landesman et al. 1965 1370 1336 138 175
7 7 Kraus et al. 1966 506 524 15 20
8 8 Tervila & Vartiainen 1971 108 103 6 2
9 9 Campbell & MacGillivray 1975 153 102 65 40Effect size
We calculate the log risk ratios and corresponding sampling variances
collins1985b_OR <- escalc(measure = "OR", ai = pre.xti, n1i = pre.nti,
ci = pre.xci, n2i = pre.nci, data = dat.collins1985b)
collins1985b_OR
id author year pre.nti pre.nci pre.xti pre.xci yi vi
1 1 Weseley & Douglas 1962 131 136 14 14 0.0418 0.1596
2 2 Flowers et al. 1962 385 134 21 17 -0.9237 0.1177
3 3 Menzies 1964 57 48 14 24 -1.1221 0.1780
4 4 Fallis et al. 1964 38 40 6 18 -1.4733 0.2989
5 5 Cuadros & Tatum 1964 1011 760 12 35 -1.3910 0.1143
6 6 Landesman et al. 1965 1370 1336 138 175 -0.2969 0.0146
7 7 Kraus et al. 1966 506 524 15 20 -0.2615 0.1207
8 8 Tervila & Vartiainen 1971 108 103 6 2 1.0888 0.6864
9 9 Campbell & MacGillivray 1975 153 102 65 40 0.1353 0.0679 Random-effects model
collins1985b_DL <- rma(yi, vi, data = collins1985b_OR, method = "DL")
summary(collins1985b_DL)
Random-Effects Model (k = 9; tau^2 estimator: DL)
logLik deviance AIC BIC AICc
-9.4686 20.8785 22.9371 23.3316 24.9371
tau^2 (estimated amount of total heterogeneity): 0.2297 (SE = 0.1976)
tau (square root of estimated tau^2 value): 0.4793
I^2 (total heterogeneity / total variability): 70.66%
H^2 (total variability / sampling variability): 3.41
Test for Heterogeneity:
Q(df = 8) = 27.2649, p-val = 0.0006
Model Results:
estimate se zval pval ci.lb ci.ub
-0.5168 0.2037 -2.5367 0.0112 -0.9160 -0.1175 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Prediction intervals
predict(collins1985b_DL, transf = exp)
pred ci.lb ci.ub pi.lb pi.ub
0.5964 0.4001 0.8891 0.2149 1.6552 Confidence intervals
confint(collins1985b_DL, digits = 2)
estimate ci.lb ci.ub
tau^2 0.23 0.07 2.20
tau 0.48 0.27 1.48
I^2(%) 70.66 43.12 95.85
H^2 3.41 1.76 24.09 Alternative estimators
collins1985b_ML <- rma(yi, vi, data = collins1985b_OR, method = "ML")
collins1985b_REML <- rma(yi, vi, data = collins1985b_OR, method = "REML")
# estimate of tau^2
round(c(DL = collins1985b_DL$tau2, ML = collins1985b_ML$tau2, REML = collins1985b_REML$tau2), 2) DL ML REML
0.23 0.24 0.30 # I^2
round(c(DL = collins1985b_DL$I2, ML = collins1985b_ML$I2, REML = collins1985b_REML$I2), 2) DL ML REML
70.66 71.44 75.92 Sensitivity analysis
Data (Hackshaw et al.)
dat.hackshaw1998
study author year country design cases or or.lb or.ub yi vi
1 1 Garfinkel 1981 USA cohort 153 1.18 0.90 1.54 0.1655 0.0188
2 2 Hirayama 1984 Japan cohort 200 1.45 1.02 2.08 0.3716 0.0330
3 3 Butler 1988 USA cohort 8 2.02 0.48 8.56 0.7031 0.5402
4 4 Cardenas 1997 USA cohort 150 1.20 0.80 1.60 0.1823 0.0313
5 5 Chan 1982 Hong Kong case-control 84 0.75 0.43 1.30 -0.2877 0.0797
6 6 Correa 1983 USA case-control 22 2.07 0.81 5.25 0.7275 0.2273
7 7 Trichopolous 1983 Greece case-control 62 2.13 1.19 3.83 0.7561 0.0889
8 8 Buffler 1984 USA case-control 41 0.80 0.34 1.90 -0.2231 0.1927
9 9 Kabat 1984 USA case-control 24 0.79 0.25 2.45 -0.2357 0.3390
10 10 Lam 1985 Hong Kong case-control 60 2.01 1.09 3.72 0.6981 0.0981
11 11 Garfinkel 1985 USA case-control 134 1.23 0.81 1.87 0.2070 0.0456
12 12 Wu 1985 USA case-control 29 1.20 0.50 3.30 0.1823 0.2317
13 13 Akiba 1986 Japan case-control 94 1.52 0.87 2.63 0.4187 0.0796
14 14 Lee 1986 UK case-control 32 1.03 0.41 2.55 0.0296 0.2174
15 15 Koo 1987 Hong Kong case-control 86 1.55 0.90 2.67 0.4383 0.0770
16 16 Pershagen 1987 Sweden case-control 70 1.03 0.61 1.74 0.0296 0.0715
17 17 Humble 1987 USA case-control 20 2.34 0.81 6.75 0.8502 0.2926
18 18 Lam 1987 Hong Kong case-control 199 1.65 1.16 2.35 0.5008 0.0324
19 19 Gao 1987 China case-control 246 1.19 0.82 1.73 0.1740 0.0363
20 20 Brownson 1987 USA case-control 19 1.52 0.39 5.96 0.4187 0.4839
21 21 Geng 1988 China case-control 54 2.16 1.08 4.29 0.7701 0.1238
22 22 Shimizu 1988 Japan case-control 90 1.08 0.64 1.82 0.0770 0.0711
23 23 Inoue 1988 Japan case-control 22 2.55 0.74 8.78 0.9361 0.3982
24 24 Kalandidi 1990 Greece case-control 90 1.62 0.90 2.91 0.4824 0.0896
25 25 Sobue 1990 Japan case-control 144 1.06 0.74 1.52 0.0583 0.0337
26 26 Wu-Williams 1990 China case-control 417 0.79 0.62 1.02 -0.2357 0.0161
27 27 Liu 1991 China case-control 54 0.74 0.32 1.69 -0.3011 0.1802
28 28 Jockel 1991 Germany case-control 23 2.27 0.75 6.82 0.8198 0.3171
29 29 Brownson 1992 USA case-control 431 0.97 0.78 1.21 -0.0305 0.0125
30 30 Stockwell 1992 USA case-control 210 1.60 0.80 3.00 0.4700 0.1137
31 31 Du 1993 China case-control 75 1.19 0.66 2.13 0.1740 0.0893
32 32 Liu 1993 China case-control 38 1.66 0.73 3.78 0.5068 0.1760
33 33 Fontham 1994 USA case-control 651 1.26 1.04 1.54 0.2311 0.0100
34 34 Kabat 1995 USA case-control 67 1.10 0.62 1.96 0.0953 0.0862
35 35 Zaridze 1995 Russia case-control 162 1.66 1.12 2.45 0.5068 0.0399
36 36 Sun 1996 China case-control 230 1.16 0.80 1.69 0.1484 0.0364
37 37 Wang 1996 China case-control 135 1.11 0.67 1.84 0.1044 0.0664 Meta-analytic model
hackshaw1998_or <- rma(yi, vi, data = dat.hackshaw1998, measure = "OR", method = "FE",
slab = paste(author, year, sep = ", "))
hackshaw1998_or
Fixed-Effects Model (k = 37)
I^2 (total heterogeneity / total variability): 24.21%
H^2 (total variability / sampling variability): 1.32
Test for Heterogeneity:
Q(df = 36) = 47.4979, p-val = 0.0952
Model Results:
estimate se zval pval ci.lb ci.ub
0.1858 0.0373 4.9797 <.0001 0.1126 0.2589 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Funnel plots
par(mfrow = c(2, 2))
funnel(hackshaw1998_or, main = "Standard Error")
funnel(hackshaw1998_or, yaxis = "vi", main = "Sampling Variance")
funnel(hackshaw1998_or, yaxis = "seinv", main = "Inverse Standard Error")
funnel(hackshaw1998_or, yaxis = "vinv", main = "Inverse Sampling Variance")
Egger test
regtest(hackshaw1998_or)
Regression Test for Funnel Plot Asymmetry
Model: fixed-effects meta-regression model
Predictor: standard error
Test for Funnel Plot Asymmetry: z = 2.5745, p = 0.0100
Limit Estimate (as sei -> 0): b = 0.0050 (CI: -0.1508, 0.1608)Trim-and-fill analysis
taf <- trimfill(hackshaw1998_or)
taf
Estimated number of missing studies on the left side: 7 (SE = 4.0405)
Fixed-Effects Model (k = 44)
I^2 (total heterogeneity / total variability): 30.44%
H^2 (total variability / sampling variability): 1.44
Test for Heterogeneity:
Q(df = 43) = 61.8207, p-val = 0.0313
Model Results:
estimate se zval pval ci.lb ci.ub
0.1556 0.0364 4.2706 <.0001 0.0842 0.2270 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1funnel(taf, legend = TRUE)
Cumulative meta-analysis
In the order of publication year.
tmp <- cumul(hackshaw1998_or, order = order(dat.hackshaw1998$year))
tmp
estimate se zval pvals ci.lb ci.ub Q Qp I2
Garfinkel, 1981 0.1655 0.1370 1.2079 0.2271 -0.1031 0.4341 0.0000 1.0000 0.0000
Chan, 1982 0.0791 0.1233 0.6414 0.5213 -0.1625 0.3207 2.0866 0.1486 52.0746
Gao, 1987 0.1071 0.1035 1.0347 0.3008 -0.0957 0.3099 2.2615 0.3228 11.5640
Wang, 1996 0.1067 0.0960 1.1111 0.2665 -0.0815 0.2949 2.2616 0.5199 0.0000
Hirayama, 1984 0.1645 0.0849 1.9372 0.0527 -0.0019 0.3309 3.9214 0.4167 0.0000
Butler, 1988 0.1716 0.0844 2.0342 0.0419 0.0063 0.3369 4.4513 0.4864 0.0000
Cardenas, 1997 0.1736 0.0761 2.2799 0.0226 0.0244 0.3228 4.4544 0.6154 0.0000
Correa, 1983 0.1873 0.0752 2.4920 0.0127 0.0400 0.3347 5.7708 0.5668 0.0000
Trichopolous, 1983 0.2213 0.0729 3.0363 0.0024 0.0785 0.3642 9.1915 0.3264 12.9626
Buffler, 1984 0.2094 0.0719 2.9120 0.0036 0.0685 0.3504 10.1893 0.3354 11.6720
Kabat, 1984 0.2027 0.0714 2.8404 0.0045 0.0628 0.3426 10.7650 0.3761 7.1061
Lam, 1985 0.2272 0.0696 3.2648 0.0011 0.0908 0.3636 13.1441 0.2840 16.3124
Garfinkel, 1985 0.2252 0.0662 3.4046 0.0007 0.0956 0.3549 13.1522 0.3581 8.7605
Wu, 1985 0.2245 0.0655 3.4245 0.0006 0.0960 0.3529 13.1600 0.4355 1.2158
Akiba, 1986 0.2344 0.0638 3.6713 0.0002 0.1093 0.3595 13.6096 0.4792 0.0000
Lee, 1986 0.2306 0.0633 3.6460 0.0003 0.1066 0.3546 13.7990 0.5408 0.0000
Koo, 1987 0.2409 0.0617 3.9060 0.0001 0.1200 0.3618 14.3315 0.5740 0.0000
Pershagen, 1987 0.2302 0.0601 3.8309 0.0001 0.1124 0.3480 14.9246 0.6009 0.0000
Humble, 1987 0.2378 0.0597 3.9810 0.0001 0.1207 0.3548 16.2222 0.5770 0.0000
Lam, 1987 0.2638 0.0567 4.6539 0.0000 0.1527 0.3749 18.1434 0.5129 0.0000
Brownson, 1987 0.2648 0.0565 4.6874 0.0000 0.1541 0.3756 18.1926 0.5747 0.0000
Geng, 1988 0.2775 0.0558 4.9751 0.0000 0.1682 0.3869 20.2026 0.5085 0.0000
Shimizu, 1988 0.2691 0.0546 4.9287 0.0000 0.1621 0.3762 20.7449 0.5365 0.0000
Inoue, 1988 0.2741 0.0544 5.0383 0.0000 0.1675 0.3807 21.8537 0.5291 0.0000
Kalandidi, 1990 0.2808 0.0535 5.2452 0.0000 0.1758 0.3857 22.3225 0.5600 0.0000
Sobue, 1990 0.2633 0.0514 5.1244 0.0000 0.1626 0.3640 23.6755 0.5382 0.0000
Wu-Williams, 1990 0.1931 0.0476 4.0542 0.0001 0.0998 0.2865 36.9441 0.0756 29.6234
Liu, 1991 0.1870 0.0473 3.9498 0.0001 0.0942 0.2798 38.2825 0.0735 29.4717
Jockel, 1991 0.1914 0.0472 4.0578 0.0000 0.0990 0.2839 39.5363 0.0727 29.1789
Brownson, 1992 0.1580 0.0435 3.6342 0.0003 0.0728 0.2432 42.8688 0.0468 32.3518
Stockwell, 1992 0.1631 0.0431 3.7826 0.0002 0.0786 0.2476 43.7111 0.0506 31.3675
Du, 1993 0.1633 0.0427 3.8270 0.0001 0.0797 0.2470 43.7124 0.0646 29.0818
Liu, 1993 0.1668 0.0425 3.9296 0.0001 0.0836 0.2500 44.3760 0.0716 27.8889
Fontham, 1994 0.1766 0.0391 4.5186 0.0000 0.1000 0.2532 44.7251 0.0837 26.2160
Kabat, 1995 0.1752 0.0387 4.5220 0.0000 0.0993 0.2512 44.8005 0.1018 24.1080
Zaridze, 1995 0.1872 0.0380 4.9226 0.0000 0.1127 0.2618 47.4581 0.0779 26.2508
Sun, 1996 0.1858 0.0373 4.9797 0.0000 0.1126 0.2589 47.4979 0.0952 24.2072
H2
Garfinkel, 1981 1.0000
Chan, 1982 2.0866
Gao, 1987 1.1308
Wang, 1996 0.7539
Hirayama, 1984 0.9803
Butler, 1988 0.8903
Cardenas, 1997 0.7424
Correa, 1983 0.8244
Trichopolous, 1983 1.1489
Buffler, 1984 1.1321
Kabat, 1984 1.0765
Lam, 1985 1.1949
Garfinkel, 1985 1.0960
Wu, 1985 1.0123
Akiba, 1986 0.9721
Lee, 1986 0.9199
Koo, 1987 0.8957
Pershagen, 1987 0.8779
Humble, 1987 0.9012
Lam, 1987 0.9549
Brownson, 1987 0.9096
Geng, 1988 0.9620
Shimizu, 1988 0.9429
Inoue, 1988 0.9502
Kalandidi, 1990 0.9301
Sobue, 1990 0.9470
Wu-Williams, 1990 1.4209
Liu, 1991 1.4179
Jockel, 1991 1.4120
Brownson, 1992 1.4782
Stockwell, 1992 1.4570
Du, 1993 1.4101
Liu, 1993 1.3867
Fontham, 1994 1.3553
Kabat, 1995 1.3177
Zaridze, 1995 1.3559
Sun, 1996 1.3194 forest(tmp, at = log(c(0.75, 1, 1.25, 1.5, 1.75)),
atransf = exp, cex = 0.75, header = "Author(s) and Year")
Dose-response meta-analysis
Data
data(coffee_mort)
head(coffee_mort, 10) id author year type dose cases n logrr se gender area
1 1 LeGrady et al. 1987 ci 0.5 57 249 0.0000000 0.0000000 M USA
2 1 LeGrady et al. 1987 ci 2.5 136 655 -0.2876821 0.1391187 M USA
3 1 LeGrady et al. 1987 ci 4.5 144 619 -0.1743534 0.1373198 M USA
4 1 LeGrady et al. 1987 ci 6.5 115 387 0.0861777 0.1401409 M USA
5 2 Rosengren et al. 1991 ci 0.0 17 192 0.0000000 0.0000000 M Europe
6 2 Rosengren et al. 1991 ci 1.5 88 1121 -0.1203576 0.2531537 M Europe
7 2 Rosengren et al. 1991 ci 3.5 155 2464 -0.3418342 0.2442559 M Europe
8 2 Rosengren et al. 1991 ci 5.5 155 1986 -0.1261707 0.2440559 M Europe
9 2 Rosengren et al. 1991 ci 7.5 39 575 -0.2665264 0.2784189 M Europe
10 2 Rosengren et al. 1991 ci 9.5 24 427 -0.5108256 0.2802582 M EuropeLinear and spline model
# linear model
lin <- dosresmeta(logrr ~ dose, id = id, se = se, type = type,
cases = cases, n = n, data = coffee_mort, method = "mm")
summary(lin)Call: dosresmeta(formula = logrr ~ dose, id = id, type = type, cases = cases,
n = n, data = coffee_mort, se = se, method = "mm")
Two-stage random-effects meta-analysis
Estimation method:
Covariance approximation: Greenland & Longnecker
Chi2 model: X2 = 44.4226 (df = 1), p-value = 0.0000
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
(Intercept) -0.0324 0.0049 -6.6650 0.0000 -0.0419 -0.0229 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Std. Dev
0.0163
Univariate Cochran Q-test for residual heterogeneity:
Q = 77.0088 (df = 21), p-value = 0.0000
I-square statistic = 72.7%
22 studies, 22 values, 1 fixed and 1 random-effects parameters# spline model
library(rms)
k <- with(coffee_mort, quantile(dose, c(.1, .5, .9)))
spl <- dosresmeta(logrr ~ rcs(dose, k), id = id, se = se, type = type,
cases = cases, n = n, data = coffee_mort, method = "mm")
summary(spl)Call: dosresmeta(formula = logrr ~ rcs(dose, k), id = id, type = type,
cases = cases, n = n, data = coffee_mort, se = se, method = "mm")
Two-stage random-effects meta-analysis
Estimation method:
Covariance approximation: Greenland & Longnecker
Chi2 model: X2 = 209.5069 (df = 2), p-value = 0.0000
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
rcs(dose, k)dose -0.0764 0.0105 -7.3101 0.0000 -0.0969 -0.0559 ***
rcs(dose, k)dose' 0.1053 0.0229 4.5935 0.0000 0.0604 0.1502 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Std. Dev Corr
rcs(dose, k)dose 0.0354 rcs(dose, k)dose
rcs(dose, k)dose' 0.0797 -1
Univariate Cochran Q-test for residual heterogeneity:
Q = 101.3912 (df = 42), p-value = 0.0000
I-square statistic = 58.6%
22 studies, 44 values, 2 fixed and 1 random-effects parametersGraphical prediction
newdata <- data.frame(dose = seq(0, 7, .2)) %>%
cbind(spl = predict(spl, newdata = ., xref = 0, expo = TRUE)) %>%
cbind(lin = predict(lin, newdata = ., xref = 0, expo = TRUE))
ggplot(newdata, aes(dose, spl.pred)) +
geom_line() +
geom_ribbon(aes(ymin = spl.ci.lb, ymax = spl.ci.ub), alpha = .1) +
geom_line(aes(y = lin.pred), linetype = "dashed") +
scale_y_continuous(trans = "log") +
labs(x = "Coffee consumption (cups/day)", y = "Relative Risk") +
theme_classic()