BMJ子刊：手把手教你用R语言实现meta分析~

Meta分析是循证医学中非常重要的工具，如果配上一个好的科学问题，它可以整合当前已有的证据，以量化的方式给出相对公正的结论，临床价值也可以非常巨大！

那如何使用R做Meta分析？

是不是觉得无从下手？ 

幸运的是，小编前几天阅读到一篇文献[1]，名为《How to perform a meta-analysis with R: a practical tutorial》，内容实用并且对初学者非常友好，对于想要进入这个领域的朋友可能会有一些帮助。

马上进入操作部分： 

---

相关R包的安装和载入

# install.packages('meta')
# install.packages('metasens')

library(meta)
library(metasens)

R包get！可能是一个伟大的开始！

---

数据的准备

此篇论文使用到的数据来自于2006年Joy等人发表在《Cochrane Database Syst Rev》上的研究，名为“Haloperidol versus placebo for schizophrenia”。此项研究一共纳入了17项临床试验，旨在探索氟哌啶醇（Haloperidol）在精神分裂症患者中的疗效。结局的评价为是否发生了临床症状的缓解（有vs无），因此是一个二分类（Binary）的结局。

其中，风险比（Risk ratio, RR）用于评价药物干预与结局的关系。如果RR大于1，说明氟哌啶醇要优于安慰剂组。

下面用R创建论文中的数据，并且命名为joy： 

joy <- data.frame(author = c('Arvanitis', 'Beasley', 'Bechelli', 'Borison', 'Chouinard', 'Durost', 
                             'Garry', 'Howard', 'Marder', 'Nishikawa', 'Nishikawa', 'Reschke', 'Selman',
                             'Serafetinides', 'Simpson', 'Spencer', 'Vichaiya'), 
                  year = c(1997, 1996, 1983, 1992, 1993, 1964, 1962, 1974, 1994, 1982, 1984, 1974, 1976, 1972, 1967, 1992, 1971), 
                  resp.h = c(25, 29, 12, 3, 10, 11, 7, 8, 19, 1, 11, 20, 17, 4, 2, 11, 9), 
                  fail.h = c(25, 18, 17, 9, 11, 8, 18, 9, 45, 9, 23, 9, 1, 10, 14, 1, 20), 
                  drop.h = c(2, 22, 1, 0, 0, 0, 1, 0, 2, 0, 3, 0, 11, 0, 0, 0, 1), 
                  resp.p = c(18, 20, 2, 0, 3, 1, 4, 3, 14, 0, 0, 2, 7, 0, 0, 1, 0), 
                  fail.p = c(33, 14, 28, 12, 19, 14, 21, 10, 50, 10, 13, 9, 4, 13, 7, 11, 29), 
                  drop.p =c(0, 34, 1, 0, 0, 0, 1, 0, 2, 0, 0, 0, 18, 1, 1, 0, 1))
joy

我们一共创建了17个临床试验的数据（行），并且含有8个变量（列），其中8个变量的解释如下： 

author: 论文的第一作者；

year: 发表的年份；

resp.h: 氟哌啶醇组中有症状改善的人数；

resp.p: 安慰剂组中有症状改善的人数；

fail.h: 氟哌啶醇组中没有改善的人数；

fail.p: 安慰剂组中没有改善的人数；

drop.h: 氟哌啶醇组中脱落的人数；

drop.p: 安慰剂组中脱落的人数。

因为后续需要进行亚组分析，所以事先创建一个新的变量miss，用于标注是否存在缺失值： 

joy$miss <- ifelse((joy$drop.h + joy$drop.p) == 0,
                  c('Without missing data'), c('With missing data'))joy


如果能重复出上图，恭喜你，数据准备这一步就搞定了！

---

固定和随机效应Meta分析

因为此项研究的结局是一个二分类变量，所以将会选择metabin()函数进行分析。而关于结局为连续变量的情况，有机会的话会在未来介绍。

结局为二分类的Meta分析代码如下：

settings.meta(digits = 2) # 显示小数点后两位

m.publ <- metabin(resp.h, resp.h + fail.h, resp.p, resp.p + fail.p, # 最重要的四个输入内容
                 data = joy, 
                 studlab = paste0(author, ' (', year, ')'), # 修饰研究的显示标签：author(year)
                 method.tau = 'PM')    # 选择“Paule and Mandel”，二分类结局的推荐方法

m.publ

如果可以重复出上图，说明meta分析已经初步完成了！

上述结果中含有非常多信息，包括每项研究的结果，RR，可信区间，固定效应，随机效应，异质性，以及meta分析方法的详细信息等。

为了更加直观的显示结果，制作一张森林图（Forest plot），并以PDF的格式保存：

pdf('figure1.pdf', width = 10, height = 6)  # 同时调整图片的宽和高

forest(m.publ,
       sortvar = year,    # 按时间排序
       prediction = TRUE, # 在森林图中显示“prediction interval”
       label.left = 'Favours placebo',      # 在森林图底部的左侧加上标签
       label.right = 'Favours haloperidol') # 在森林图底部的右侧加上标签

dev.off() 

代码运行后，在工作路径中 会出现一个文件: 'figure1.pdf'[如何查找工作路径：如何将Excel或csv文件导入R？]，双击之后可得下图：

从上图可知，在固定效应和随机效应模型中，“钻石” （图片底部的菱形）指代的是RR以及可信区间，并没有“碰到”RR=1的直线，提示氟哌啶醇的疗效显著的优于安慰剂。

但是！

将异质性问题（图片左下角，p < 0.04）也考虑在内的话，故事就不一样了！

因为上图中的prediction interval(红线)碰到了RR=1的直线，提示在未来的研究中，氟哌啶醇可能并不会优于安慰剂。

关于森林图的详细介绍，请询问R：?forest.meta

---

评估缺失数据的影响

为了研究缺失数据对结果的影响，进一步做亚组分析（subgroup analysis）： 

m.publ.sub <- update(m.publ, 
                     byvar = miss,        # 缺失数据的分组变量
                     print.byvar = FALSE) # 不显示标签名字
m.publ.sub

同样，也可以将上述的结果进行作图（森林图），并保存到“figure2.pdf”：

pdf('figure2.pdf', width = 10, height = 7.05)

forest(m.publ.sub,
       sortvar = year,
       xlim = c(0.1, 100),                   # 设定x轴的范围
       at = c(0.1, 0.3, 1, 3, 10, 30, 100))  # x轴上标注具体的刻度

dev.off()

双击“figure2.pdf”可得下图： 

从上图可知，两个亚组中（with missing data vs without missing data）的结果都提示氟哌啶醇要优于安慰剂，但是在without missing data组中，氟哌啶醇的效果更佳（RR值更大）。

关于缺失值的机制到底是什么这个问题，永远是一个迷，但我们可以构建几个关于缺失值的假设（Assumptions）去处理。

其中，函数metamiss()含有一些Imputation methods，有助于处理缺失值的问题，详细了解请通过?metamiss查看，下图是针对结局为分类变量的几种方法： 

下面，将所有不同的缺失值处理方法全部跑一遍，从而可以对比不同方法之间的差别，代码如下： 

# 使用不同的缺失值处理方法 
mmiss.0 <- metamiss(m.publ, drop.h, drop.p)  # 默认的方法为“miss.0”
mmiss.1 <- metamiss(m.publ, drop.h, drop.p, method = '1')  
mmiss.pc <- metamiss(m.publ, drop.h, drop.p, method = 'pc')
mmiss.pe <- metamiss(m.publ, drop.h, drop.p, method = 'pe')
mmiss.p <- metamiss(m.publ, drop.h, drop.p, method = 'p')
mmiss.b <- metamiss(m.publ, drop.h, drop.p, method = 'b', small.values = 'bad')
mmiss.w <- metamiss(m.publ, drop.h, drop.p, method = 'w', small.values = 'bad')
mmiss.gh <- metamiss(m.publ, drop.h, drop.p, method = 'GH')
mmiss.imor2 <- metamiss(m.publ, drop.h, drop.p, method = 'IMOR', IMOR.e = 2)
mmiss.imor0.5 <- metamiss(m.publ, drop.h, drop.p, method = 'IMOR', IMOR.e = 0.5)

# 给各个方法注明标签，用于下方森林图
labels <- c('Available case analysis (ACA)',
          'Impute no events (ICA-0)', 'Impute events (ICA-1)',
          'Observed risk in control group (ICA-pc)',
          'Observed risk in experimental group (ICA-pe)',
          'Observed group-specific risks (ICA-p)',
          'Best-case scenario (ICA-b)', 'Worst-case scenario (ICA-w)',
          'Gamble-Hollis analysis',
          'IMOR.e = 2, IMOR.c = 2', 'IMOR.e = 0.5, IMOR.c = 0.5')

# 使用inverse-variance method进行汇总
m.publ.iv <- update(m.publ, method = 'Inverse')

# 将结果整合
mbr = metabind(m.publ.iv,
               mmiss.0, mmiss.1,
               mmiss.pc, mmiss.pe, mmiss.p,
               mmiss.b, mmiss.w, mmiss.gh,
               mmiss.imor2, mmiss.imor0.5,
               name = labels, pooled = 'random')

# 制作森林图
pdf('figure3.pdf', width = 7, height = 8)

forest(mbr, xlim = c(0.25, 4),
       label.left = 'Favours placebo', 
       label.right = 'Favours haloperidol',
       leftcols = 'studlab', 
       leftlab = 'Meta-Analysis Method',
       type.study = 'diamond',
       hetlab = '', 
       print.Q = TRUE, 
       fs.study = 10)

dev.off()

双击工作路径中的“figure3.pdf”，出图：

从上图可知，大体上来说，不同方法所计算出的RR比较类似，范围从1.90到2.64之间。所有的敏感性分析都提示氟哌啶醇要优于安慰剂，因此，在这项研究中，缺失数据对于结果并没有造成非常严重的问题。

关于Meta分析的缺失值问题，感兴趣的朋友可以参考下方文献深入了解[2]。

小研究效应（samll-study effects）的评估

有时候，相比于大样本研究，一些小样本研究的结果可能会非常不一样（比如，疗效明显优于大样本研究的结果），因此需要评估并且考虑这些结果对meta分析的影响。

可以通过漏斗图（Funnel plot）用于判断是否存在小研究效应，代码如下： 

funnel(m.publ)

从上图可知，横坐标为RR；纵坐标为它的SE（y轴是倒置，下大上小），SE越大则表示越不精确。

如果漏斗图中的散点看上去不对称（上图这个样子），则提示小研究效应的存在。

如果存在发表偏移（Publication bias），漏斗图也会呈现出不对称。

关于如何识别是否存在发表偏移，有一个更加专业的漏斗图可帮助我们判断，名为轮廓增强版漏斗图（Contour-enhanced funnel plot）。

现在就来画一个：

funnel(m.publ, 
       contour.levels = c(0.9, 0.95, 0.99), 
       col.contour = c ('darkgray', 'gray', 'lightgray'))

如上图，发表偏移可能不是导致不对称的主要因素，因为大部分小研究（大SE值的点，即位于漏斗底部的点）分布在白色区域（代表没有统计学意义的区域）。

如果确实存在小研究效应，则需要对其进行校正。

其中，一种比较常见的方法为剪补法（trim-and-fill method），代码如下： 

tf.publ <- trimfill(m.publ)
summary(tf.publ)

如上述结果所示，经过“剪补”，校正后的随机效应量RR = 1.4 （0.83-2.38）， p值为0.2，提示氟哌啶醇并不优于安慰剂。

再作出相应的漏斗图： 

funnel(tf.publ)

另外一种常见校正的方法为“limit meta-analysis”，此方法通过回归进行校正。代码如下： 

l1.publ = limitmeta(m.publ)
summary(l1.publ)

与剪补法的结果类似，校正后的随机效应量RR为1.29(0.93-1.79)， p值为0.12， 提示提示氟哌啶醇并不优于安慰剂。

再作漏斗图： 

funnel(l1.publ)

正如作者所言[1]，这篇文章仅仅只能作为Meta分析的入门读物，想要更加深入的了解，请查看R语言官网的meta分析专栏：

https://cran./view=MetaAnalysis

好啦，今天的内容就到这里。

如果有帮助，记得分享给需要的人！

参考文献

[1]. Balduzzi S, et al. How to perform a meta-analysis with R: a practical tutorial. Evid Based Ment Health 2019. 

[2]. Higgins JPT, White IR, Wood AM. Imputation methods for missing outcome data in meta-analysis of clinical trials. Clin Trials 2008;5:225–39

▌声明：本文由R语言和统计首发

▌编辑：June

▌我们的宗旨是：让R语言和统计变得简单！