R语言 | metafor的Meta分析及森林图示例

Posted 2021-04-18 生物空间站

R包metafor的Meta分析及森林图示例

森林图（ forest plots ）是 Meta 分析中最常用的结果综合表达形式，它以一条垂直于 X 轴的无效线（通常坐标 X=1 或 0 ）为中心，使用平行于横轴的多条线段表示了每个被纳入研究的效应量和可信区间。通过森林图，我们可以很直观地从中获得变量的综合效应。

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')), '%)')))

我就知道你在看