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

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

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

1 方法介紹

1.1 普通Poisson回歸

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

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

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

(1)

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

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

1.1.1 案例1

對45 例來自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ù)庫結(jié)構(gòu)見表1。

表1 案例1、案例2、案例3的數(shù)據(jù)庫結(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é)果分析和解讀

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

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

1.2 修正Poisson回歸

Poisson回歸通常適用于處理罕見結(jié)局事件的前瞻性研究資料，即服從Poisson分布的資料。當(dāng)將其應(yīng)用于服從二項(xiàng)分布的資料時(shí)，對RR的估計(jì)誤差便會(huì)增大，但是這個(gè)問題可以通過一種穩(wěn)健的誤差估計(jì)法即“三明治法”(sandwich variance estimator)得到校正，被稱為修正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過程中 REPEATED 語句估計(jì)得到更為穩(wěn)健的誤差方差，解決了普通Poisson回歸估計(jì)參數(shù)區(qū)間過于保守的問題。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過程分析結(jié)果中，X1的RRa=3.06，RRa的95%CI為(1.11，8.49)，又因?yàn)閄1代表“是否吸煙(1=吸煙，0 =不吸煙)”，說明調(diào)整混雜因素后吸煙者冠心病復(fù)發(fā)的概率是不吸煙者的3.06倍，見表2?？梢钥闯鲈谂c普通Poisson回歸相同的情況下，得到了更為精確的參數(shù)區(qū)間估計(jì)范圍。

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

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

在醫(yī)學(xué)研究中，很多事件的發(fā)生是非獨(dú)立性的。例如疾病的聚集性或家族性，或傳染性疾病。修正Poisson回歸是在獨(dú)立數(shù)據(jù)的背景下提出的，并通過分析和模擬證明在這種情況下是可以使用的[6-8]。廣義估計(jì)方程作為一種證據(jù)性較強(qiáng)的方差估計(jì)方法并考慮了數(shù)據(jù)聚集性，因此可以通過使用廣義估計(jì)方程來校正標(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ú)立性二分類數(shù)據(jù)中，使用穩(wěn)健誤差方差估計(jì)法解釋聚類效應(yīng)及Poisson模型作為二分類數(shù)據(jù)的工作模型。使用Zou[4]描述修正Poisson的SAS代碼來完成計(jì)算，其中將SAS PROC GENMOD的重復(fù)語句中單個(gè)體標(biāo)識符更改為聚集性標(biāo)識符[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 =不吸煙)”，說明吸煙者冠心病復(fù)發(fā)的概率是不吸煙者的3.06倍，見表2?？梢钥闯鲈谙嗤那闆r下與修正Poisson回歸相比，得到了更為精確的參數(shù)區(qū)間估計(jì)范圍。將聚集性變量“社區(qū)”納入回歸模型后，不僅可以改善原本回歸分析中面臨的殘差不獨(dú)立性的問題，而且可以進(jìn)一步通過聚集性變量(社區(qū))減少殘差，提高分析效率。

2 實(shí)例分析

本研究基于中國健康與養(yǎng)老調(diào)查2011—2018(CHARLS 2011—2018)的數(shù)據(jù)，對13 283名45～100歲的中老年人進(jìn)行分析研究，以探討我國中老年人腹型肥胖與死亡的關(guān)系。該研究是多中心的前瞻性隊(duì)列研究，影響因素包括：X1，是否腹型肥胖(1=腹型肥胖，0=非腹型肥胖)；X2，性別(1=男，0=女)；X3，年齡(1=60歲以上，0=45～60歲)；X4，戶籍(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倍。見表3。

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

3 討論

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

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

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

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