网状meta必备技能之6-利用R中的meta包实现meta分析

其明技术专家 meta分析, 博客, 基础科研 2017年5月3日

2.74W 0 1

我们仅仅是代码的编辑者、整合者、搬运工，仅免费传授方法，下文数据和代码取自于网络和免费软件“R语言说明书”，如果您觉得我们侵犯了您的版权，请通知我们撤稿。请大家谅解，谢谢！

相信大家已经对R软件和Rstudio有了初步的了解，“我们学习的太复杂了，希望代码能简单些，大家能一起飞……”（创始人语录）。

#第一部分感谢“加勒比海带@Shanghai”老师的代码和说明 #

#1 代码来源####

http://3g.dxy.cn/bbs/topic/29474540#!_id=29474558

#第二部分 meta分析的实现#

#2.1总得来说，在R中实现传统的meta分析的包中，个人感觉meta包是比较好的一个包。

在这里将数据读入R等步骤略去，最为新手来说，数据导入最简单的方法是利用R+Rstudio的组合，先利用txt或者CSV格式在你的电脑上将数据整理好，然后利用Rstudio中的import dataset导入。
在这里主要利用meta包中现成的数据进行示例。例子是心梗后服用阿司匹林能否降低死亡率，具体数据如下表所示。其中event.e表示治疗组的死亡人数，n.e表示治疗组的总人数，event.c表示对照组的总人数，n.c表示对照组的总人数，非常简单的一个例子。

study	year	event.e	n.e	event.c	n.c
MRC-1	1974	49	615	67	624
CDP	1976	44	758	64	771
MRC-2	1979	102	832	126	850
GASP	1979	32	317	38	309
PARIS	1980	85	810	52	406
AMIS	1980	246	2267	219	2257
ISIS-2	1988	1570	8587	1720	8600

#2.2 代码

利用meta包中的metabin函数可以顺利地实现效应量的合并，非常简单，代码如下：
library("meta")
data(Fleiss93)
metabin(event.e, n.e, event.c,n.c,data=Fleiss93,sm="OR")
跑出来的结果如下图所示：

第一部分是各个研究的OR值及95%可信区间，以及分别在固定效应模型和随机效应模型下的权重。第二部分是合并的结果，分别列出了固定效应模型和随机效应模型合并的OR值及95%可信区间，以及检验的结果（z值和相应的p值）。第三部分列出了一致性检验的结果。最后一部分是本meta分析所用的具体方法，包括合并效应值的方法和计算研究间方差的方法。
对于该meta分析，首先需要看的是异质性检验的结果，这里研究间方差tau^2=0.0096，I^2=40%左右，经过检验异质性无统计学意义（p=0.127>0.10）,因此可选用固定效应模型进行效应量的合并，结果是OR=0.897, 95%CI: 0.841-0.957,有统计学意义。一般来说，我们只需要看OR值和95%可信区间就可以判断是否有统计学意义，不需要再看相应的p值，OR值的可信区间不包含1等价于p<0.05。比较一下固定效应模型和随机效应模型的结果，我们看到随机效应模型的可信区间更宽，因此结果相对于固定效应模型来说更保守，主要原因是相对于固定效应模型，随机效应模型在计算标准误的时候引入了研究间方差。
接下来的问题是怎么将forest plot（森林图）做出来，meta分析森林图一般来说不能少。代码如下：
metaresult<-metabin(event.e, n.e, event.c,n.c,data=Fleiss93,sm="OR",
studlab=paste(study, year),comb.random=FALSE)
forest(metaresult)
在上段代码中，多了comb.random=FALSE，代表本例不用随机效应模型进行效应量合并，而是采用固定效应模型。同理，如果需要用随机效应模型进行合并只要在后面改成comb.fixed=FALSE即可。在这里我们将meta分析的结果赋给metaresult这个列表，然后用forest函数调用metaresult即可。得到森林图如下图所示。

各研究的OR值和95%可信区间、权重、效应量合并值以及异质性检验结果等都在该图中。
这类研究一般是RCT研究，因此选择RR值可能更好，一般事件发生率比较低的情况下，RR约等于OR值。如果要使用RR值，只需将sm="OR"改成sm="RR"即可。
metaresult<-metabin(event.e, n.e, event.c,n.c,data=Fleiss93,sm="RR",
studlab=paste(study, year),comb.random=FALSE)

到这里基本上将meta分析的主要部分做完了，接下来需要做一个发表偏倚检验和敏感性分析。
发表偏倚的检验由于纳入的研究个数小于10个，在cochrane手册中不建议做test，只要求做一个漏斗图。当研究个数大于等于10各的时候，对于这种四格表资料不建议使用egger检验或begg检验（具体可见我的另外一个陈年老贴），可以使用Peters检验，power稍高（所有的发表偏倚检验方法在研究个数不多的条件下，power都不高），且不需要AS类方法那样做反正弦变换。发表偏倚代码如下：
funnel(metaresult)
metabias(metaresult, method.bias="peters")
如果纳入研究个数小于10例，提示如下：
In print.metabias(list(k = 7L, k.min = 10, version ="3.6-0")) :
Number of studies (k=7) toosmall to test for small study effects (k.min=10). Change parameter 'k.min' ifappropriate.
这里的small study effects是更广泛意义上的发表偏倚。漏斗图如下：

从该漏斗图中大致可以得出，该meta分析存在发表偏倚的可能性。在这种情况下，可以使用trim and filled或者copas模型等进行校正。以trim and filled为例进行校正，代码如下：
tf1 <- trimfill(metaresult, comb.fixed=TRUE)
summary(tf1)
funnel(tf1)
三行代码分别表示利用trimand filled方法进行校正，并将结果赋予tf1列表中，概括tf1结果，做出填补后的漏斗图。结果如下：

结果表明在漏斗图右侧，填补了3个研究，异质性检验还是没有统计学意义，因此还是采用固定效应模型，OR=0.914，95%CI：0.859-0.973, 同时我们可以看到，随机效应模型的可信区间比固定效应模型要宽，且无统计学意义。经过校正后，合并的效应值还是有统计学意义，在做conclusion的时候更能说明干预的效果。
至于敏感性分析，可以分为很多种，比如比较常见的把risk of bias比较大的（也就是常说的研究质量较低的）研究排除掉做一次敏感性分析，或者将每个研究一次排除掉做敏感性分析。后者可以用meta包中的metainf实现，代码如下：
metainf(metaresult, pooled="fixed")
forest(metainf(metaresult), comb.fixed=TRUE)
如果用随机效应模型，只要将pooled="fixed"改成pooled="random"即可，结果如下：

结果表明，如果排除掉ISIS-2这个研究后，余下的研究合并在一起的效应值将没有统计学意义，OR=0.90, 95%CI: 0.80-1.02，主要原因是该研究由于样本量大，所占的权重也是最大的，而且结果是有统计学意义的，如果将这么大权重的一个有统计学意义的研究排除在外，而剩下的都是些小样本的没有统计学意义的研究，很有可能得出上述的结论，须在discussion中对这点进行强调。
一篇普通meta分析的套路就基本这些，还有就是加几个亚组分析结果或者做meta-regression探索一下异质性的来源，由于数据限制，meta回归的例子由于缺少相应的协变量这里无法举例，等下次完善例子后再实现。

接上周末，主要做一下如何用metareg函数来实现meta回归。实际上在meta分析中，很多地方都可以归结到线性模型的范畴。Meta回归是一种探索异质性来源的方法（本例子的异质性是没有统计学意义的，在这里只是为了方便起见），当然，也可以通过其来校正协变量如年龄等对合并效应量的影响。还是用上述Fleiss93的数据集，在这里自己虚构了一个年龄的协变量。
Fleiss93$age<-c(55, 48, 50, 75, 60, 70, 65)
该语句表示在Fleiss93加入变量age，对7各研究分别赋值为55,48, 50, 75, 60, 70, 65岁的平均年龄。
然后，我们就利用metareg这个函数进行meta回归。代码如下：
library("meta")
data(Fleiss93)
Fleiss93$age<-c(55, 48, 50, 75, 60, 70, 65)
metaresult<-metabin(event.e, n.e, event.c,n.c,data=Fleiss93,sm="OR",
studlab=paste(study, year),comb.random=FALSE)
metareg(metaresult,~age)
基本代码的步骤和上述的发表偏倚等差不多，~age表示协变量，如果有多个协变量的话，可依次加在后面，比如~age+quality+contury等等。需要注意的是，这部分涉及到回归的时候，都需要在取对数以后进行，其实在metaresult列表中，各研究的效应值都已经在log的刻度上了，可以用metaresult$TE查看。
回归的结果如下：

第一部分主要的结果是经过meta回归校正以后的异质性情况，可以看到I^2已经减小到了0.00%，检验也是没有统计学意义的（p=0.5043）。第二部分是关于协变量和校正以后的效应值，该部分结果中所有的效应值都是在log刻度上表示的，如果要转化为OR值，需要使用exp函数，直接在R中用exp（***）表示即可。从上图中可以看出，回归中的年龄是有统计学意义的（p=0.0177）,。如果需要得到校正以后各年龄的合并效应预测值可以使用metafor包中的predict函数，在此不再赘述。
可以做一个bubble图来显示meta回归的结果。代码如下：
reg_age<-metareg(metaresult,~age)
bubble(reg_age)
bubble图如下所示