稀疏二分類數(shù)據(jù)貝葉斯Meta分析方法及實現(xiàn)

張?zhí)灬?/p>

Meta分析可以定量、科學(xué)地合成研究結(jié)果，已在許多科學(xué)領(lǐng)域取得革命性的成果，有助于建立循證實踐、解決相互矛盾的研究結(jié)果問題[1]。眾多醫(yī)學(xué)Meta分析中涉及稀疏二分類數(shù)據(jù)的現(xiàn)象十分常見，即感興趣的測量結(jié)局(如某種干預(yù)措施的不良事件)為二分類數(shù)據(jù)且十分稀疏，特別是在納入Meta分析的單個研究中有1個或2個臂的事件發(fā)生數(shù)為0(分別稱為“單零研究”或 “雙零研究”)，一項對500項Cochrane系統(tǒng)評價隨機抽樣調(diào)查的研究表明，有30%的系統(tǒng)評價中至少有1個研究1個臂的事件發(fā)生數(shù)為0[2]。經(jīng)典的Meta分析方法不適用于稀疏數(shù)據(jù)(sparse data)，特別是納入Meta分析中含有多個單零研究或雙零研究時，在數(shù)據(jù)分析方法學(xué)方面面臨著眾多挑戰(zhàn)，已引起國內(nèi)外研究者的關(guān)注，相繼提出多種模型和方法[3-6]。本文介紹二項式-正態(tài)層次模型(BNHM)和貝塔-二項式模型(BBM)，以及基于這兩種模型框架下的貝葉斯Meta分析方法，并通過實例來介紹軟件實現(xiàn)過程。

1 方法

1.1 研究數(shù)據(jù) 數(shù)據(jù)來源于1項Cochrane系統(tǒng)評價[7]，其主要觀察長效β2受體拮抗劑(LABA)合用糖皮質(zhì)激素吸入劑(ICS)對兒童慢性哮喘的干預(yù)作用，對照組為單純應(yīng)用ICS。本文選取測量結(jié)局為嚴(yán)重不良事件且FEV1平均基線值≥ 80% 預(yù)計值亞組的數(shù)據(jù)，共納入13個研究，其中含有3個單零研究、4個雙零研究，表1顯示，數(shù)據(jù)中變量名study、rtrt、ntrt、rctrl、nctrl分別表示研究名稱、治療組事件發(fā)生人數(shù)和總?cè)藬?shù)、對照組事件發(fā)生人數(shù)和總?cè)藬?shù)。

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

1.2 模型與方法

1.2.1 BNHM 針對二分類數(shù)據(jù)，經(jīng)典Meta分析基于異質(zhì)性假設(shè)可分為固定效應(yīng)和隨機效應(yīng)2種模型，隨機效應(yīng)采用正態(tài)-正態(tài)層次模型(NNHM)最為流行[8],但在某些情況下，該模型假設(shè)難以成立，因此，更多的研究者建議采用BNHM。模擬研究表明，在大多數(shù)情況下，BNHM表現(xiàn)優(yōu)于NNHM，可以獲得無偏倚點估計及合適的置信區(qū)間[4，9]。本文采用的模型由Smith等提出[10]：假設(shè)納入Meta分析第i(i=1,2,...,N)個研究第k(k=0,1)個臂的事件發(fā)生人數(shù)和總?cè)藬?shù)分別為rik和nik，每個臂的事件發(fā)生率為pik，則有：

rik～Bin(pik,nik)

logit(pik)=μi+θixik

θi～N(θ,τ2)

模型中，μi表示第i個研究中事件發(fā)生的基線風(fēng)險，是比值(odds)的對數(shù)尺度，可視為經(jīng)典logistic回歸的截距參數(shù)，為固定效應(yīng)；θ表示平均干預(yù)效應(yīng)，是比值比(OR)的對數(shù)尺度；τ2表示研究間干預(yù)效應(yīng)異質(zhì)性；xik為干預(yù)指示器，如果干預(yù)為治療組(k=1)則其為+0.5，如果干預(yù)為對照組(k=0)則其為-0.5?？梢园l(fā)現(xiàn)，實際上BNHM為廣義線性混合效應(yīng)模型(GLMMs)家族中的一員，不依賴于正態(tài)近似，可以對二項式數(shù)直接建模，非常靈活。

1.2.2 BBM 針對二分類數(shù)據(jù)Meta分析的BBM，是真正的隨機效應(yīng)模型(REM)，實質(zhì)是兩步(two-stage)回歸模型[11，12]。

假設(shè)納入Meta分析第i(i=1,2,...,N)個研究第k(k=0,1)個臂的事件發(fā)生人數(shù)和總?cè)藬?shù)分別為rik和nik，每個臂的事件發(fā)生率為pik，則rik服從帶有3參數(shù)(n,α,β)貝塔-二項式分布(第一階段)：rik～Bin(pik,nik)，pik～Beta(αk,βk)，其中貝塔分布的均數(shù)為μk=αk∕(αk+βk)，可知期望E(rik)=nikμk，相應(yīng)方差為Var(rik)=nikμk(1-μk)[nik+αk+βk)∕(αk+βk+1)]，假定所有研究中治療組和對照組相關(guān)系數(shù)均相等，均為ρ=1∕(αk+βk+1)。

第二階段，Meta分析模型為：g(g)=b0+θxk。式中，g(g)為廣義線性模型家族常用的連接函數(shù)，根據(jù)效應(yīng)量不同選擇不同連接函數(shù)。如效應(yīng)量為OR，則采用logit連接函數(shù)；如效應(yīng)量為RR，則采用對數(shù)連接函數(shù)。θ表示平均干預(yù)效應(yīng)；xt為干預(yù)指示器，xt=0表示干預(yù)為對照組，xt=1表示干預(yù)為治療組。

因該模型所具有的閉型對數(shù)似然函數(shù)有利于參數(shù)估計，可以比較輕松地在頻率學(xué)框架或貝葉斯框架的軟件實現(xiàn)，如前者可以采用SAS軟件的NLMIXED程序?qū)崿F(xiàn)[11]，后者可以采用R軟件的MetaStan包實現(xiàn)。

1.3 模型擬合 本文采用完全貝葉斯策略擬合BNHM(Smith模型)[10]，需要對μi、θ、τ的3個參數(shù)指定先驗分布。對于μi指定為無信息先驗分布N(0,102)[13]，對于θ分別指定為無信息先驗分布N(0,1002)[14]及弱信息先驗分布N(0,2.822)[15]，對于τ分別指定為均勻分布uniform(0,5)[14]、半正態(tài)分布HN(0.5)[15,16]及半柯西分布HC(0.5)等，共6種模型組合，具體如表2所示。

表2 不同參數(shù)的先驗分布設(shè)定

BNHM和BBM均通過R軟件(ver3.4.4)的MetaStan擴展包(ver0.1.0)[17]擬合進行貝葉斯Meta分析，迭代50 000次，前20 000次用于退火以消除初始值的影響；采用收縮因子(shrink factor)進行收斂性判斷，如果每個參數(shù)的收縮因子接近于1，則認(rèn)為是馬爾可夫鏈?zhǔn)諗縖18]。具體代碼及簡要解釋見表3和表4。

表3 基于BNHM稀疏二分類數(shù)據(jù)貝葉斯Meta分析方法代碼(R軟件)及簡要解釋

表4 基于BBM稀疏二分類數(shù)據(jù)貝葉斯Meta分析方法代碼(R軟件)及簡要解釋

2 結(jié)果

2.1 BNHM擬合結(jié)果 基于貝葉斯分析框架下，擬合BNHM主要結(jié)果如表5所示，其中3個參數(shù)的收縮因子(結(jié)果中Rhat值顯示)均為1，表明馬爾可夫鏈已收斂；基于BNHM框架，對參數(shù)μi、θ、τ指定相應(yīng)的先驗分布，擬合不同的模型結(jié)果如表2所示?？梢园l(fā)現(xiàn)，不同模型中最主要的2個參數(shù)干預(yù)效應(yīng)參數(shù)θ和異質(zhì)性參數(shù)τ的后驗均數(shù)及95%CI十分相近，但將τ指定為均勻分布時獲得的后驗均數(shù)估計值較其他先驗分布更大且95%CI更寬。

2.2 BBM擬合結(jié)果 基于貝葉斯分析框架下，擬合BBM主要結(jié)果如表6所示，模型中9個參數(shù)的收縮因子(結(jié)果中Rhat值顯示)均為1，表明馬爾可夫鏈已收斂。主要感興趣的參數(shù)是平均干預(yù)效應(yīng)θ(lnOR)，后驗分布均數(shù)點估計及95%CI為0.12(-0.83，1.09)，返回OR及其95%CI為1.13(0.44，2.97)。

表5 基于不同先驗分布稀疏二分類數(shù)據(jù)貝葉斯Meta分析結(jié)果

表6 擬合BBM主要參數(shù)的后驗分布結(jié)果

3 討論

在醫(yī)學(xué)領(lǐng)域內(nèi)，雖然將發(fā)生概率<1/1 000的事件肯定為稀有事件，但更多的是將發(fā)生概率< 1/100定義為“稀有”[5，6]，如觀察某些藥物少見或罕見的不良反應(yīng)，此類數(shù)據(jù)稱為稀疏數(shù)據(jù)或稀有數(shù)據(jù)(rare data)。對于稀有事件，單個研究結(jié)果通常不足以發(fā)現(xiàn)不同干預(yù)措施效果的差異，而 Meta 分析則可以通過匯總來自多個臨床試驗的數(shù)據(jù)，增加樣本量和統(tǒng)計效能，可能是獲得證據(jù)的唯一方法[6]，但經(jīng)典Meta分析方法學(xué)并不一定適用于稀疏數(shù)據(jù)，需要進一步探索和研究相關(guān)的分析方法。

目前常用于稀疏二分類數(shù)據(jù)Meta分析的方法主要有Peto法(Peto's method)、Mantel-Haenszel法、干預(yù)臂連續(xù)校正(treatment-arm continuity correction)法、Logistic 回歸、貝葉斯Meta分析、BBM等，并有文獻就相關(guān)方法進行了優(yōu)劣比較，有興趣的讀者可以研讀相關(guān)文獻[5]。簡而言之，經(jīng)典方法(如Peto法、Mantel-Haenszel法等)和通用Meta分析軟件一般是對單零研究進行連續(xù)性校正(如對四格表數(shù)據(jù)每格的數(shù)據(jù)各加0.5)，而將雙零研究排除在外，不進行合并；Logistic 回歸雖然不需要對單零研究進行連續(xù)性校正，但也會將雙零研究除外，不納入分析。這些處理方法存在的主要問題有[4，5]：①采用“連續(xù)性校正”對零事件進行修正，對稀疏數(shù)據(jù)的 Meta 分析會造成結(jié)果偏倚；②忽略了雙零研究通過樣本量而提供的相關(guān)信息；③雙零研究不納入分析可能不符合“倫理”。

貝葉斯統(tǒng)計學(xué)和經(jīng)典統(tǒng)計學(xué)是當(dāng)前兩大主要統(tǒng)計學(xué)學(xué)派，Meta分析也主要有頻率學(xué)和貝葉斯兩大策略，尤其是貝葉斯策略，隨著現(xiàn)代計算機技術(shù)的出現(xiàn)和發(fā)展，特別是馬爾可夫鏈蒙特卡洛(MCMC)方法的建立，成功解決了高維積分運算等制約貝葉斯統(tǒng)計發(fā)展問題，近年來在Meta分析特別是處理復(fù)雜或特殊數(shù)據(jù)如稀疏數(shù)據(jù)時應(yīng)用越來越廣泛。針對二分類數(shù)據(jù)，貝葉斯策略主要有NNHM和BHNM兩種建?？蚣?，后者更為合適[19]?；贐HNM框架下有多個隨機效應(yīng)模型，模擬研究表明以Smith模型為優(yōu)[10，20]。貝葉斯統(tǒng)計方法綜合未知參數(shù)的總體信息、樣本信息及先驗信息，根據(jù)貝葉斯定理，獲得未知參數(shù)的后驗分布，進而對未知參數(shù)進行統(tǒng)計推斷，在實踐中對先驗信息的利用和指定非常重要。在BHNM框架下，Smith模型中[10]不同參數(shù)的先驗分布可以允許利用先驗信息確定，本研究參考既往文獻，一般將基線風(fēng)險參數(shù)μi指定為無信息先驗正態(tài)分布N(0,102)[13]。對干預(yù)效應(yīng)參數(shù)θ一般可指定為無信息先驗(non-informative prior)，如均數(shù)為0、具有很大方差的正態(tài)分布N(0,1002)[14]；也有研究者基于Cochrane系統(tǒng)評價計算出1個弱信息先驗(WIP)，如均數(shù)為0、方差約為8的正態(tài)分布N(0,2.822)[15]。對于異質(zhì)性參數(shù)τ一般指定為均勻分布uniform(0,5)[14]，但基于此的先驗信息獲得的τ后驗均數(shù)較大、95%CI更寬；鑒于τ值為0.25、 0.5、1和 2分別表示異質(zhì)性中等、較大、大、很大，而模擬結(jié)果顯示半正態(tài)分布HN(0.5)獲得的τ后驗均數(shù)與眾多Meta分析結(jié)果更為常見與接近，所以HN(0.5)是一個非常好的先驗選擇[15,16]，而半柯西分布HC(0.5)則可以作為敏感性分析，這些參數(shù)的先驗設(shè)定對今后稀疏二分類數(shù)據(jù)Meta分析實踐具有一定的指導(dǎo)意義。

BNHM、BBM等幾種新的模型和方法為合并稀疏二分類數(shù)據(jù)提供了更為合理的分析策略。一項模型研究表明[11]，測量結(jié)局為稀疏事件的Meta分析，在干預(yù)效應(yīng)估計偏倚和準(zhǔn)確性方面，BBM較其他模型為優(yōu)。實際上，BNHM和BBM均為兩階段REM，但對分布假設(shè)有異同，對事件兩種模型均假定服從二項式分布，而對事件概率前者認(rèn)為服從正態(tài)分布(經(jīng)數(shù)據(jù)轉(zhuǎn)換)，后者認(rèn)為服從貝塔分布，正因為兩種模型對隨機效應(yīng)的分布假設(shè)不同，所合并的結(jié)果有可能會有差異。在BBM框架下，可以選擇OR、危險比(RR)以及危險差(RD)等作為干預(yù)效應(yīng)量，而Bakbergenuly等[21]建議BBM采用OR為干預(yù)效應(yīng)量。

BNHM和BBM等模型均可以通過WinBUGS、OpenBUGS、JAGS、Stan等貝葉斯軟件包來實現(xiàn)，但需要系統(tǒng)評價員編寫相關(guān)代碼，以及設(shè)定不同馬爾科夫鏈的初始值設(shè)定等，而且應(yīng)用MCMC方法等均需要具備一定的專業(yè)知識(如收斂性診斷)和實踐經(jīng)驗等，這對于非統(tǒng)計學(xué)專業(yè)人員而言有一定的難度。有鑒于此，Günhan 編寫的MetaStan包，可以通過R軟件調(diào)用Stan軟件，從而避免編寫Stan軟件實現(xiàn)模型代碼(較R軟件代碼編寫難度更高)，輕松地擬合BNHM及BBM等模型而獲得相應(yīng)的結(jié)果。此外，收斂性診斷是MCMC算法非常重要的問題[22]，MetaStan包主要采用的是收縮因子法，在實踐中，還可以與shinystan包聯(lián)用，通過繪制和觀察Gibbs動態(tài)抽樣圖、迭代歷史圖、抽樣自相關(guān)圖，以及參數(shù)后驗分布的核密度估計圖等，判斷馬爾可夫鏈?zhǔn)欠袷諗縖10]，鑒于篇幅，本研究不再展示，有興趣的讀者可以自行研究和運用。

基于以下優(yōu)點，針對稀疏二分類數(shù)據(jù)進行Meta分析，推薦使用基于BNHM和BBM等模型框架下的貝葉斯方法：①建模靈活，并在貝葉斯分析軟件中可實現(xiàn)性強，如采用R軟件的MetaStan包可以輕松擬合BNHM和BBM；②充分考慮了模型中參數(shù)如方差等的不確定性，合并效應(yīng)量點估計及95%CI來源于參數(shù)的后驗分布，結(jié)果可靠；③二項式似然可以直接對稀疏數(shù)據(jù)，特別是單零和雙零研究數(shù)據(jù)直接建模，不需要對數(shù)據(jù)進行連續(xù)性校正，降低了經(jīng)典Meta分析因采用連續(xù)性校正造成的偏倚。

值得注意的是，雖然目前有眾多的稀疏二分類數(shù)據(jù)Meta分析模型和方法，但尚無定論明確回答哪種方法是最佳的稀疏數(shù)據(jù)Meta分析策略，不同的方法基于不同的假說，其效度難以評估。因此，系統(tǒng)評價員在制定系統(tǒng)評價/Meta分析計劃時就要考慮到，如果遇到稀疏數(shù)據(jù)時，應(yīng)選擇多種分析方法進行敏感性分析，以確定結(jié)果是否穩(wěn)健。