基于二項(xiàng)式-正態(tài)層次模型框架下比例的貝葉斯Meta分析方法及實(shí)現(xiàn)

張?zhí)灬?/p>

Meta分析可以定量、科學(xué)地合成研究結(jié)果，已在許多科學(xué)領(lǐng)域取得革命性的成果，有助于建立循證實(shí)踐、解決相互矛盾的研究結(jié)果問(wèn)題[1]。經(jīng)典的Meta分析是針對(duì)兩種干預(yù)措施基于頭對(duì)頭(head to head)的直接比較，主要觀察兩種干預(yù)措施的效果差異，相關(guān)方法已較成熟。近年來(lái)，國(guó)內(nèi)外出現(xiàn)大量的系統(tǒng)評(píng)價(jià)采用Meta分析方法合并比例(proportions)或率(rates)[2-6]，但這些文獻(xiàn)均是基于經(jīng)典的正態(tài)-正態(tài)層次模型(normal-normal hierarchical model，NNHM)，獲得的結(jié)果可能會(huì)存在偏倚，實(shí)際情況下，某些比例或率等數(shù)據(jù)(如接近于0或1時(shí))并不適用于該模型。本文介紹一種更適用于比例的Meta分析模型及方法，即二項(xiàng)式-正態(tài)層次模型(binomial-normal hierarchical model，BNHM)和基于該模型框架下貝葉斯Meta分析方法，并通過(guò)實(shí)例來(lái)介紹軟件實(shí)現(xiàn)過(guò)程。

1 方法

1.1 研究數(shù)據(jù) 數(shù)據(jù)來(lái)源于2個(gè)研究。數(shù)據(jù)1：來(lái)自Saber等的Meta分析文獻(xiàn)[4]，共納入29個(gè)研究(所有研究比例均<0.3，5個(gè)研究比例為0)，觀察Pipeline栓塞裝置(Pipeline embolization device，PED)治療顱內(nèi)動(dòng)脈瘤的不良反應(yīng)，本文選取擬分析的數(shù)據(jù)為缺血性并發(fā)癥測(cè)量結(jié)果。數(shù)據(jù)2：來(lái)自Pritz等的Meta分析文獻(xiàn)[7]，共納入14個(gè)研究(10個(gè)研究比例>0.7,2個(gè)研究比例為1)，觀察高動(dòng)力療法(hyperdynamic therapy)治療血管痙攣的療效，以神經(jīng)功能缺失的臨床癥狀改善為治療成功的標(biāo)準(zhǔn)。表1中“n”和“r”分別表示每個(gè)研究總?cè)藬?shù)和事件發(fā)性(缺血并發(fā)癥或治療成功)人數(shù)。

表1 分析數(shù)據(jù)

該模型可以SAS、Stata和R等軟件來(lái)擬合。但在實(shí)踐中會(huì)遇到很多問(wèn)題[9,10]：①如果某些研究的事件發(fā)生比例很小(如接近于0)、或很大(如接近于1)，則抽樣模型正態(tài)分布假設(shè)是否滿足值得懷疑，且因研究估計(jì)值的方差會(huì)趨向于0，則該研究會(huì)得到一個(gè)很大的權(quán)重，計(jì)算所得的研究估計(jì)值95%CI可能會(huì)在[0,1]之外，此時(shí)需要采用logit轉(zhuǎn)換法或雙反正弦法等方法對(duì)數(shù)據(jù)進(jìn)行轉(zhuǎn)換；②當(dāng)事件發(fā)生比例為0或1時(shí)，則方差為0，無(wú)法使用IV法計(jì)算，要對(duì)數(shù)據(jù)進(jìn)行連續(xù)性校正，如對(duì)0格子數(shù)據(jù)加0.5。

在實(shí)踐中，上述數(shù)據(jù)轉(zhuǎn)換或連續(xù)性校正等方法均會(huì)導(dǎo)致結(jié)果偏倚。因此，有學(xué)者提出BNHM用于合并比例數(shù)據(jù)[8,9]，并給出了貝葉斯隨機(jī)效應(yīng)模型[11，12]，具體如下。

ri～Binomial(pi, ni)

1.3 模型擬合 BNHM可由貝葉斯分析專用軟件包來(lái)擬合，前期研究曾經(jīng)給出了傳統(tǒng)的貝葉斯分析軟件WinBUGS實(shí)現(xiàn)該模型的代碼[11，13]。針對(duì)實(shí)例數(shù)據(jù)，本文介紹新近比較流行的貝葉斯分析軟件JAGS(ver4.3.0)聯(lián)合R(ver3.4.4)實(shí)現(xiàn)該模型的擬合[14]，具體代碼及簡(jiǎn)要解釋見(jiàn)表2。通過(guò)R軟件的R2jags擴(kuò)展包來(lái)調(diào)用JAGS軟件進(jìn)行貝葉斯Meta分析，迭代50 000次，前20 000次用于退火以消除初始值的影響；采用收縮因子(shrink factor)進(jìn)行收斂性判斷，如果每個(gè)參數(shù)的收縮因子接近于1則認(rèn)為馬爾可夫鏈?zhǔn)諗縖15]；mcmcplots擴(kuò)展包繪制森林圖(在BUGS語(yǔ)言軟件中稱為caterpillar圖，實(shí)質(zhì)上與Meta分析中森林圖無(wú)異)；為了便于比較，本文同時(shí)采用R軟件的meta擴(kuò)展包來(lái)實(shí)現(xiàn)在頻率學(xué)框架下擬合NNHM和BNHM，具體代碼及簡(jiǎn)要解釋見(jiàn)表3。

表2 比例的貝葉斯Meta分析代碼(R軟件)及簡(jiǎn)要解釋

表3 比例的頻率學(xué)Meta分析代碼(R軟件)及簡(jiǎn)要解釋

2 結(jié)果

2個(gè)研究數(shù)據(jù)貝葉斯分析結(jié)果顯示，每個(gè)參數(shù)的收縮因子(結(jié)果中Rhat值顯示)都接近于1(結(jié)果從略)，表明馬爾可夫鏈已收斂；不同數(shù)據(jù)貝葉斯方法獲得的偏離信息準(zhǔn)則(deviance information criterion,DIC)值分別為146.4和69.5，基于不同模型、采用不同分析框架和方法，隨機(jī)效應(yīng)模型獲得的數(shù)字化結(jié)果見(jiàn)表4，基于BNHM，不論采用貝葉斯方法還是頻率學(xué)方法，總的合并比例結(jié)果點(diǎn)估計(jì)及95%CI相差不大，但貝葉斯方法估計(jì)的研究間方差(τ2值)較大；基于NNHM獲得合并效應(yīng)量(點(diǎn)估計(jì)及95%CI)較BNHM結(jié)果相差比較大，特別是事件發(fā)生率小的數(shù)據(jù)(如數(shù)據(jù)1)；基于NNHM的算法獲得的τ2值均較BNHM低，該模型低估了研究間方差。繪制森林圖(圖1和2)，圖中所示細(xì)線表示各個(gè)研究和總的合并結(jié)果的95%CI，粗線表示后驗(yàn)分布的上、下四分位數(shù)間距離。

表4 基于不同模型、不同分析框架下Meta分析結(jié)果

圖1 數(shù)據(jù)1森林圖(caterpillar圖)

圖2 數(shù)據(jù)2森林圖(caterpillar圖)

3 討論

單臂研究在醫(yī)學(xué)研究中較為常見(jiàn)，所獲得的觀察指標(biāo)按資料性質(zhì)可以分為定量資料和定性資料，其中比例即是常用的定性資料指標(biāo)。所謂比例是指某事物內(nèi)部各組成部分的觀察單位數(shù)與所有組成部分的總觀察單位數(shù)之比，可以表示分布結(jié)構(gòu)的比例，也可以表示現(xiàn)象發(fā)生的頻率，與觀察時(shí)間單位無(wú)關(guān)，醫(yī)學(xué)統(tǒng)計(jì)學(xué)中有很多相對(duì)指標(biāo)被稱為“率”，但實(shí)質(zhì)上是“比例”，如患病率就是符合比例定義，但名稱為率[11]。

目前發(fā)表的眾多關(guān)于比例的系統(tǒng)評(píng)價(jià)文獻(xiàn)多是基于NNHM，在頻率學(xué)框架下采用IV法合并。一般認(rèn)為，如果每個(gè)研究總樣本量ni足夠大(如>100)，nipi與ni(1-pi)不太小(如均>5)，則可以直接采用IV法合并；如果pi>0.7或pi<0.3則需要對(duì)數(shù)據(jù)進(jìn)行轉(zhuǎn)換[11]，常用的方法有Freeman-Tukey轉(zhuǎn)換[16]及l(fā)ogit轉(zhuǎn)換等方法，這些方法均可由Stata軟件的metaprop命令及R軟件的meta包和metafor包等輕松實(shí)現(xiàn)。如果比例為0或1的研究，則方差為0，采用直接合并數(shù)據(jù)或進(jìn)行l(wèi)ogit數(shù)據(jù)轉(zhuǎn)換后再合并，均需要對(duì)數(shù)據(jù)進(jìn)行連續(xù)校正，一般情況下各種軟件默認(rèn)加0.5校正，這樣會(huì)高估或低估合并結(jié)果；而雙反正弦數(shù)據(jù)轉(zhuǎn)換(也稱為Freeman-Tukey轉(zhuǎn)換)不需要校正，但如果納入Meta分析的研究總?cè)藬?shù)<10個(gè)(如數(shù)據(jù)2)，雙反正弦數(shù)據(jù)轉(zhuǎn)換后合并的結(jié)果再返回比例時(shí)，則有可能出現(xiàn)誤導(dǎo)性結(jié)果。本研究發(fā)現(xiàn)，無(wú)論是否進(jìn)行數(shù)據(jù)轉(zhuǎn)換，NNHM會(huì)明顯低估研究間方差；且如果研究比例接近甚至達(dá)到0或1，該模型估算可信區(qū)間會(huì)出現(xiàn)偏倚，這與Hamza等[10]模擬研究結(jié)果完全一致。

模擬研究表明，在大多數(shù)情況下，BNHM表現(xiàn)優(yōu)于NNHM，可以獲得無(wú)偏倚點(diǎn)估計(jì)及合適的置信區(qū)間[10]。該模型實(shí)質(zhì)上是廣義線性模型，可以在頻率學(xué)和貝葉斯框架下實(shí)現(xiàn)，即使比例為0或1，也不需要進(jìn)行連續(xù)校正，前者可以使用R軟件的meta包和metafor包等頻率學(xué)軟件包來(lái)實(shí)現(xiàn)，后者可以通過(guò)WinBUGS、OpenBUGS、JAGS和Stan等貝葉斯軟件包來(lái)實(shí)現(xiàn)。R軟件中有不少擴(kuò)展包，如R2WinBUGS、R2OpenBUGS、R2jags、rstan等，可以分別調(diào)用WinBUGS等貝葉斯軟件包，則操作更靈活、顯示結(jié)果更豐富。貝葉斯統(tǒng)計(jì)方法綜合未知參數(shù)的總體信息、樣本信息及先驗(yàn)信息，根據(jù)貝葉斯定理，獲得未知參數(shù)的后驗(yàn)分布，進(jìn)而對(duì)未知參數(shù)進(jìn)行統(tǒng)計(jì)推斷，目前已廣泛應(yīng)用于復(fù)雜數(shù)據(jù)的Meta分析(如網(wǎng)狀Meta分析等)。本文是在綜合了前人研究工作基礎(chǔ)上，就如何合并比例數(shù)據(jù)，提供了基于BNHM框架下的貝葉斯Meta分析方法，因?yàn)镸eta分析中固定效應(yīng)模型假設(shè)研究之間具有同質(zhì)性，但實(shí)際上一般很難成立[17]，本文使用的研究數(shù)據(jù)也證明了這一點(diǎn)，表4和圖1、2顯示，無(wú)論貝葉斯還是頻率學(xué)框架，各種不同模型的分析結(jié)果均說(shuō)明，納入Meta分析的研究間均存在較大的異質(zhì)性，因此本文是基于更為合適的隨機(jī)效應(yīng)來(lái)建模，并以2個(gè)特殊、但在實(shí)踐中卻是常見(jiàn)的數(shù)據(jù)類型為例，提出R軟件實(shí)現(xiàn)方法。事實(shí)上BNHM模型總優(yōu)于NNHM模型，NNHM在很多情況下不適用[10]，因此，在Meta分析時(shí)應(yīng)優(yōu)先使用BNHM。在實(shí)踐中，貝葉斯方法對(duì)先驗(yàn)信息的利用和指定非常重要，貝葉斯Meta分析的先驗(yàn)分布可以允許利用先驗(yàn)信息確定，但在先驗(yàn)信息或知識(shí)缺乏的情況下，可以按下列方法處理：對(duì)絕對(duì)或相對(duì)效應(yīng)測(cè)量(位置參數(shù))選擇指定無(wú)信息或模糊正態(tài)先驗(yàn)分布[18]；對(duì)研究精度指定為伽馬先驗(yàn)分布，或?qū)ρ芯块g質(zhì)異質(zhì)性指定為均勻先驗(yàn)分布。請(qǐng)注意，任何基于馬爾可夫鏈蒙特卡羅(Markov Chain Monte Carlo，MCMC)的推斷都是在假定馬爾可夫鏈已經(jīng)到達(dá)穩(wěn)定狀態(tài)(收斂)下進(jìn)行的，因此，收斂性診斷是MCMC算法非常重要的問(wèn)題[19]，本文采用的是收縮因子法， 還可以通過(guò)觀察Gibbs動(dòng)態(tài)抽樣圖、迭代歷史圖、抽樣自相關(guān)圖以及參數(shù)后驗(yàn)分布的核密度估計(jì)圖等，判斷馬爾可夫鏈否收斂。

需要指出的是，在頻率學(xué)框架下，比例的置信區(qū)間估算有簡(jiǎn)單近似法、連續(xù)相關(guān)校正簡(jiǎn)單近似法、得分法、連續(xù)校正得分法、精確法等方法[11]，多數(shù)教科書(shū)多推薦精確法，但鑒于精確法比較保守，有研究模擬顯示在某些時(shí)候近似法較精確法為優(yōu)[20]；而貝葉斯框架下，比例的可信區(qū)間則來(lái)源于未知參數(shù)的后驗(yàn)分布。

基于以下優(yōu)點(diǎn)，針對(duì)比例進(jìn)行Meta分析，本文推薦使用基于BNHM框架下的貝葉斯方法，其優(yōu)點(diǎn)包括：①建模靈活，并在貝葉斯分析軟件中可實(shí)現(xiàn)性強(qiáng)；②充分考慮了模型中參數(shù)(如方差等)的不確定性，合并效應(yīng)量點(diǎn)估計(jì)及95%可信區(qū)間來(lái)源于參數(shù)的后驗(yàn)分布，結(jié)果可靠；③二項(xiàng)式似然允許有0格子數(shù)據(jù)，所以即使有研究的比例為0或1，也不需要與經(jīng)典頻率學(xué)方法一樣對(duì)數(shù)據(jù)進(jìn)行連續(xù)性校正。