分享

第6章:利用R中的meta包实现meta分析

 chengxiaoli115 2017-07-28

 

 

相信大家已经对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中的importdataset导入。


在这里主要利用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.0096I^2=40%左右,经过检验异质性无统计学意义(p=0.127>0.10,因此可选用固定效应模型进行效应量的合并,结果是OR=0.897, 95%CI: 0.841-0.957,有统计学意义。一般来说,我们只需要看OR值和95%可信区间就可以判断是否有统计学意义,不需要再看相应的p值,OR值的可信区间不包含1等价于p<>。比较一下固定效应模型和随机效应模型的结果,我们看到随机效应模型的可信区间更宽,因此结果相对于固定效应模型来说更保守,主要原因是相对于固定效应模型,随机效应模型在计算标准误的时候引入了研究间方差。
接下来的问题是怎么将forest plot(森林图)做出来,meta分析森林图一般来说不能少。


代码如下:
metaresult<-metabin(event.e, n.e,event.c,n.c,data="">
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="">
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 forsmall study effects (k.min=10). Change parameter 'k.min' ifappropriate.


这里的small study effects是更广泛意义上的发表偏倚。漏斗图如下:


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



结果表明在漏斗图右侧,填补了3个研究,异质性检验还是没有统计学意义,因此还是采用固定效应模型,OR=0.91495%CI0.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,="">
该语句表示在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,="">
metaresult<-metabin(event.e, n.e,event.c,n.c,data="">
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<>
bubble(reg_age)
bubble
图如下所示


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

 

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



META专栏主编


猴哥:Freescience公众号meta分析栏目现任主编。副主任医师,武汉大学肿瘤学博士,专注于胃肠道肿瘤分子生物学机制、系统评价/Meta分析、数据挖掘、临床统计研究。



科研路,不孤单!

Freescience医学科研联盟全国火热招募ing

50家高校及医院的小伙伴已经加入啦

具体点这里

想围观一篇SCI论文是怎样写成的吗?

Freescience线上沙龙之SCI论文合作呼唤你

具体点这里

FS科研软件库,集合50+医学科研必备神器,现在统统打包分享点这里



还有Freescience科研交流群



    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息,谨防诈骗。如发现有害或侵权内容,请点击一键举报。
    转藏 全屏 打印 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多