單組率研究稀疏數(shù)據(jù)的Meta分析方法

劉 曼 陳文松 劉玉秀,3 陸夢潔 劉雅琦 袁陽丹

近年來，醫(yī)學(xué)研究中對單組率研究的Meta分析呈增長趨勢[1-3]，然而當(dāng)遇到稀疏數(shù)據(jù)時，不顧方法應(yīng)用條件盲目選擇的情況時有發(fā)生，導(dǎo)致合并結(jié)果的偏倚，甚至產(chǎn)生誤導(dǎo)的結(jié)論[4]。有學(xué)者指出，單組率研究的Meta分析若出現(xiàn)稀疏數(shù)據(jù)，尤其是零事件的研究較多時，常用的Meta分析方法統(tǒng)計性能不佳，應(yīng)謹(jǐn)慎選擇[5]。本文針對單組率研究遇到稀疏數(shù)據(jù)時進(jìn)行Meta分析的方法，擇優(yōu)選出目前認(rèn)為具有較好性能的基于Freeman-Tukey轉(zhuǎn)換、反正弦轉(zhuǎn)換的倒方差法和廣義線性混合模型(GLMM)方法，通過Monte-Carlo模擬比較其統(tǒng)計性能，為Meta分析的方法選擇提供依據(jù)、提出建議。

1 方法

Meta分析有3個重要方面，一是分析的框架，二是模型的選擇，三是異質(zhì)性評價方法的選擇。

1.1 框架的選擇 單組率研究的常用Meta分析方法原理通常是基于正態(tài)分布的。對于發(fā)生率為p的總體，例數(shù)為n的樣本中某一事件的發(fā)生數(shù)e服從二項(xiàng)分布，若n足夠大且p不接近0或1，可用二項(xiàng)分布近似正態(tài)分布。然而，當(dāng)n較小或p接近0或1時，由于分布的離散性，則不宜用二項(xiàng)分布近似正態(tài)分布。因此，二項(xiàng)數(shù)據(jù)建模最流行的框架是通過變換使之近似于正態(tài)分布，在進(jìn)行Meta分析的倒方差法中，提出了幾種轉(zhuǎn)換方法，包括ln轉(zhuǎn)換(自然對數(shù)轉(zhuǎn)換)、logit轉(zhuǎn)換[6]、反正弦轉(zhuǎn)換[7]和Freeman-Tukey雙重反正弦轉(zhuǎn)換(FT轉(zhuǎn)換)[8-10]。前2種轉(zhuǎn)換在計算過程中，遇到零事件可能會面臨除零問題，盡管可以做連續(xù)性校正，但研究證明這些校正仍會產(chǎn)生偏倚，難以保證結(jié)果的穩(wěn)定性，而后2種轉(zhuǎn)換則無需校正且可以保持較好的統(tǒng)計性能，因此不少文獻(xiàn)推薦選擇后兩者[6, 7,11]。

單組率的Meta分析還有另外1種分析框架即GLMM，是Stijnen等[10]提出的確切研究內(nèi)似然模型(EWLM)，包括二項(xiàng)-正態(tài)模型(BN)和泊松-正態(tài)模型(PN)等[12]，可方便地實(shí)現(xiàn)稀疏二分類數(shù)據(jù)的Meta分析，他提出以各研究真實(shí)的分布代替其正態(tài)分布的假設(shè)，在單組率Meta分析中，最常用的是BN[13]，該模型包括2個部分，第1部分用二項(xiàng)分布對數(shù)據(jù)進(jìn)行建模，假設(shè)事件數(shù)ei服從參數(shù)為pi，樣本量為ni的二項(xiàng)分布，即

第2部分是logit轉(zhuǎn)換后使用正態(tài)分布來模擬研究之間的異質(zhì)性，即用正態(tài)分布對隨機(jī)效應(yīng)進(jìn)行建模

logit(pi)～Normal(μ,τ2)(2)

其中，μ是總體logit率的均值，τ2是logit尺度下的研究間方差(τ為標(biāo)準(zhǔn)差)。τ為0時則簡化為固定效應(yīng)模型。令 logit(pi)=θi，因此事件數(shù)ei的真實(shí)分布為[5]，

ei～Binomial(ni, exp (θi)1+exp (θi))(3)

1.2 模型的選擇[14]在固定效應(yīng)模型下，分析中的所有研究都具有共同的真實(shí)效應(yīng)，合并效應(yīng)是對共同效應(yīng)大小的估計，其原假設(shè)是共同效應(yīng)為0(對于差值)或1(對于比值)。在隨機(jī)效應(yīng)模型下，假設(shè)研究中的真實(shí)效應(yīng)是從真實(shí)效應(yīng)的分布中抽樣的，原假設(shè)是這些效應(yīng)的平均值為0(對于差值)或1(對于比值)。2種模型的選擇文獻(xiàn)已有詳細(xì)說明[15]，在此不贅述。固定效應(yīng)模型下，若出現(xiàn)0事件，經(jīng)模擬證明，F(xiàn)T轉(zhuǎn)換是優(yōu)選方法，而在隨機(jī)效應(yīng)模型下，因研究間變異差異較大，F(xiàn)T轉(zhuǎn)換可能難以保持穩(wěn)定的統(tǒng)計性能，而在實(shí)際應(yīng)用中，逆方差法應(yīng)用較多，忽視了GLMM。

1.3 異質(zhì)性評價方法的選擇[14]對于研究間方差τ2，目前常用的參數(shù)估計方法有3種。1種由Dersimonian和Laird提出的非迭代法[16]，另外2種是極大似然(ML)方法和限制極大似然(REML)方法。在倒方差法中，這些方法都可以應(yīng)用，對于隨機(jī)效應(yīng)模型，由于極大似然估計會導(dǎo)致方差的低估，故首選REML方法。對于二項(xiàng)式模型的GLMM，由于隨機(jī)效應(yīng)的高維積分的密集計算，無法使用REML方法，因此大多數(shù)軟件都使用ML方法估計研究間方差。

2 Monte-Carlo模擬研究

考慮采用3個度量值評估FT轉(zhuǎn)換、反正弦轉(zhuǎn)換和GLMM在Meta分析應(yīng)用中的統(tǒng)計性能[17]，包括點(diǎn)估計的相對偏倚、置信區(qū)間覆蓋率和平均寬度。所有模擬在SAS 9.4系統(tǒng)環(huán)境下實(shí)現(xiàn)。將每次Meta分析模擬的研究個數(shù)(k)設(shè)定為30，取常用的95%CI固定不變。考慮到樣本量和數(shù)據(jù)結(jié)構(gòu)對結(jié)果的影響，設(shè)置了2個不同的總樣本量(6 000和60 000)，各樣本量下又分平衡和不平衡2種結(jié)構(gòu)。總樣本量為6 000時，平衡結(jié)構(gòu)的每個研究的樣本量設(shè)定為200，不平衡結(jié)構(gòu)的最小樣本量為50，其他樣本以10為增量依次遞增至340。對總樣本量為60 000的情形，相應(yīng)設(shè)置擴(kuò)大10倍即可。模擬研究平均總體率(p)從非常罕見的事件(0.001、0.005、0.01)到相對常見的事件(0.05、0.10、0.30、0.50)變化，研究間標(biāo)準(zhǔn)差(τ)設(shè)置為0.005、0.3、1(logit尺度下)。基于以上模擬參數(shù)設(shè)置，共需模擬2×2×3×3×7=252種場景?；贛eta分析隨機(jī)效應(yīng)模型框架，在不同模擬場景下，進(jìn)行10 000次模擬，模擬分為2步，首先從給定的logit率[logit(p)]和研究間方差(τ2)的正態(tài)分布中模擬真實(shí)的logit率(ηi)，再從二項(xiàng)分布中模擬率為pi=eηi/(1+eηi)、樣本量為ni的Meta分析數(shù)據(jù)集。分別采用不同的變化參數(shù)設(shè)定進(jìn)行Meta分析，獲得相應(yīng)的點(diǎn)估計和95%CI估計。該隨機(jī)模擬過程反復(fù)多次進(jìn)行，直至達(dá)到規(guī)定的模擬次數(shù)，再根據(jù)各次的模擬結(jié)果計算各性能評價指標(biāo)。

3 結(jié)果

在不同的模擬場景下，通過Monte-Carlo模擬Meta分析獲得的不同方法的相對偏倚(圖1)、置信區(qū)間覆蓋率(圖2)和置信區(qū)間平均寬度(圖3)結(jié)果?？傮w來看，在總樣本量相同的情況下，平衡和不平衡2種樣本量結(jié)構(gòu)下的各種方法的統(tǒng)計性能行為都非常接近；在總樣本量不同的情況下，不同方法間的行為表現(xiàn)存在明顯差別。

圖1顯示，總樣本量為6 000時，當(dāng)研究間方差較小(τ<1)時，平均總體率>5%的各種方法的點(diǎn)估計均接近無偏；當(dāng)研究間方差大時(τ=1)，F(xiàn)T轉(zhuǎn)換和反正弦轉(zhuǎn)換只有在率為50%時才接近無偏；而總體率<5%時，只有GLMM繼續(xù)保持無偏?？倶颖玖繛?0 000時，當(dāng)研究間方差較小(τ<1)時，平均總體率>0.5%的各種方法的點(diǎn)估計均接近無偏；當(dāng)研究間方差大時(τ=1)，F(xiàn)T轉(zhuǎn)換和反正弦轉(zhuǎn)換依舊表現(xiàn)較差，而GLMM無論在何種情況都接近無偏。

圖2顯示，樣本量為6 000時，當(dāng)研究間方差較小(τ<1)時，平均總體率>5%的各種方法的置信區(qū)間的覆蓋率均接近名義水平；當(dāng)研究間方差大時(τ=1)，F(xiàn)T轉(zhuǎn)換和反正弦轉(zhuǎn)換只有在總體率>30%時才接近名義水平，而GLMM整體上保持在名義水平附近。樣本量為60 000時，當(dāng)研究間方差較小(τ<1)時，平均總體率>0.5%的各種方法的置信區(qū)間均接近名義水平；當(dāng)研究間方差大時(τ=1)，F(xiàn)T轉(zhuǎn)換和反正弦轉(zhuǎn)換的覆蓋率轉(zhuǎn)換均較差，而GLMM的覆蓋率依舊保持在名義水平。圖3顯示，各方法間的置信區(qū)間的寬度無明顯差別，而隨著平均總體率和研究間方差的增大，各方法的區(qū)間寬度均明顯增大。

圖1 單組率研究Meta分析不同方法Monte-Carlo模擬的相對偏倚結(jié)果矩陣圖圖2 單組率研究Meta分析不同方法Monte-Carlo模擬的置信區(qū)間覆蓋率結(jié)果矩陣圖注 每行對應(yīng)不同樣本量和樣本量結(jié)構(gòu);每列對應(yīng)不同的研究間標(biāo)準(zhǔn)差;連點(diǎn)線對應(yīng)合并方法;橫實(shí)線表示0參照線;橫軸表示不同總體率注 每行對應(yīng)不同樣本量和樣本量結(jié)構(gòu);每列對應(yīng)不同的研究間方差;連點(diǎn)線對應(yīng)不同合并方法;橫實(shí)線表示95%CI參照線;橫軸表示不同總體率

圖3 單組率研究Meta分析不同方法Monte-Carlo模擬的置信區(qū)間平均寬度結(jié)果矩陣圖

4 GLMM的軟件實(shí)現(xiàn)

4.1 SAS實(shí)現(xiàn) SAS中的GLMM是通過PROC NLMIXED實(shí)現(xiàn)的[18]，NLMIXED過程適合包括固定和隨機(jī)效應(yīng)模型的非線性混合模型，考慮到隨機(jī)效應(yīng)，可以為數(shù)據(jù)指定條件分布，包括正態(tài)、二項(xiàng)和泊松分布等，通過NLMIXED過程擬合的模型可以看作是通過MIXED過程擬合的隨機(jī)模型的擴(kuò)展，在MIXED過程中可以通過ML和REML估計參數(shù)，而NLMIXED過程僅能實(shí)現(xiàn)ML。另外由于標(biāo)準(zhǔn)Wald方法沒有考慮τ估計的不確定性，NLMIXED過程的統(tǒng)計基礎(chǔ)是基于t分布的，其自由度df=N-l[19]。

用NLMIXED過程做Meta分析可分別估計固定效應(yīng)和隨機(jī)效應(yīng)模型。固定效應(yīng)模型中，需指定總體率的初始值及取值范圍，并指定事件數(shù)的分布類型。在隨機(jī)效應(yīng)模型中，除了總體率的初始值，還需指定研究間方差的初始值，并指定隨機(jī)項(xiàng)。Meta分析的GLMM的NLMIXED詳細(xì)過程見本文附錄。

4.2 R語言實(shí)現(xiàn) R語言的“meta”包中的“metaprop”函數(shù)可用來做單組率的Meta分析，通過指定“metaprop”函數(shù)中的參數(shù)“method=“GLMM”，R語言就會自動調(diào)用“rma.glmm”函數(shù)擬合GLMM。默認(rèn)情況下，模型中各個系數(shù)的檢驗(yàn)統(tǒng)計量以及相應(yīng)的置信區(qū)間是基于標(biāo)準(zhǔn)正態(tài)分布的，作為一種替代方法，可以設(shè)置“rma.glmm”函數(shù)中的參數(shù)“test ='t'”，該方法基于t分布獲得各個系數(shù)和置信區(qū)間。對于異質(zhì)性檢驗(yàn)，該函數(shù)會自動執(zhí)行2種設(shè)定。①Wald型檢驗(yàn)，用于檢驗(yàn)與飽和模型中添加的啞變量相對應(yīng)的系數(shù)的顯著性。②似然比檢驗(yàn)，檢驗(yàn)相同的系數(shù)集，但是通過計算-2倍固定效應(yīng)和飽和模型的對數(shù)似然差來進(jìn)行。對于由rma.glmm函數(shù)擬合的模型類型，這2個檢驗(yàn)并不相同，甚至可能導(dǎo)致相反的結(jié)論。

GLMM模型不需要計算單個研究的觀察結(jié)果，因此，當(dāng)事件數(shù)為零時，不必在事件計數(shù)中添加常數(shù)，直接令p=0不是問題(所有研究均為零事件的極端情況除外)。原則上，單個研究的權(quán)重可以從每個單獨(dú)研究的似然貢獻(xiàn)中得出，然而，該信息目前在R軟件中不可用，所以在R中無法獲得各研究權(quán)重，且只適用于logit轉(zhuǎn)換，對于單組率的其他轉(zhuǎn)換則只能使用倒方差法。

4.3 STATA實(shí)現(xiàn) STATA中的“metan”命令可用來做Meta分析[20]，在“metan”中，使用基于漸近方差的正態(tài)分布來計算置信區(qū)間。對于率，此區(qū)間可能包含不允許的值，尤其是當(dāng)率接近0或1時。當(dāng)率為0或1時，由于估計的標(biāo)準(zhǔn)誤為零，因此無法計算置信區(qū)間，“metan”命令會自動從合并計算中排除率等于0或1的研究。

Nyaga等[21]在此基礎(chǔ)上開發(fā)出了“metaprop”命令，以補(bǔ)充“metan”，允許使用Wilson得分法和Clooper-Pearson確切法計算95%CI并結(jié)合了FT轉(zhuǎn)換，該命令還可以使用二項(xiàng)分布對研究內(nèi)變異進(jìn)行建模。需要STATA 10或更高的版本，通過鍵入“ssc install metaprop”命令安裝，還可以下載“metaprop”的更新“metaprop_one”，該命令可用來擬合GLMM，更新后的命令需要STATA 13版本。具體語句見附錄。

5 案例分析

一項(xiàng)發(fā)熱嬰兒單核細(xì)胞增多性李斯特菌和腸球菌檢出率的Meta分析[4]就遇到了稀疏數(shù)據(jù)的問題。該研究為了評估由單核細(xì)胞增多性李斯特菌和腸球菌引起嚴(yán)重細(xì)菌感染的患病率，納入了15項(xiàng)研究，共計20 703份血樣，分別對李斯特菌和腸球菌引起菌血癥的患病率做Meta分析，兩者的患病率均不高，因此該Meta分析中出現(xiàn)了大量的零事件，尤其是李斯特菌，除了2項(xiàng)研究各出現(xiàn)1例事件外，其余研究的患病率均為0，用R語言中“meta”包中的“metaprop”函數(shù)對李斯特菌菌血癥發(fā)生率做Meta分析的森林圖見圖4，其中方法選擇GLMM，各研究95%CI方法選擇Clopper-Pearson法。用GLMM擬合的結(jié)果為0.01%(95%CI：0～0.04%)，F(xiàn)T轉(zhuǎn)換的結(jié)果為0.02%(95%CI：0～0.07%)，反正弦轉(zhuǎn)換的結(jié)果為0(95%CI：0～0.01%)。

原文中首先對本組資料的零事件進(jìn)行校正，得到的結(jié)果為0.03%(95%CI：0～0.06%)，該結(jié)果與GLMM的結(jié)果有明顯差別。

圖4 單核細(xì)胞增多性李斯特菌引起嚴(yán)重細(xì)菌感染的患病率森林圖

6 討論

本文討論了單組率研究出現(xiàn)稀疏數(shù)據(jù)時常用的Meta分析方法，包括基于FT轉(zhuǎn)換和反正弦轉(zhuǎn)換的倒方差法和GLMM方法，通過Monte-Carlo模擬，設(shè)置多種模擬分析場景，比較了3種方法的統(tǒng)計性能，結(jié)果發(fā)現(xiàn)GLMM方法具有較好的統(tǒng)計性能且對稀疏數(shù)據(jù)具有穩(wěn)健性，應(yīng)作為單組率研究稀疏數(shù)據(jù)Meta分析的首選方法。同時介紹了隨機(jī)效應(yīng)模型下GLMM方法的原理及其在不同軟件中的實(shí)現(xiàn)，為醫(yī)學(xué)研究者提供了實(shí)用參考。

盡管FT轉(zhuǎn)換可以保留所有研究而無需對零事件做校正，且在固定效應(yīng)模型下，大多數(shù)情況下的FT轉(zhuǎn)換都可以保持優(yōu)良的統(tǒng)計性能，這在之前的模擬中也得到了驗(yàn)證[11]，但當(dāng)研究數(shù)增大且研究中零事件增多時，F(xiàn)T轉(zhuǎn)換的統(tǒng)計性能不佳；在隨機(jī)效應(yīng)模型下，隨著研究間方差的增大，F(xiàn)T轉(zhuǎn)換的性能難以滿足統(tǒng)計學(xué)要求，因此應(yīng)謹(jǐn)慎使用。無疑，GLMM不需要對零事件進(jìn)行校正，其優(yōu)異的統(tǒng)計性能成為單組率研究Meta分析隨機(jī)效應(yīng)模型下的首選。但該方法也并非完美無缺，對于較小的樣本量、極低的事件發(fā)生率和較大的研究間方差，各研究會出現(xiàn)較多的零事件，此時GLMM的點(diǎn)估計會出現(xiàn)偏倚，好在這一問題在極低事件率時并不重要。

事實(shí)上，GLMM在很多軟件系統(tǒng)中都有設(shè)置，例如通過SAS中的PROC NLMIXED過程、R語言的“rma.glmm“函數(shù)和STATA的metaprop_one命令等，均不難實(shí)現(xiàn)。因此，當(dāng)面臨單組率研究稀疏數(shù)據(jù)的Meta分析時，不應(yīng)再受限于軟件實(shí)現(xiàn)的困難，應(yīng)首選GLMM方法。

附錄

1.采用SAS實(shí)現(xiàn)固定效應(yīng)模型Meta分析的方法

proc nlmixed data=Meta_data df=1e6;

parms p=0.5;

bounds 0

model e～binomial(n,p);

estimate "logit(p)" log(p/(1-p));

run;

2.采用SAS實(shí)現(xiàn)隨機(jī)效應(yīng)模型Meta分析的方法

proc nlmixed data=Meta_data DF=1E6 qpoints=100;

parms theta=-4.5 to -0.5 by 0.5 sigmasq=0.3;

pi=exp(logodds)/(1+exp(logodds));

model e ～ binomial(n,pi);

random logodds ～ normal(theta,sigmasq) subject=study;

run;

3.采用R軟件實(shí)現(xiàn)Meta分析的方法

library("meta")

4.采用STATA軟件實(shí)現(xiàn)Meta分析的方法

ssc install metaprop_one

metaprop_one e n, random logit