Qin Liu et al.(2009), A Two-Stage Hierarchical Regression Model for Meta-Analysis of Epidemiologic Nonlinear Dose-Response Data Dos-Response Data
Qin Liu et al. (2009) describe how to perform a two-stage analysis to estimate estimate the summarized nonlinear dose-response relation.
The aim of the first-stage analysis is to estimate in each study the (same) dose-response association between the adjusted log relative risks and the levels of exposure, as described previously by Greenland and Longnecker (1992). Practical examples can be found here. In the second stage analysis the study-specific estimates are combined using extension of the generalized least square method and of the multivariate maximum likelihood approaches, as described by Berkey et al. (1998).
The authors describe in details how to construct the design matrix. As the dose-specific relative risks are estimated as contrasts to their reference exposure, also the design matrices need to be constructed in a similar way. In the dosresmeta function this is done internally by the default option center = TRUE. The argument is particularly important when the reference exposure levels varies across studies and/or for non-zero reference exposures. In addition, the dose-response model typically does not include the intercept, since the log relative risk is 0 by definition for the referent value (Liu et al. 2009).
The methodology is illustrated reanalyzing two meta-analyses, where the non-linearity is investigated adopting quadratic models.
Example 1. Studies with zero exposure dose as reference
The first example includes six studies investigating the dose-response relationship between alcohol intake and vascular disease. We present the code to reproduce main tables and figures.
require("dosresmeta")
alcohol_cvd <- read.table("http://alessiocrippa.altervista.org/data/alcohol_cvd.txt")
alcohol_cvd id author type dose cases control n logrr se
1 1 Bianchi cc 0.000 126 288 414 0.0000000 NA
2 1 Bianchi cc 9.060 61 200 261 -0.2231435 0.2233380
3 1 Bianchi cc 27.000 69 159 228 -0.0001000 0.2337519
4 1 Bianchi cc 45.000 22 22 44 0.5306283 0.3765137
5 1 Bianchi cc 64.800 19 15 34 0.8754687 0.4440046
6 2 Bobak cc 0.000 77 181 258 0.0000000 NA
7 2 Bobak cc 16.050 88 325 413 -0.4307829 0.2213052
8 2 Bobak cc 46.425 24 178 202 -1.0788100 0.2975654
9 2 Bobak cc 77.160 13 51 64 -0.6161861 0.3870792
10 3 Malarcher-wine cc 0.000 83 125 208 0.0000000 NA
11 3 Malarcher-wine cc 1.176 46 129 175 -0.5798185 0.3171470
12 3 Malarcher-wine cc 8.916 17 41 58 -0.5621189 0.4186140
13 3 Malarcher-wine cc 18.720 4 7 11 0.6151856 0.9091006
14 4 Malarcher-beer cc 0.000 83 125 208 0.0000000 NA
15 4 Malarcher-beer cc 0.955 29 88 117 -0.2876821 0.3332336
16 4 Malarcher-beer cc 7.430 32 49 81 0.5128236 0.3383726
17 4 Malarcher-beer cc 15.600 18 21 39 -0.3147107 0.4713454
18 5 Vliegenthart-wine ci 0.000 159 321 480 0.0000000 NA
19 5 Vliegenthart-wine ci 6.000 229 746 975 -0.4155155 0.1442153
20 5 Vliegenthart-wine ci 18.000 38 169 207 -0.6348783 0.2383486
21 5 Vliegenthart-wine ci 28.800 39 94 133 -0.0833816 0.2738854
22 6 Vliegenthart-beer ci 0.000 302 967 1269 0.0000000 NA
23 6 Vliegenthart-beer ci 6.250 129 307 436 -0.2231435 0.1565872
24 6 Vliegenthart-beer ci 18.750 19 30 49 -0.0100503 0.3510891
25 6 Vliegenthart-beer ci 30.000 15 26 41 -0.2613648 0.4165114Fig. 1. Summarized data on the relation between alcohol intake and vascular disease with estimated trend of ln RR according to the amount of alcohol intake.
require("ggplot2")
ggplot(alcohol_cvd, aes(dose, logrr, group = id, shape = author)) +
geom_line() + geom_point() + theme_classic() + ylab("LnRR") +
xlab("Alcohol intake (mL/day)") + xlim(0, 90) + ylim(-1.5, 1.5)
The authors consider two models: a) nonlinear dose-response relation that includes both the linear and quadratic term of alcohol intake; b) as in a) but including the type of study as a covariate.
N.B. meta regression models are not yet implemented alone in the dosresmeta package, but can be achieved using the mvmeta package.
## model a)
quadrAlc <- dosresmeta(formula = logrr ~ dose + I(dose^2), id = id,
type = type, se = se, cases = cases, n = n,
data = alcohol_cvd, method = "mm")
summary(quadrAlc)Call: dosresmeta(formula = logrr ~ dose + I(dose^2), id = id, type = type,
cases = cases, n = n, data = alcohol_cvd, se = se, method = "mm")
Two-stage random-effects meta-analysis
Estimation method:
Covariance approximation: Greenland & Longnecker
Chi2 model: X2 = 10.2626 (df = 2), p-value = 0.0059
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
dose.(Intercept) -0.0307 0.0121 -2.5298 0.0114 -0.0545 -0.0069 *
I(dose^2).(Intercept) 0.0007 0.0005 1.5604 0.1187 -0.0002 0.0017
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Std. Dev Corr
dose 0.0175 dose
I(dose^2) 0.0008 -0.9883
(Note: Truncated estimate - 1 negative eigenvalues set to 0)
Univariate Cochran Q-test for residual heterogeneity:
Q = 35.1177 (df = 10), p-value = 0.0001
I-square statistic = 71.5%
6 studies, 12 values, 2 fixed and 1 random-effects parameters## model b)
type <- as.factor(with(alcohol_cvd,
tapply(type, id, function(x) unique(x))))
metaMod <- mvmeta(quadrAlc$bi ~ type, S = quadrAlc$Si, method = "mm")
summary(metaMod)Call: mvmeta(formula = quadrAlc$bi ~ type, S = quadrAlc$Si, method = "mm")
Multivariate random-effects meta-regression
Dimension: 2
Estimation method: Method of moments
Fixed-effects coefficients
dose :
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
(Intercept) -0.0173 0.0157 -1.1051 0.2691 -0.0481 0.0134
type2 -0.0439 0.0283 -1.5507 0.1210 -0.0994 0.0116
I(dose^2) :
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
(Intercept) 0.0002 0.0005 0.4461 0.6555 -0.0008 0.0013
type2 0.0018 0.0010 1.7161 0.0861 -0.0003 0.0038 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Structure: General positive-definite
Std. Dev Corr
dose 0.0202 dose
I(dose^2) 0.0008 -0.9872
(Note: Truncated estimate - 1 negative eigenvalues set to 0)
Multivariate Cochran Q-test for residual heterogeneity:
Q = 27.5794 (df = 8), p-value = 0.0006
I-square statistic = 71.0%
6 studies, 12 observations, 4 fixed and 1 random-effects parametersTable 3. Parameter estimates using data on alcohol and vascular disease based on Model a).
pi <- do.call("rbind", mapply(function(b, s)
2 * (1 - pnorm(abs(b/(diag(s)^.5)))),
b = with(quadrAlc, split(bi, seq(nrow(bi)))),
s = quadrAlc$Si, SIMPLIFY = F))
require("Gmisc")
tab3 <- cbind(matrix(paste0(quadrAlc$bi, " (", round(pi, 2), ")"),
nrow(quadrAlc$bi)),
do.call("rbind", lapply(quadrAlc$Si, function(x)
round(x[!lower.tri(x)], 8))))
rownames(tab3) <- unique(alcohol_cvd$author)
htmlTable(tab3, align = 'llccc', header = c("b1i (p-value)", "b2i (p-value)", "Si11","Si12", "Si22"), rowlabel = "Author")| Author | b1i (p-value) | b2i (p-value) | Si11 | Si12 | Si22 |
|---|---|---|---|---|---|
| Bianchi | -0.00977694352242094 (0.48) | 0.00040324089054791 (0.11) | 0.00019568 | -3.32e-06 | 6e-08 |
| Bobak | -0.0395170268874542 (0) | 0.000389663934540105 (0.04) | 0.00017888 | -2.35e-06 | 3e-08 |
| Malarcher-wine | -0.168084963150333 (0.08) | 0.0116306998871314 (0.07) | 0.00902574 | -0.0005647 | 4.085e-05 |
| Malarcher-beer | 0.152523788983611 (0.07) | -0.0104395378224687 (0.08) | 0.00732659 | -0.00048898 | 3.61e-05 |
| Vliegenthart-wine | -0.0876806755214546 (0) | 0.00293590183557997 (0) | 0.00061985 | -2.106e-05 | 8.1e-07 |
| Vliegenthart-beer | -0.029984443070979 (0.29) | 0.000853856264988567 (0.45) | 0.00080749 | -2.97e-05 | 1.27e-06 |
It is not included in the original article, but we think it is important to present both a graphical and tabular presentation of the pooled results.
newdata <- data.frame(dose = seq(0, 60, .5))
prediction <- predict(quadrAlc, newdata, expo = T)
## Trick to "fix" the "bug" in geom_ribbon
prediction$ci.ub[prediction$ci.ub>=8] <- 8
prediction$ci.lb[prediction$ci.lb<=.3] <- .3
ggplot(prediction,
aes(dose, pred)) + geom_line() + theme_classic() +
geom_ribbon(aes(ymin = ci.lb, ymax = ci.ub), alpha = .05) +
scale_y_continuous("Relative risk", trans="log",
breaks = c(.5, 1, 2, 3, 4, 7), limits = (c(.3, 8))) +
scale_x_continuous("Alcohol intake, grams/day", breaks = seq(0, 60, 12))
newdata <- data.frame(dose = seq(0, 60, 12))
predict(quadrAlc, newdata, expo = T) dose I(dose^2) pred ci.lb ci.ub
1 0 0 1.0000000 1.0000000 1.0000000
2 12 144 0.7698940 0.6495077 0.9125938
3 24 576 0.7346081 0.5873315 0.9188151
4 36 1296 0.8687090 0.5060480 1.4912723
5 48 2304 1.2731712 0.3868184 4.1905056
6 60 3600 2.3125607 0.2751283 19.4379728Example 2. Studies with non-zero exposure dose as reference
The second example include 15 studies published between 1966 and 1998 investigated the relation between obesity and renal cell cancer. We base our analysis only on the data presented in table 6 of the original paper, so the results will differ from those presented by the authors.
require("dosresmeta")
bmi_rc <- read.table("http://alessiocrippa.altervista.org/data/bmi_rc.txt")
bmi_rc id author type interval bmi cases control n logrr se
1 1 McLaughlin 1984 M cc <=23.6 21.60 80 102 182 0.0000000 0.0000000
2 1 McLaughlin 1984 M cc 23.7-25.6 24.60 66 115 181 -0.1053605 0.2161514
3 1 McLaughlin 1984 M cc 25.6-28 26.80 74 112 186 -0.1053605 0.2161514
4 1 McLaughlin 1984 M cc 28-32.9 30.50 90 97 187 0.4054651 0.2233380
5 2 McLaughlin 1984 F cc <=21.6 19.70 35 74 109 0.0000000 0.0000000
6 2 McLaughlin 1984 F cc 21.7-23.5 22.55 37 72 109 0.0953102 0.3195882
7 2 McLaughlin 1984 F cc 23.5-26.2 24.85 48 65 113 0.4054651 0.3195882
8 2 McLaughlin 1984 F cc 26.2-31.7 29.00 58 54 112 0.7419373 0.3006828
9 3 Goodman 1986 M cc <24 20.00 28 49 77 0.0000000 0.0000000
10 3 Goodman 1986 M cc 24-27 25.50 80 77 157 0.6418539 0.3162535
11 3 Goodman 1986 M cc 27-34 31.00 65 47 112 0.9932518 0.3493654
12 4 Goodman 1986 F cc <24 20.00 30 37 67 0.0000000 0.0000000
13 4 Goodman 1986 F cc 24-27 25.50 13 20 33 -0.2231435 0.4570899
14 4 Goodman 1986 F cc 27-34 31.00 28 14 42 0.8754688 0.4462326
15 5 Asal 1988 M cc <23.8 22.00 40 55 95 0.0000000 0.0000000
16 5 Asal 1988 M cc 23.8-25.6 24.70 46 51 97 0.2623642 0.3034709
17 5 Asal 1988 M cc 25.6-29.4 27.50 50 60 110 0.1823216 0.2802634
18 5 Asal 1988 M cc 29.4-37 33.20 73 29 102 1.1939224 0.3113583
19 6 Asal 1988 F cc <22.5 19.30 24 36 60 0.0000000 0.0000000
20 6 Asal 1988 F cc 22.5-25.7 24.10 24 40 64 -0.2231435 0.3536530
21 6 Asal 1988 F cc 25.7-30.1 27.90 24 32 56 0.1823216 0.3640670
22 6 Asal 1988 F cc 30.1-38.9 34.50 28 30 58 0.1823216 0.3740724
23 7 Yuan 1998 M cc <22 20.00 62 112 174 0.0000000 0.0000000
24 7 Yuan 1998 M cc 22-23 22.50 152 179 331 0.5306283 0.2094377
25 7 Yuan 1998 M cc 24-25 24.50 200 229 429 0.4700036 0.1990237
26 7 Yuan 1998 M cc 26-27 26.50 129 128 257 0.6931472 0.2216974
27 7 Yuan 1998 M cc 28-29 28.50 107 79 186 0.9932518 0.2367357
28 7 Yuan 1998 M cc 30-32 31.00 131 57 188 1.5260563 0.2424004
29 8 Yuan 1998 F cc <22 20.00 126 177 303 0.0000000 0.0000000
30 8 Yuan 1998 F cc 22-23 22.50 95 90 185 0.5306283 0.2094377
31 8 Yuan 1998 F cc 24-25 24.50 66 68 134 0.4054651 0.2124807
32 8 Yuan 1998 F cc 26-27 26.50 35 42 77 0.2623642 0.2921310
33 8 Yuan 1998 F cc 28-29 28.50 32 20 52 0.8329091 0.3195882Fig. 2. Summarized data on the relation between BMI and renal cell cancer with estimated trend of ln RR according to the level of BMI.
ggplot(bmi_rc, aes(bmi, logrr, group = author, shape = author)) +
scale_shape_manual(values = seq_along(unique(bmi_rc$id))) +
labs(shape = unique(bmi_rc$author)) +
geom_line() + geom_point() + theme_classic() + ylab("LnRR") +
xlab("BMI") + xlim(15, 35) + ylim(-.5, 2)
quadrBmi <- dosresmeta(formula = logrr ~ bmi + I(bmi^2), id = id,
type = type, se = se, cases = cases, n = n,
data = bmi_rc, method = "mm")
summary(quadrBmi)Call: dosresmeta(formula = logrr ~ bmi + I(bmi^2), id = id, type = type,
cases = cases, n = n, data = bmi_rc, se = se, method = "mm")
Two-stage random-effects meta-analysis
Estimation method:
Covariance approximation: Greenland & Longnecker
Chi2 model: X2 = 28.6197 (df = 2), p-value = 0.0000
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
bmi.(Intercept) -0.1094 0.1375 -0.7951 0.4265 -0.3789 0.1602
I(bmi^2).(Intercept) 0.0036 0.0027 1.3289 0.1839 -0.0017 0.0089
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Std. Dev Corr
bmi 0.1612 bmi
I(bmi^2) 0.0033 -0.9865
Univariate Cochran Q-test for residual heterogeneity:
Q = 23.4445 (df = 14), p-value = 0.0534
I-square statistic = 40.3%
8 studies, 16 values, 2 fixed and 1 random-effects parametersnewdata <- data.frame(bmi = seq(18, 35, .5))
ggplot(predict(quadrBmi, newdata, expo = T),
aes(bmi, pred)) + geom_line() + theme_classic() +
scale_y_continuous("Relative risk" , trans="log",
breaks = c(.5, 1, 2, 3, 4, 6)) +
geom_ribbon(aes(ymin = ci.lb, ymax = ci.ub), alpha = .05)
newdata <- data.frame(bmi = seq(20, 35, 5))
predict(quadrBmi, newdata, expo = T) bmi I(bmi^2) pred ci.lb ci.ub
1 20 400 1.000000 1.000000 1.000000
2 25 625 1.295536 1.055363 1.590366
3 30 900 2.007518 1.521618 2.648581
4 35 1225 3.720741 2.189746 6.322153Information about R session:
R version 3.2.1 (2015-06-18)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows 7 (build 7600)
locale:
[1] LC_COLLATE=Swedish_Sweden.1252 LC_CTYPE=Swedish_Sweden.1252 LC_MONETARY=Swedish_Sweden.1252
[4] LC_NUMERIC=C LC_TIME=Swedish_Sweden.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggplot2_1.0.1 dosresmeta_2.0.0 mvmeta_0.4.5 Gmisc_1.1 htmlTable_1.3 Rcpp_0.11.5
[7] kfigr_1.1.0 knitr_1.10.5
loaded via a namespace (and not attached):
[1] Formula_1.2-1 cluster_2.0.1 magrittr_1.5 MASS_7.3-40 splines_3.2.1
[6] munsell_0.4.2 colorspace_1.2-6 lattice_0.20-31 stringr_0.6.2 plyr_1.8.2
[11] tools_3.2.1 nnet_7.3-9 grid_3.2.1 gtable_0.1.2 latticeExtra_0.6-26
[16] htmltools_0.2.6 forestplot_1.1 abind_1.4-3 yaml_2.1.13 survival_2.38-1
[21] digest_0.6.8 RColorBrewer_1.1-2 reshape2_1.4.1 formatR_1.2 acepack_1.3-3.3
[26] rpart_4.1-9 evaluate_0.7 rmarkdown_0.6.1 labeling_0.3 scales_0.2.4
[31] Hmisc_3.15-0 foreign_0.8-63 proto_0.3-10 References
Berkey, C. S., D. C. Hoaglin, A. Antczak-Bouckoms, F. Mosteller, and G. A. Colditz. 1998. “Meta-Analysis of Multiple Outcomes by Regression with Random Effects.” Statistics in Medicine 17 (22): 2537–50.
Greenland, Sander, and Matthew P. Longnecker. 1992. “Methods for Trend Estimation from Summarized Dose-Response Data, with Applications to Meta-Analysis.” American Journal of Epidemiology 135 (11): 1301–9. http://aje.oxfordjournals.org/content/135/11/1301.
Liu, Qin, Nancy R. Cook, Anna Bergström, and Chung-Cheng Hsieh. 2009. “A Two-Stage Hierarchical Regression Model for Meta-Analysis of Epidemiologic Nonlinear Dose–response Data.” Computational Statistics & Data Analysis 53 (12): 4157–67. doi:10.1016/j.csda.2009.05.001.