meta分析

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

我们仅仅是代码的编辑者、整合者、搬运工,仅免费传授方法,下文数据和代码取自于网络和免费软件“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图如下所示

对于结果的解释如果需要再深入了解,可以去看一下metafor这个包中关于mixed effect model部分的内容(网址:http://www.jstatsoft.org/v36/i03/),写得很通俗易懂。

再次感谢“加勒比海带@Shanghai”老师,谢谢大家###

 

拓展知识:

SNPmeta

SNPMETA1

SNPMETA2_检索式

SNPMETA3_索要原始文献+数据提取

 

教你如何做一篇meta分析1_

教你如何做一篇meta分析2_

教你如何做一篇meta分析3_哈温平衡的计算

教你如何做一篇meta分析4_NOS量表的制作

 

网状emta分析必备技能

网状emta分析必备技能1

网状meta分析必备技能2_R软件gemtc包

网状Meta分析必备技能3_数据初始值的设定

网状Meta分析必备技能4_Gemtc实现生存分析

网状meta分析必备技能5_HR风险比

网状meta分析必备技能6_R+Rstudio运用meta包做简单meta分析

网状meta分析必备技能7_使用R、GeMTC和STATA软件实现连续变量的网状Meta分析

 

网状meta分析stata简易教程

NMA(网状meta分析)stata简易教程(1)

NMA(网状meta分析)stata简易教程(2)网状图

NMA(网状meta分析)stata简易教程(3)贡献图

NMA(网状meta分析)stata简易教程(4)漏斗图的制作

NMA(网状meta分析)stata简易教程(5)排序图的制作

 

诊断性meta分析简单实现

诊断性meta分析简单实现(1)之revman

诊断性meta分析简单实现(2)

诊断性meta分析简单实现(3)

诊断性meta分析简单实现(4)_meta-disc

诊断性meta分析简单实现(5)_结果的展示和解读

 

其他

meta分析统计方法-随机对照试验的风险评估rob

生物标记物(基因)联合诊断模型的stata实现ROC(AUC)

生物标记物(基因)联合诊断模型的R实现ROC(AUC)

(1)

关键词:

热评文章

发表评论