【原创】用R软件做剂量反应关系的meta分析(之二)

文章来源：
Orsini N, Ruifeng L, Wolk A, Khudyakov P, Spiegelman D. Meta-analysis for linear and non-linear dose-response relationships: examples, an evaluation of approximations, and software. American Journal of Epidemiology. 2012; 175(1):66-73.
Summarized data about the relation between total alcohol intake and lung cancer risk in 4 prospective cohort studies participating in the Pooling Project of Prospective Studies of Diet and Cancer.
第一步：安装包：dosresmeta mvmeta，加载包、读取数据
install.packages(c("mvmeta","dosresmeta"))
library("dosresmeta")
alcohol_lc <- read.table(" http://alessiocrippa.altervista.org/data/ex_alcohol_lc.txt ")
alcohol_lc

第二步 线性的固定效应模型的结果（这是是用的R软件的一个专门做剂量反应关系模型的命令dosresmeta）先看看线性模型的结果怎么样Fixed-effect dose-response model assuming linearity
lin.fixed <- dosresmeta(formula = logrr ~ dose, type = type, id = id, se = se, cases = cases, n = peryears, data = alcohol_lc, method = "fixed")
summary(lin.fixed)

线性模型拟合是有意义的，
CHI2 model:x2=9.4983,p=0.0021.,logLik=15.7177,AIC=-29.4353,BIC=-30.090
predict(lin.fixed, delta = 12)

增加12个单位的RR值是1.068842
第三步：线性的随机效应模型的情况
## Random-effect dose-response model assuming linearity
lin <- dosresmeta(formula = logrr ~ dose, type = type, id = id, se = se, cases = cases, n = peryears, data = alcohol_lc)
summary(lin)
predict(lin, delta = 12)

和固定效应模型相比，同样是有意义的，且模型的拟合更好，
CHI2 model:x2=4.3753,p=0.00365.logLik=10.8524,AIC=-17.7048,BIC=-19.5076
第四步，我们再来看看，非线性的处理方法吧，这里需要安装和加载RMS软件包和一系列配套的软件包。
## Non-linearity (spline) using random-effect
install.packages("rms")
library("Hmisc")
library("survival")
library("SparseM")
library("rms")
knots <- quantile(alcohol_lc$dose, c(.05, .35, .65, .95))
knots

（这个是设定百分位数点的，可以自己指定，这里设定5%，35%，65%和95%）
spl <- dosresmeta(formula = logrr ~ rcs(dose, knots), type = type, id = id, se = se, cases = cases, n = peryears, data = alcohol_lc)
summary(spl)

模型拟合P=0.1707,无统计学意义。
## Multivariate Wald test
wald.test(b = coef(spl), Sigma = vcov(spl), Terms = 2:3)

模型拟合检验，P=0.2,无统计学意义。
## Tabulate result
pred <- predict(spl, data.frame(dose = seq(0, 60, 12)), xref = 0)
print(pred, digits = 2)

## Figure 1 A of the paper
newdata = data.frame(dose = seq(0, 60, 1))
with(predict(spl, newdata, xref = 0),{
plot(get("rcs(dose, knots)dose"), pred, type = "l", log = "y", ylab = "Relative risk", las = 1,
xlab = "Alcohol intake, grams/day", ylim = c(.8, 1.8), bty = "l")
lines(get("rcs(dose, knots)dose"), ci.lb, lty = "dashed")
lines(get("rcs(dose, knots)dose"), ci.ub, lty = "dashed")
})
rug(alcohol_lc$dose)

注:程序代码来自于:R软件剂量反应包"dosresmeta"的作者Alessio Crippa 特此表示感谢!