相信大家已經(jīng)對(duì)R軟件和Rstudio有了初步的了解，“我們學(xué)習(xí)的太復(fù)雜了，希望代碼能簡單些，大家能一起飛……”（創(chuàng)始人語錄）。
#第一部分 感謝“加勒比海帶@Shanghai”老師的代碼和說明 #
#1 代碼來源####
#第二部分 meta分析的實(shí)現(xiàn)#
#2.1總得來說，在R中實(shí)現(xiàn)傳統(tǒng)的meta分析的包中，個(gè)人感覺meta包是比較好的一個(gè)包。在這里將數(shù)據(jù)讀入R等步驟略去，最為新手來說，數(shù)據(jù)導(dǎo)入最簡單的方法是利用R+Rstudio的組合，先利用txt或者CSV格式在你的電腦上將數(shù)據(jù)整理好，然后利用Rstudio中的importdataset導(dǎo)入。
在這里主要利用meta包中現(xiàn)成的數(shù)據(jù)進(jìn)行示例。例子是心梗后服用阿司匹林能否降低死亡率，具體數(shù)據(jù)如下表所示。其中event.e表示治療組的死亡人數(shù)，n.e表示治療組的總?cè)藬?shù)，event.c表示對(duì)照組的總?cè)藬?shù)，n.c表示對(duì)照組的總?cè)藬?shù)，非常簡單的一個(gè)例子。
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函數(shù)可以順利地實(shí)現(xiàn)效應(yīng)量的合并，非常簡單，代碼如下：
library('meta')
data(Fleiss93)
metabin(event.e, n.e,event.c,n.c,data=Fleiss93,sm='OR')
跑出來的結(jié)果如下圖所示：
第一部分是各個(gè)研究的OR值及95%可信區(qū)間，以及分別在固定效應(yīng)模型和隨機(jī)效應(yīng)模型下的權(quán)重。第二部分是合并的結(jié)果，分別列出了固定效應(yīng)模型和隨機(jī)效應(yīng)模型合并的OR值及95%可信區(qū)間，以及檢驗(yàn)的結(jié)果（z值和相應(yīng)的p值），第三部分列出了一致性檢驗(yàn)的結(jié)果，最后一部分是本meta分析所用的具體方法，包括合并效應(yīng)值的方法和計(jì)算研究間方差的方法。
對(duì)于該meta分析，首先需要看的是異質(zhì)性檢驗(yàn)的結(jié)果，這里研究間方差tau^2=0.0096，I^2=40%左右，經(jīng)過檢驗(yàn)異質(zhì)性無統(tǒng)計(jì)學(xué)意義（p=0.127>0.10）,因此可選用固定效應(yīng)模型進(jìn)行效應(yīng)量的合并，結(jié)果是OR=0.897, 95%CI: 0.841-0.957,有統(tǒng)計(jì)學(xué)意義。一般來說，我們只需要看OR值和95%可信區(qū)間就可以判斷是否有統(tǒng)計(jì)學(xué)意義，不需要再看相應(yīng)的p值，OR值的可信區(qū)間不包含1等價(jià)于p。比較一下固定效應(yīng)模型和隨機(jī)效應(yīng)模型的結(jié)果，我們看到隨機(jī)效應(yīng)模型的可信區(qū)間更寬，因此結(jié)果相對(duì)于固定效應(yīng)模型來說更保守，主要原因是相對(duì)于固定效應(yīng)模型，隨機(jī)效應(yīng)模型在計(jì)算標(biāo)準(zhǔn)誤的時(shí)候引入了研究間方差。
接下來的問題是怎么將forest plot（森林圖）做出來，meta分析森林圖一般來說不能少。
代碼如下：
metaresult
studlab=paste(study, year),comb.random=FALSE)
forest(metaresult)
在上段代碼中，多了comb.random=FALSE，代表本例不用隨機(jī)效應(yīng)模型進(jìn)行效應(yīng)量合并，而是采用固定效應(yīng)模型。同理，如果需要用隨機(jī)效應(yīng)模型進(jìn)行合并只要在后面改成comb.fixed=FALSE即可。在這里我們將meta分析的結(jié)果賦給metaresult這個(gè)列表，然后用forest函數(shù)調(diào)用metaresult即可。得到森林圖如下圖所示。
各研究的OR值和95%可信區(qū)間、權(quán)重、效應(yīng)量合并值以及異質(zhì)性檢驗(yàn)結(jié)果等都在該圖中。
這類研究一般是RCT研究，因此選擇RR值可能更好，一般事件發(fā)生率比較低的情況下，RR約等于OR值。如果要使用RR值，只需將sm='OR'改成sm='RR'即可。
metaresult
studlab=paste(study, year),comb.random=FALSE)
到這里基本上將meta分析的主要部分做完了，接下來需要做一個(gè)發(fā)表偏倚檢驗(yàn)和敏感性分析。
發(fā)表偏倚的檢驗(yàn)由于納入的研究個(gè)數(shù)小于10個(gè)，在cochrane手冊中不建議做test，只要求做一個(gè)漏斗圖。當(dāng)研究個(gè)數(shù)大于等于10各的時(shí)候，對(duì)于這種四格表資料不建議使用egger檢驗(yàn)或begg檢驗(yàn)（具體可見我的另外一個(gè)陳年老貼），可以使用Peters檢驗(yàn)，power稍高（所有的發(fā)表偏倚檢驗(yàn)方法在研究個(gè)數(shù)不多的條件下，power都不高），且不需要AS類方法那樣做反正弦變換。
發(fā)表偏倚代碼如下：
funnel(metaresult)
metabias(metaresult,method.bias='peters')
如果納入研究個(gè)數(shù)小于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是更廣泛意義上的發(fā)表偏倚。漏斗圖如下：
從該漏斗圖中大致可以得出，該meta分析存在發(fā)表偏倚的可能性。在這種情況下，可以使用trim and filled或者copas模型等進(jìn)行校正。以trim andfilled為例進(jìn)行校正，代碼如下：
tf1
summary(tf1)
funnel(tf1)
三行代碼分別表示利用trimand filled方法進(jìn)行校正，并將結(jié)果賦予tf1列表中，概括tf1結(jié)果，做出填補(bǔ)后的漏斗圖。結(jié)果如下：
結(jié)果表明在漏斗圖右側(cè)，填補(bǔ)了3個(gè)研究，異質(zhì)性檢驗(yàn)還是沒有統(tǒng)計(jì)學(xué)意義，因此還是采用固定效應(yīng)模型，OR=0.914，95%CI：0.859-0.973, 同時(shí)我們可以看到，隨機(jī)效應(yīng)模型的可信區(qū)間比固定效應(yīng)模型要寬，且無統(tǒng)計(jì)學(xué)意義。經(jīng)過校正后，合并的效應(yīng)值還是有統(tǒng)計(jì)學(xué)意義，在做conclusion的時(shí)候更能說明干預(yù)的效果。
至于敏感性分析，可以分為很多種，比如比較常見的把risk of bias比較大的（也就是常說的研究質(zhì)量較低的）研究排除掉做一次敏感性分析，或者將每個(gè)研究一次排除掉做敏感性分析。后者可以用meta包中的metainf實(shí)現(xiàn)，代碼如下：
metainf(metaresult, pooled='fixed')
forest(metainf(metaresult), comb.fixed=TRUE)
如果用隨機(jī)效應(yīng)模型，只要將pooled='fixed'改成pooled='random'即可，結(jié)果如下：
結(jié)果表明，如果排除掉ISIS-2這個(gè)研究后，余下的研究合并在一起的效應(yīng)值將沒有統(tǒng)計(jì)學(xué)意義，OR=0.90, 95%CI: 0.80-1.02，主要原因是該研究由于樣本量大，所占的權(quán)重也是最大的，而且結(jié)果是有統(tǒng)計(jì)學(xué)意義的，如果將這么大權(quán)重的一個(gè)有統(tǒng)計(jì)學(xué)意義的研究排除在外，而剩下的都是些小樣本的沒有統(tǒng)計(jì)學(xué)意義的研究，很有可能得出上述的結(jié)論，須在discussion中對(duì)這點(diǎn)進(jìn)行強(qiáng)調(diào)。
一篇普通meta分析的套路就基本這些，還有就是加幾個(gè)亞組分析結(jié)果或者做meta-regression探索一下異質(zhì)性的來源，由于數(shù)據(jù)限制，meta回歸的例子由于缺少相應(yīng)的協(xié)變量這里無法舉例，等下次完善例子后再實(shí)現(xiàn)。
接上周末，主要做一下如何用metareg函數(shù)來實(shí)現(xiàn)meta回歸。實(shí)際上在meta分析中，很多地方都可以歸結(jié)到線性模型的范疇。Meta回歸是一種探索異質(zhì)性來源的方法（本例子的異質(zhì)性是沒有統(tǒng)計(jì)學(xué)意義的，在這里只是為了方便起見），當(dāng)然，也可以通過其來校正協(xié)變量如年齡等對(duì)合并效應(yīng)量的影響。還是用上述Fleiss93的數(shù)據(jù)集，在這里自己虛構(gòu)了一個(gè)年齡的協(xié)變量。
Fleiss93$age
該語句表示在Fleiss93加入變量age，對(duì)7各研究分別賦值為55,48, 50, 75, 60, 70, 65歲的平均年齡。
然后，我們就利用metareg這個(gè)函數(shù)進(jìn)行meta回歸。代碼如下：
library('meta')
data(Fleiss93)
Fleiss93$age
metaresult
studlab=paste(study, year),comb.random=FALSE)
metareg(metaresult,~age)
基本代碼的步驟和上述的發(fā)表偏倚等差不多，~age表示協(xié)變量，如果有多個(gè)協(xié)變量的話，可依次加在后面，比如~age+quality+contury等等。需要注意的是，這部分涉及到回歸的時(shí)候，都需要在取對(duì)數(shù)以后進(jìn)行，其實(shí)在metaresult列表中，各研究的效應(yīng)值都已經(jīng)在log的刻度上了，可以用metaresult$TE查看。
回歸的結(jié)果如下：
第一部分主要的結(jié)果是經(jīng)過meta回歸校正以后的異質(zhì)性情況，可以看到I^2已經(jīng)減小到了0.00%，檢驗(yàn)也是沒有統(tǒng)計(jì)學(xué)意義的（p=0.5043）。第二部分是關(guān)于協(xié)變量和校正以后的效應(yīng)值，該部分結(jié)果中所有的效應(yīng)值都是在log刻度上表示的，如果要轉(zhuǎn)化為OR值，需要使用exp函數(shù)，直接在R中用exp（***）表示即可。從上圖中可以看出，回歸中的年齡是有統(tǒng)計(jì)學(xué)意義的（p=0.0177）,。如果需要得到校正以后各年齡的合并效應(yīng)預(yù)測值可以使用metafor包中的predict函數(shù)，在此不再贅述。
可以做一個(gè)bubble圖來顯示meta回歸的結(jié)果。代碼如下：
reg_age
bubble(reg_age)
bubble圖如下所示
對(duì)于結(jié)果的解釋如果需要再深入了解，可以去看一下metafor這個(gè)包中關(guān)于mixedeffect model部分的內(nèi)容（網(wǎng)址：http://www.jstatsoft.org/v36/i03/），寫得很通俗易懂。
###再次感謝“加勒比海帶@Shanghai”老師，謝謝大家###
META專欄主編
猴哥：Freescience公眾號(hào)meta分析欄目現(xiàn)任主編。副主任醫(yī)師，武漢大學(xué)腫瘤學(xué)博士，專注于胃腸道腫瘤分子生物學(xué)機(jī)制、系統(tǒng)評(píng)價(jià)/Meta分析、數(shù)據(jù)挖掘、臨床統(tǒng)計(jì)研究。
科研路，不孤單！
Freescience醫(yī)學(xué)科研聯(lián)盟全國火熱招募ing
50家高校及醫(yī)院的小伙伴已經(jīng)加入啦
具體點(diǎn)這里
想圍觀一篇SCI論文是怎樣寫成的嗎？
Freescience線上沙龍之SCI論文合作呼喚你
具體
點(diǎn)這里FS科研軟件庫，集合50+醫(yī)學(xué)科研必備神器，現(xiàn)在統(tǒng)統(tǒng)打包分享，點(diǎn)這里
還有Freescience科研交流群
