王俊蘋(píng) 沈靈智 盧肇駿 寇 碩 鄭衛(wèi)軍

1 浙江中醫(yī)藥大學(xué)公共衛(wèi)生學(xué)院，310000 浙江 杭州；2 浙江省疾病預(yù)防控制中心，310000 浙江 杭州

前瞻性研究以二分類(lèi)變量作為研究結(jié)局是常見(jiàn)的現(xiàn)象，研究者習(xí)慣采用logistic回歸獲得的優(yōu)勢(shì)比(odds ratio，OR)反映暴露對(duì)結(jié)果的關(guān)聯(lián)性；但是在隊(duì)列研究中更推薦直接計(jì)算相對(duì)危險(xiǎn)度(relative risk，RR)，采用logistic回歸計(jì)算的OR值往往會(huì)高估效應(yīng)，尤其在結(jié)局事件比較常見(jiàn)時(shí)[1-3]。2004年Zou[4]提出利用修正Poisson回歸法計(jì)算RR值，該法采用“三明治法”(sandwich variance estimator)矯正標(biāo)準(zhǔn)誤差，從而獲得了比較穩(wěn)健的誤差估計(jì)值，目前這種方法在諸多隊(duì)列研究案例中得到廣泛應(yīng)用[5]。運(yùn)用修正Poisson回歸模型對(duì)服從二項(xiàng)分布的資料進(jìn)行分析能給實(shí)際研究帶來(lái)很多便利，這主要包括：通過(guò)一種穩(wěn)健的誤差方差估計(jì)法(sandwich variance estimator)校正RR值以解決普通Poisson回歸對(duì)RR的估計(jì)誤差較大的問(wèn)題，有實(shí)用的標(biāo)準(zhǔn)軟件(如SAS的GENMOD過(guò)程)等。與標(biāo)準(zhǔn)logistic回歸模型相似，修正Poisson回歸要求觀測(cè)值之間必須是相互獨(dú)立的，否則會(huì)導(dǎo)致統(tǒng)計(jì)推斷的系統(tǒng)性偏差。但在實(shí)際研究中，樣本之間許多都是具有相互關(guān)系的觀測(cè)值，即觀測(cè)數(shù)據(jù)并非來(lái)自完全獨(dú)立的隨機(jī)樣本。傳統(tǒng)的修正Poisson方法無(wú)法應(yīng)對(duì)非獨(dú)立性結(jié)局事件的回歸建模，比如多中心隊(duì)列研究(此類(lèi)結(jié)局事件通常存在聚集性)。因此，本研究基于國(guó)內(nèi)外研究進(jìn)展，并結(jié)合SAS軟件，探討如何運(yùn)用SAS軟件結(jié)合修正Poisson回歸分析非獨(dú)立性數(shù)據(jù)。

1 方法介紹

1.1 普通Poisson回歸

Poisson分布指的是在一個(gè)極大人群、空間和時(shí)間范圍內(nèi)，觀察對(duì)象某種現(xiàn)象發(fā)生數(shù)的分布。常用于稀有事件的發(fā)生次數(shù)的概率分析，以發(fā)生數(shù)作為因變量，構(gòu)建回歸模型，來(lái)探討影響事件發(fā)生的因素。

Poisson回歸是一種將對(duì)數(shù)和二項(xiàng)分布連接起來(lái)的廣義線性模型。該模型可以被寫(xiě)成：

log[λ(Xi)]=β0+β1X1i+β2X2i+...+βKXKi

(1)

在這種情況下，若要計(jì)算變量和因變量的數(shù)量依存關(guān)系，可以基于公式(1)計(jì)算得到

因此，Poisson回歸可以計(jì)算RR值來(lái)反映暴露因素對(duì)結(jié)局事件發(fā)生的影響。

1.1.1 案例1

對(duì)45 例來(lái)自3個(gè)社區(qū)(community)的非器質(zhì)性心臟病且僅有胸悶癥狀就診者進(jìn)行分析研究，以探討吸煙與24 h早搏次數(shù)的關(guān)系。影響因素是X1，是否吸煙(1=吸煙，0=不吸煙)；結(jié)局變量是Y1，24 h早搏次數(shù)(離散型定量數(shù)據(jù)，呈Poisson分布)。其數(shù)據(jù)庫(kù)結(jié)構(gòu)見(jiàn)表1。

表1 案例1、案例2、案例3的數(shù)據(jù)庫(kù)結(jié)構(gòu)

1.1.2 24 h早搏數(shù)的Poisson回歸SAS程序

首先，以24 h早搏數(shù)作為因變量構(gòu)建Poisson回歸模型，自變量包括X1，利用SAS“PROC GENMOD”進(jìn)行模型的構(gòu)建，SAS程序如下：

PROC GENMOD DATA =A;

MODEL Y1 =X1/LINK = log dist = Poisson;

ESTIMATE ′adjusted RR for X1′ X1 1/EXP;

RUN;

1.1.3 結(jié)果分析和解讀

通過(guò)SAS程序，計(jì)算得到X1的RR=1.31，RR的95%CI為(1.04，1.64)，又因?yàn)閄1代表“是否吸煙(1=吸煙，0 =不吸煙)”，說(shuō)明吸煙者出現(xiàn)早搏的風(fēng)險(xiǎn)是不吸煙者的1.31倍。

傳統(tǒng)的Poisson回歸可以廣泛用于呈Poisson分布離散型結(jié)局事件，特別是在偏態(tài)分布的情況下，可以代替線性回歸進(jìn)行數(shù)據(jù)分析。

1.2 修正Poisson回歸

Poisson回歸通常適用于處理罕見(jiàn)結(jié)局事件的前瞻性研究資料，即服從Poisson分布的資料。當(dāng)將其應(yīng)用于服從二項(xiàng)分布的資料時(shí)，對(duì)RR的估計(jì)誤差便會(huì)增大，但是這個(gè)問(wèn)題可以通過(guò)一種穩(wěn)健的誤差估計(jì)法即“三明治法”(sandwich variance estimator)得到校正，被稱(chēng)為修正Poisson回歸，這種方法由Zou[4]在2004年提出。

1.2.1 案例2

同樣基于3個(gè)社區(qū)的非器質(zhì)性心臟病且僅有胸悶癥狀就診者，調(diào)查的暴露因素是X1，是否吸煙(1=吸煙，0=不吸煙)；協(xié)變量是X2，是否喝咖啡(1=喝、0=不喝)；結(jié)局變量是Y2，冠心病是否復(fù)發(fā)(1=復(fù)發(fā)，0=未復(fù)發(fā))。

1.2.2 修正Poisson回歸的SAS程序

修正Poisson回歸法基于廣義估計(jì)方程原理， 利用SAS的GENMOD過(guò)程中 REPEATED 語(yǔ)句估計(jì)得到更為穩(wěn)健的誤差方差，解決了普通Poisson回歸估計(jì)參數(shù)區(qū)間過(guò)于保守的問(wèn)題。SAS程序如下：

PROC GENMOD DATA =A;

CLASS ID;

MODEL Y2 =X1 X2/DIST = Poisson LINK =log;

REPEATED SUBJECT =ID;

ESTIMATE ′adjusted RR for X1′ X1 1/EXP;

RUN;

1.2.3 結(jié)果分析和解讀

GENMOD過(guò)程分析結(jié)果中，X1的RRa=3.06，RRa的95%CI為(1.11，8.49)，又因?yàn)閄1代表“是否吸煙(1=吸煙，0 =不吸煙)”，說(shuō)明調(diào)整混雜因素后吸煙者冠心病復(fù)發(fā)的概率是不吸煙者的3.06倍，見(jiàn)表2?？梢钥闯鲈谂c普通Poisson回歸相同的情況下，得到了更為精確的參數(shù)區(qū)間估計(jì)范圍。

表2 3個(gè)社區(qū)冠心病患者復(fù)發(fā)與吸煙的關(guān)系研究

1.3 針對(duì)非獨(dú)立性數(shù)據(jù)Poisson回歸

在醫(yī)學(xué)研究中，很多事件的發(fā)生是非獨(dú)立性的。例如疾病的聚集性或家族性，或傳染性疾病。修正Poisson回歸是在獨(dú)立數(shù)據(jù)的背景下提出的，并通過(guò)分析和模擬證明在這種情況下是可以使用的[6-8]。廣義估計(jì)方程作為一種證據(jù)性較強(qiáng)的方差估計(jì)方法并考慮了數(shù)據(jù)聚集性，因此可以通過(guò)使用廣義估計(jì)方程來(lái)校正標(biāo)準(zhǔn)誤差，而不是采用通常應(yīng)用于獨(dú)立數(shù)據(jù)的方差估計(jì)方法。

1.3.1 案例3

同樣基于3個(gè)社區(qū)的非器質(zhì)性心臟病且僅有胸悶癥狀就診者，調(diào)查的暴露因素X1，是否吸煙(1=吸煙，0=不吸煙)；協(xié)變量X2，是否喝咖啡(1=喝、0 =不喝)；聚集性變量社區(qū)，community(1、2、3)；結(jié)局變量Y2，冠心病是否復(fù)發(fā)(1=復(fù)發(fā)，0=未復(fù)發(fā))。

1.3.2 非獨(dú)立Poisson回歸的SAS程序

非獨(dú)立Poisson回歸將穩(wěn)健誤差方差估計(jì)法擴(kuò)展應(yīng)用到非獨(dú)立性二分類(lèi)數(shù)據(jù)中，使用穩(wěn)健誤差方差估計(jì)法解釋聚類(lèi)效應(yīng)及Poisson模型作為二分類(lèi)數(shù)據(jù)的工作模型。使用Zou[4]描述修正Poisson的SAS代碼來(lái)完成計(jì)算，其中將SAS PROC GENMOD的重復(fù)語(yǔ)句中單個(gè)體標(biāo)識(shí)符更改為聚集性標(biāo)識(shí)符[9]。SAS程序如下：

PROC GENMOD DATA =A;

CLASS city;

MODEL Y2 =X1 X2/DIST = Poisson LINK =log;

REPEATED SUBJECT =community;

ESTIMATE ′adjusted RR for X1′ X1 1/EXP

RUN;

1.3.3 結(jié)果分析和解讀

在本案例中，X1的RRa=3.06，RRa的95%CI為(1.24，7.58)，又因?yàn)閄1代表“是否吸煙(1=吸煙，0 =不吸煙)”，說(shuō)明吸煙者冠心病復(fù)發(fā)的概率是不吸煙者的3.06倍，見(jiàn)表2?？梢钥闯鲈谙嗤那闆r下與修正Poisson回歸相比，得到了更為精確的參數(shù)區(qū)間估計(jì)范圍。將聚集性變量“社區(qū)”納入回歸模型后，不僅可以改善原本回歸分析中面臨的殘差不獨(dú)立性的問(wèn)題，而且可以進(jìn)一步通過(guò)聚集性變量(社區(qū))減少殘差，提高分析效率。

2 實(shí)例分析

本研究基于中國(guó)健康與養(yǎng)老調(diào)查2011—2018(CHARLS 2011—2018)的數(shù)據(jù)，對(duì)13 283名45～100歲的中老年人進(jìn)行分析研究，以探討我國(guó)中老年人腹型肥胖與死亡的關(guān)系。該研究是多中心的前瞻性隊(duì)列研究，影響因素包括：X1，是否腹型肥胖(1=腹型肥胖，0=非腹型肥胖)；X2，性別(1=男，0=女)；X3，年齡(1=60歲以上，0=45～60歲)；X4，戶(hù)籍(1=非農(nóng)業(yè)，0=農(nóng)業(yè))；X5，婚姻(1=未婚，2=結(jié)婚)。聚集性變量為城市(1=東部城市，2=中部城市，3=西部城市)。

研究分別采用logistic回歸、普通的Poisson回歸、修正Poisson回歸、非獨(dú)立Poisson回歸進(jìn)行分析，調(diào)整混雜因素和中心效應(yīng)后，腹型肥胖死亡的概率是非腹型肥胖者的0.89倍。見(jiàn)表3。

表3 不同的分析方法中老年人腹型肥胖(X1)與死亡的關(guān)系

3 討論

本研究使用具有可交換相關(guān)結(jié)構(gòu)的廣義估計(jì)方程來(lái)解釋聚集性，通過(guò)模擬及實(shí)證非獨(dú)立的前瞻性數(shù)據(jù)，研究非獨(dú)立Poisson回歸方法估計(jì)相對(duì)風(fēng)險(xiǎn)的性能，結(jié)果顯示該方法在少量或大量的集群條件下均表現(xiàn)較好。

與logistic回歸法比較，普通Poisson回歸法適用的資料類(lèi)型范圍更廣，除二分類(lèi)結(jié)局外還可應(yīng)用于處理結(jié)局為離散型定量數(shù)據(jù)。當(dāng)結(jié)局事件的發(fā)生率較為常見(jiàn)(>10%)時(shí)，OR值往往會(huì)明顯高估或低估真實(shí)的RR值，進(jìn)而對(duì)臨床和公共衛(wèi)生的正確決策產(chǎn)生影響[10-11]，這時(shí)直接計(jì)算RR值較為恰當(dāng)。與普通Poisson回歸法比較， 修正Poisson回歸法在參數(shù)點(diǎn)估計(jì)值相同的情況下，得到了更為精確的參數(shù)區(qū)間估計(jì)范圍，從而解決了普通Poisson回歸法對(duì)參數(shù)區(qū)間估計(jì)過(guò)于保守的問(wèn)題。與修正Poisson回歸法比較，非獨(dú)立Poisson回歸法在參數(shù)點(diǎn)估計(jì)值相同的情況下，得到了更為精確的參數(shù)區(qū)間估計(jì)范圍，從而解決了非獨(dú)立性數(shù)據(jù)間具有相關(guān)性的問(wèn)題。Yelland等[12]通過(guò) SAS 軟件模擬數(shù)據(jù)研究也證實(shí)了這一點(diǎn)。

本研究結(jié)果表明，在非獨(dú)立Poisson回歸方法中將廣義估計(jì)方程應(yīng)用于處理非獨(dú)立的前瞻性數(shù)據(jù)是合適的。修正Poisson回歸法作為負(fù)二項(xiàng)回歸的替代方法被提出，用以估計(jì)獨(dú)立數(shù)據(jù)背景下的相對(duì)風(fēng)險(xiǎn)，其中將其原理應(yīng)用在聚類(lèi)數(shù)據(jù)背景下的性能目前才被研究，并通常使用廣義估計(jì)方程解釋聚集性[13-14]。除了使用廣義估計(jì)方程來(lái)解釋數(shù)據(jù)聚集性外，另一種替代方法是擬合具有隨機(jī)聚類(lèi)效應(yīng)的混合模型。區(qū)別于廣義估計(jì)方程，第二種方法必須假設(shè)隨機(jī)效應(yīng)的分布比較難以驗(yàn)證，并且對(duì)其錯(cuò)誤的解釋可能對(duì)結(jié)果產(chǎn)生重大影響。

因此，對(duì)于常見(jiàn)結(jié)局事件的非獨(dú)立前瞻性研究，使用非獨(dú)立Poisson回歸法來(lái)計(jì)算暴露因素的RR值是一種較為簡(jiǎn)單準(zhǔn)確的分析方法，并可利用SAS軟件包中的PROC GENMOD程序來(lái)實(shí)現(xiàn)。