R语言 | metafor的Meta分析及森林图示例
metafor程序包是一个用于Meta分析的优秀R包，整合了非常多的Meta分析方法，也提供了优秀的森林图可视化方案。以上示例图就是用metafor包获得的。这两天在学习metafor包的森林图方法时，偶然间发现了好东西，首先给同学们分享下。
链接：http://www.metafor-project.org/doku.php/analyses
该界面中，作者收集了多篇文献中使用到的Meta分析模型、方法和技术，并提供了分析代码，以及将原始数据包装在metafor包中，真的是太赞了。我们可以很方便地从中学习、模拟这些文献中作者的分析思路。
这些Meta分析的细节部分，同学们有兴趣可自行学习，作者的文档写的非常详细。
然后本篇展示使用metafor包执行某种Meta分析并绘制森林图的一个简单示例。森林图的风格有很多种，下文搬运文档中的一个带亚组的森林图：
http://www.metafor-project.org/doku.php/plots:forest_plot_with_subgroups
使用metafor包的示例数据集“dat.bcg”作展示，该数据集记录了13项研究结果证明卡介苗芽孢杆菌（BCG）疫苗对预防结核病的有效性。
library(metafor)#示例数据集，详情 ?dat.bcgdata(dat.bcg)head(dat.bcg)
其中：
tpos，接种疫苗组的结核病阳性病例数；
tneg，接种疫苗组的结核病阴性病例数；
cpos，对照组的结核病阳性病例数；
cneg，对照组的结核病阴性病例数。
现在使用随机效应模型，评估接种疫苗是否有助于人群提升对疾病的免疫（或者说降低风险）。
#拟合随机效应模型（random-effects model），详情 ?rmares <- rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg,measure = 'RR', slab = paste(author, year, sep = ', '), method = 'REML')res#查看结果中包含的各项统计量概要，并可从中提取想要的信息#详见 ?rma 中对各项统计量的描述names(res)res$pval #如 p 值
Meta分析结果的细节部分参阅帮助文档即可。
接下来使用森林图展示上述回归中拟合的各变量的效应、显著性以及置信区间等。
#绘制森林图#可自定义设置坐标轴刻度、文本范围、字体大小等，详情 ?forestforest(res, xlim = c(-16, 6), at = log(c(0.05, 0.25, 1, 4)), atransf = exp,ilab = cbind(dat.bcg$tpos, dat.bcg$tneg, dat.bcg$cpos, dat.bcg$cneg),ilab.xpos = c(-9.5, -8, -6, -4.5), cex = 0.75, ylim = c(-1, 27),order = order(dat.bcg$alloc), rows = c(3:4, 9:15, 20:23),xlab = 'Risk Ratio', mlab = '', psize = 1, header = 'Author(s) and Year')#添加 Q 值、自由度、p 值等统计量text(-16, -1, pos = 4, cex = 0.75, bquote(paste('RE Model for All Studies (Q = ',.(formatC(res$QE, digits = 2, format = 'f')), ', df = ', .(res$k - res$p),', p = ', .(formatC(res$QEp, digits = 2, format = 'f')), '; ', I^2, ' = ',.(formatC(res$I2, digits = 1, format = 'f')), '%)')))
#继续为子组添加文本op <- par(cex = 0.75, font = 4) #定义字体尺寸text(-16, c(24,16,5), pos = 4, c('Systematic Allocation', 'Random Allocation', 'Alternate Allocation'))#将列标题添加到绘图中par(font = 2) #定义字体加粗text(c(-9.5, -8, -6, -4.5), 26, c('TB+', 'TB-', 'TB+', 'TB-'))text(c(-8.75, -5.25), 27, c('Vaccinated', 'Control'))
#继续拟合三个亚组的随机效应模型（random-effects model），并将结果添加到森林图中res.s <- rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg, measure = 'RR',subset = (alloc == 'systematic'), method = 'REML')res.r <- rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg, measure = 'RR',subset = (alloc == 'random'), method = 'REML')res.a <- rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg, measure = 'RR',subset = (alloc == 'alternate'), method = 'REML')#将计算结果继续添加到上述森林图中par(op) #返回默认字体样式addpoly(res.s, row = 18.5, cex = 0.75, atransf = exp, mlab = '')addpoly(res.r, row = 7.5, cex = 0.75, atransf = exp, mlab = '')addpoly(res.a, row = 1.5, cex = 0.75, atransf = exp, mlab = '')#添加 Q 值、自由度、p 值等统计量text(-16, 18.5, pos = 4, cex = 0.75, bquote(paste('RE Model for Subgroup (Q = ',.(formatC(res.s$QE, digits = 2, format = 'f')), ', df = ', .(res.s$k - res.s$p),', p = ', .(formatC(res.s$QEp, digits = 2, format = 'f')), '; ', I^2, ' = ',.(formatC(res.s$I2, digits = 1, format = 'f')), '%)')))text(-16, 7.5, pos = 4, cex = 0.75, bquote(paste('RE Model for Subgroup (Q = ',.(formatC(res.r$QE, digits = 2, format = 'f')), ', df = ', .(res.r$k - res.r$p),', p = ', .(formatC(res.r$QEp, digits = 2, format = 'f')), '; ', I^2, ' = ',.(formatC(res.r$I2, digits = 1, format = 'f')), '%)')))text(-16, 1.5, pos = 4, cex = 0.75, bquote(paste('RE Model for Subgroup (Q = ',.(formatC(res.a$QE, digits = 2, format = 'f')), ', df = ', .(res.a$k - res.a$p),', p = ', .(formatC(res.a$QEp, digits = 2, format = 'f')), '; ', I^2, ' = ',.(formatC(res.a$I2, digits = 1, format = 'f')), '%)')))