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

欢迎关注"R语言和统计"~~


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.r-project.org/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语言和统计变得简单！

< 往期合集 >



2023-02-13 10:44
首页  ꄲ  公众号  ꄲ  BMJ子刊：手把手教你用R语言实现meta分析~