稀有事件logistic回歸在醫(yī)學(xué)研究中的應(yīng)用*

趙晉芳 羅天娥 范月玲 曾 平 仇麗霞 劉桂芬△

如果疾病的發(fā)生水平很低，極為不常見，病例在人群中所占比重就非常小，那么稱這個(gè)醫(yī)學(xué)事件為稀有事件。如果我們采用常見的現(xiàn)況流行病調(diào)查方法或隊(duì)列研究研究這種疾病，就會(huì)導(dǎo)致收集的數(shù)據(jù)中病例數(shù)與非病例數(shù)很不均衡。比如要探索研究該疾病的影響因素，通常的做法是對(duì)病例和非病例的兩類人群建立logistic回歸模型，然而由于資料中的病例所占的比例遠(yuǎn)遠(yuǎn)低于非病例的比重，這就給稀有事件的統(tǒng)計(jì)分析帶來一系列問題，在這種情況下仍采用常規(guī)的logistic回歸方法就不適合了。本文將主要探討一種適用于解決醫(yī)學(xué)現(xiàn)象研究中稀有事件的logistic回歸模型，它校正了普通logistic回歸在參數(shù)估計(jì)、統(tǒng)計(jì)推斷和概率預(yù)測方面都有可能存在的缺陷。

原理和方法

1．稀有事件logistic回歸

醫(yī)學(xué)研究中，當(dāng)隨機(jī)反應(yīng)變量Y的結(jié)果表現(xiàn)為二分類變量時(shí)，如發(fā)病(Y=1)和不發(fā)病(Y=0)，感染(Y=1)和不感染(Y=0)，若感染(Y=1)的概率P受到因素影響時(shí)，可表示為

x'為暴露因素，α，β分別為截距項(xiàng)和回歸參數(shù)向量。logistic回歸系數(shù)的極大似然估計(jì)值^β具有一致性、漸近有效性和漸近正態(tài)性的性質(zhì)，并且在結(jié)果變量Y兩類取值頻率相等時(shí)的檢驗(yàn)效率是最高的〔1－3〕。

但在稀有事件分析中，由于存在大量Y=0的記錄，而Y=1的例數(shù)卻很少，這就會(huì)導(dǎo)致一般的logistic回歸在參數(shù)估計(jì)、統(tǒng)計(jì)推斷和概率預(yù)測方面都可能存在一定的缺陷。下面介紹一種適合醫(yī)學(xué)中稀有事件的logistic回歸(rare event logistic，re-logistic)，其基本思想是在普通logistic回歸結(jié)果基礎(chǔ)上給予適當(dāng)?shù)男Ｕ?/p>

(1)先驗(yàn)校正

先驗(yàn)校正(prior correction)是在普通logistic回歸最大似然估計(jì)值的基礎(chǔ)上，結(jié)合總體中Y=1的概率τ，以及樣本中Y=1的比例(或叫抽樣概率)ˉy對(duì)回歸系數(shù)的最大似然估計(jì)值進(jìn)行校正〔4〕。~

α為經(jīng)過先驗(yàn)校正的截距項(xiàng)。先驗(yàn)校正的思想最初源于 Prentice和 Pyke(1979)，Manski和 Lerman(1977)，以及Daniel McFadden尚沒有公開發(fā)表的一篇文獻(xiàn)〔5－7〕。先驗(yàn)校正需要已知總體率τ，關(guān)于總體中Y=1的概率τ的先驗(yàn)信息可以從普查、大樣本的隨機(jī)抽樣研究或病例－隊(duì)列研究中得到。

(2)加權(quán)校正

研究中可能存在由于樣本選擇的原因而導(dǎo)致總體概率τ和樣本概率ˉy之間有差異，而加權(quán)校正(weight correction)正是要對(duì)樣本觀察單位給予合適的權(quán)重來補(bǔ)償因選擇偏倚造成的影響。對(duì)樣本中Y=1的觀察單位給予權(quán)重w1=τ/ˉy，Y=0的觀察單位給予權(quán)重w0=(1－τ)/(1－ˉy)。則logistic回歸有以下的加權(quán)對(duì)數(shù)似然函數(shù):

最大化(3)式即可得到參數(shù)的最大似然估計(jì)值。研究表明，加權(quán)校正在大樣本和模型指定有誤時(shí)要優(yōu)于先驗(yàn)校正〔8〕，而在小樣本時(shí)，先驗(yàn)校正要優(yōu)于加權(quán)校正，但這種差別不是很大〔9，10〕。

(3)稀有事件回歸系數(shù)的MCN校正

α和β的最大似然估計(jì)值在小樣本時(shí)是有偏的，而且稀有事件會(huì)進(jìn)一步放大這種偏倚。在小樣本稀有事件中，先驗(yàn)校正和加權(quán)校正仍存在一定的偏差，尚需要進(jìn)行進(jìn)一步的校正。小樣本的稀有事件回歸系數(shù)的偏倚量(bias)〔4，11，12〕:)式中ξi=0.5Qii［(1+w1)^πi－w1］，Qii為矩陣Q=X(X'WX)－1X'的對(duì)角元素，W=diag{^πi(1－^π)wi}。從式(4)可見，實(shí)際上偏倚量bias(^β)就是以X為自變量，ξ為應(yīng)變量，W為權(quán)重的回歸方程的系數(shù)的加權(quán)最小二乘估計(jì)值。校正的參數(shù)估計(jì)值為:

校正的參數(shù)方差矩陣為:樣本的稀有事件回歸系數(shù)的校正不但得到了無偏的參數(shù)估計(jì)量，而且還降低了方差，其統(tǒng)計(jì)性質(zhì)優(yōu)于前者。這種校正方法又被Gary King和Lang che Zeng簡稱為MCN校正(McCullagh Nelder Correction)。

(4)稀有事件概率估計(jì)

稀有事件回歸系數(shù)的最大似然估計(jì)值^β本身是有偏估計(jì)值，因此個(gè)體Y=1的概率估計(jì)也是有偏的;即使^β是無偏估計(jì)值，也并不能保證概率估計(jì)值就是最優(yōu)的?？梢赃x擇下面的公式估計(jì)稀有事件中Y=1的概率:β*為結(jié)合啞變量(integration dummy)。式(7)可以看做 ~β 抽樣分布下 ~P(Y=1|~β)的期望值，而 ~P(Y=1|~β)是Pi=P(Y=1|β)的點(diǎn)估計(jì)值。式(7)亦可以近似表示為:P(Yi=1)≈~Pi+Ci或 P(Yi=1)≈~Pi－Ci。其中，Ci稱為校正因子(correction factor)，計(jì)算公式為:

在滿足一定條件下，~Pi－Ci是近似無偏的，但是模擬研究顯示~Pi+Ci有更小的均方誤。~Pi－Ci稱為Pi=P(Y=1|β)的近似無偏估計(jì)值(approximate unbiased estimator)，~Pi+Ci稱為近似 Bayes估計(jì)值(approximate Bayesian estimator)〔4〕。有研究顯示，除了某些特殊情況，如多個(gè)小樣本的meta分析中，近似無偏估計(jì)值要好于Bayes估計(jì)值外，多數(shù)應(yīng)用中，Bayes估計(jì)值要優(yōu)于近似無偏估計(jì)值。

2．非嵌套模型Vuong檢驗(yàn)

采用Vuong(1989)提出的非嵌套模型檢驗(yàn)(nonnested models test)〔13－15〕來檢驗(yàn) logistic 回歸和稀有事件logistic回歸的非嵌套關(guān)系。

式中，^Prelogit(yi|xi，wi)和 ^Plogit(yi|xi)分別為稀有事件logistic回歸和普通 logistic回歸預(yù)測概率。根據(jù)Vuong，模型1相對(duì)于模型2的非嵌套模型檢驗(yàn)的統(tǒng)計(jì)量為

實(shí)例分析

通過與山西省疾病預(yù)防控制中心聯(lián)合，對(duì)山西省運(yùn)城市五個(gè)項(xiàng)目防治縣的316例HIV/AIDS患者進(jìn)行結(jié)核病的篩查。欲對(duì)HIV/AIDS患者是否患結(jié)核病進(jìn)行分析，結(jié)果如下。Min Max

表1 HIV/AIDS患者資料簡單描述ˉ

n x±s

表2 HIV/AIDS患者資料變量編碼及構(gòu)成

HIV/AIDS患者中大多是初中文化程度，占總患者的66.14%(209/316);患者的平均年齡為41.6歲，以壯年為主;HIV/AIDS患者 CD4計(jì)數(shù)的均值為317.85(個(gè)/μl)，低于正常人水平，其最大值為1125(個(gè)/μl)，最小值為 1(個(gè)/μl)，標(biāo)準(zhǔn)差為 183.80(個(gè)/μl)，變異較大，因此對(duì)CD4作自然對(duì)數(shù)轉(zhuǎn)換，并在以后的分析中代替CD4作為自變量，且仍用CD4作為其變量名。

調(diào)查的316例HIV/AIDS患者中僅有11人是結(jié)核感染者，感染率大約為3.48%。因此我們認(rèn)為分析樣本中HIV/AIDS患者感染結(jié)核是稀有事件。

表3 普通logistic回歸參數(shù)估計(jì)

普通logistic回歸參數(shù)估計(jì)顯示HIV/AIDS患者CD4計(jì)數(shù)有統(tǒng)計(jì)學(xué)意義，與是否感染結(jié)核病有關(guān)系，CD4計(jì)數(shù)對(duì)數(shù)值每增加一個(gè)單位，HIV/AIDS患者感染結(jié)核的危險(xiǎn)性降低71.1%［1－exp(－1.240412)］，即CD4計(jì)數(shù)水平越高HIV/AIDS患者患結(jié)核病的可能性越小。

表4 logistic回歸先驗(yàn)校正

logistic回歸先驗(yàn)校正是在普通logistic回歸參數(shù)估計(jì)的基礎(chǔ)上對(duì)截距項(xiàng)做了校正，其他回歸系數(shù)估計(jì)值和標(biāo)準(zhǔn)誤均未發(fā)生改變。以往報(bào)道表明，感染了結(jié)核菌的HIV/AIDS患者每年發(fā)展為結(jié)核病的機(jī)會(huì)為7%，據(jù)山西省權(quán)威機(jī)構(gòu)2008年底提供的數(shù)字顯示，HIV/AIDS患者的結(jié)核發(fā)病率估計(jì)在5%左右，故本次研究τ取值為0.05。校正后的截距為:

表5 logistic回歸MCN先驗(yàn)校正

logistic回歸MCN先驗(yàn)校正參數(shù)估計(jì)CD4計(jì)數(shù)對(duì)數(shù)值每增加一個(gè)單位，HIV/AIDS患者感染結(jié)核的危險(xiǎn)性降低68.0%。

表6 logistic回歸加權(quán)校正

在普通logistic回歸的基礎(chǔ)上對(duì)樣本中的每個(gè)觀察單位進(jìn)行加權(quán)校正，其中:w1=τ/ˉy=0.05/0.03481013，w0=(1－τ)/(1－ˉy)=(1－0.05)/(1－0.03481013)。

加權(quán)校正logistic回歸表明CD4計(jì)數(shù)對(duì)數(shù)值每增加一個(gè)單位，HIV/AIDS患者感染結(jié)核的危險(xiǎn)性降低70.9%。

加權(quán)l(xiāng)ogistic回歸MCN校正結(jié)果表明CD4計(jì)數(shù)對(duì)數(shù)值每增加一個(gè)單位，HIV/AIDS患者感染結(jié)核的危險(xiǎn)性降低67.8%。

表7 logistic回歸MCN加權(quán)校正

表8 HIV/AIDS患者結(jié)核感染概率估計(jì)

我們采用以上的幾種模型估計(jì)本次調(diào)查樣本的HIV/AIDS患者結(jié)核感染率，不同方法的HIV/AIDS患者結(jié)核感染概率估計(jì)顯示，最大似然估計(jì)的感染概率最小，即普通logistic回歸低估了感染概率，其他三種估計(jì)方法的感染概率估計(jì)值有所提高，彌補(bǔ)了稀有事件中傳統(tǒng)的估計(jì)方法可能會(huì)低估事件Y=1的預(yù)測概率的缺憾，即校正后感染概率估計(jì)偏倚減少。加權(quán)估計(jì)和近似Bayes估計(jì)的感染概率估計(jì)值接近5.00%。

表9 不同校正方法的logistic回歸和普通logistic回歸的Vuong檢驗(yàn)

PPLUS表示近似Bayes估計(jì)概率，PJIAN表示近似無偏估計(jì)概率，PROB表示加權(quán)估計(jì)概率，P0為普通logistic回歸預(yù)測概率。

Vuong檢驗(yàn)顯示加權(quán)校正和近似Bayes估計(jì)都要優(yōu)于普通logistic回歸、近似無偏估計(jì);近似Bayes估計(jì)優(yōu)于加權(quán)校正;而近似無偏估計(jì)不如普通logistic回歸。

綜上，Vuong檢驗(yàn)和概率預(yù)測結(jié)果顯示近似Bayes估計(jì)得到的結(jié)果最優(yōu)。

討 論

醫(yī)學(xué)研究中經(jīng)常遇到二分類反應(yīng)變量資料采用logistic回歸分析，若應(yīng)變量兩類取值頻率相差特別懸殊時(shí)，普通logistic回歸不僅參數(shù)估計(jì)有偏，并且可低估稀有事件的發(fā)生概率。通過稀有事件logistic回歸校正參數(shù)和概率估計(jì)值來解決這個(gè)問題，效果較好。實(shí)例分析結(jié)果表明，在稀有事件的分析中，不管是在模型的整體表現(xiàn)或者是模型的預(yù)測預(yù)報(bào)方面，稀有事件的logistic回歸確實(shí)要更優(yōu)于普通的logistic回歸。因此對(duì)于醫(yī)學(xué)中很多不常見疾病的研究，稀有事件的logistic回歸是一種值得推廣應(yīng)用的統(tǒng)計(jì)模型。當(dāng)然，對(duì)于某一個(gè)醫(yī)學(xué)事件要根據(jù)具體的情況從專業(yè)的角度判斷其是否是稀有或罕見事件，從實(shí)際應(yīng)用看確定這一點(diǎn)并不難。

稀有事件logistic的校正方法可以在Gary King和Langche Zeng 2001年推出的 STATA程序——relogit中實(shí)現(xiàn)，relogit是一種非官方的程序需要下載安裝后才可以使用。不同模型間的Vuong檢驗(yàn)?zāi)壳吧袩o專門的程序，本文通過SAS9.1軟件編程實(shí)現(xiàn)非嵌套模型的比較。

1．Manski CF，McFadden D．Structural analysis of discrete data with econometric applications．MA:Cambridge:MIT Press，1981:51-111．

2．Imbens GW．An efficient method of moments estimator for discrete choice models with choice-based sampling．Econometrica，1992，60(5):1187-1214．

3．陳峰．醫(yī)用多元統(tǒng)計(jì)分析方法．北京:中國統(tǒng)計(jì)出版社，2001．

4．King G，Zeng LC．Logistic regression in rare events data．Political analysis，2001，9(2):137-163．

5．Prentice RL，Pyke R．Logistic disease incidence models and case-control studies．Biometrika，1979，66(3):403-411．

6．Manski CF，Lerman SR．The estimation of choice probabilities from choice based samples．Econometrica，1977，45(8):1977-1988．

7．Manski CF，McFadden D．Structural analysis of discrete data with econometric applications．MA:Cambridge:MIT Press，1981:2-50．

8．Xie Y，Manski CF．The logit model and response-based samples．Sociological Methods Research，1989，17(3):283-302．

9．Amemiya T，Vuong QH．A comparison of two consistent estimators in the choice-based sampling qualitative response model．Econometrica，1987，55(3):699-702．

10．Scott AJ，Wild CJ．Fitting logistic models under tuberculosis-control or choice based sampling．J R Statist Soc B，1986，48(2):170-182．

11．McCullagh P，Nelder JA．Generalized Linear Models．New York:Chapman ＆ Hall，CRC，1999．

12．King G，Zeng LC．Explaining rare events in international relations．International Organization，2001，55(3):693-715．

13．Strazzera E，Genius M．Evaluation of likelihood based tests for non-nested dichotomous choice contingent valuation models．Working Paper CRENoS，2000:2-28．

14．Greene WH．Accounting for excess zeros and sample selection in Poisson and negative binomial regression models．Working Paper No．EC-94-10，Department of Economics，Stern School of Business，New York University，1994．

15．Vuong QH．Likelihood ratio tests for model selection and non-nested hypotheses．Econometrica，1989，57(2):307-333．