N∶M條件logistic回歸分析在統(tǒng)計軟件上的實現(xiàn)*

廣州南方醫(yī)科大學公共衛(wèi)生與熱帶醫(yī)學學院生物統(tǒng)計學系(510515) 鄭一男 曹佩華 歐春泉

本文在STATA、R、SAS常用統(tǒng)計軟件上實現(xiàn)N∶M條件logistic回歸分析，并指出Cox風險比例模型進行條件logistic回歸分析時結(jié)點調(diào)整方法的適用情況。

問題的提出

在有關(guān)疾病危險因素的流行病學調(diào)查研究中，常將病例與對照按年齡、性別等因素進行匹配，以消除某些混雜因素對研究結(jié)果的干擾作用〔1－3〕。大家熟悉的是1∶1或1∶M匹配，即給每個病例選擇1個或M個對照與之匹配。匹配后往往再按匹配因素進行分層分析，這樣可使每一個匹配層中都有一定數(shù)目的病例與對照?；蛘咧苯硬捎肗∶M匹配，每一匹配層中有N個病例和M個對照，而且各層中病例數(shù)和對照數(shù)的比例可以不固定，研究實施更為靈活。配對的病例-對照研究數(shù)據(jù)均采用條件logistic回歸分析。關(guān)于1∶1和1∶M條件logistic回歸分析的軟件應用有不少討論，但對N∶M條件logistic回歸分析卻鮮有介紹，筆者在實際工作中遇到不少非統(tǒng)計專業(yè)的研究人員錯誤地運用SPSS軟件的Cox模塊進行N∶M條件logistic回歸分析，導致結(jié)果出現(xiàn)較大偏差。本文將全面介紹STATA、R和SAS統(tǒng)計軟件如何對N∶M配對病例-對照研究資料進行條件logistic回歸分析，并詳細闡述在使用Cox模塊進行條件logistic回歸分析時的注意事項。讀者可根據(jù)自身的使用偏好選擇任一軟件進行分析。

原理及方法

*:廣東省科技計劃項目(項目號:0903109)△通訊作者:歐春泉，E-mail:ocq@fimmu．com

2．Cox比例風險模型

運用Cox比例風險模型可實現(xiàn)條件logistic回歸分析〔5〕。由于同一層內(nèi)的生存時間均相同，故生存時間的結(jié)點(ties)問題非常突出，如何處理結(jié)點是分析的一個關(guān)鍵。目前主流統(tǒng)計軟件中有四種處理方法:Breslow、Efron、Exact marginal likelihood(精確邊際似然法)與 Exact partial likelihood(精確偏似然法)。Breslow是最簡單快捷的，Efron比Breslow更加精確，但計算速度略慢，兩種方法在處理大量結(jié)點數(shù)據(jù)時均會使參數(shù)估計結(jié)果出現(xiàn)嚴重偏差。精確邊際似然法實際上是利用15-point Gauss-Laguerre quadrature近似計算得到，對于結(jié)點數(shù)較少的數(shù)據(jù)，該法比Efron慢;隨著結(jié)點數(shù)增加，其效率優(yōu)勢能迅速顯露出來，但當平均結(jié)點數(shù)大于30時〔6〕，誤差將逐漸加大。對于生存時間為離散變量的資料，由于產(chǎn)生的結(jié)點較多，應使用精確偏似然法進行處理。

在1∶1和1∶M資料中，每層中僅有一個結(jié)點，因此四種結(jié)點處理方法均等價可行。對于N∶M資料(N為病例數(shù))，每層中有N個結(jié)點，使用精確偏似然法最為合適。

下面給出Breslow與精確偏似然法的對數(shù)似然函數(shù) ln L〔6〕:

在時間區(qū)間(t0i，ti］中的第 i個觀測(i=1，…，N)

有:

j為死亡時間 t(j)的順序編號(j=1，…，D);dj，rj分別為t(j)時刻的死亡數(shù)與存活數(shù);Dj，Rj分別為在t(j)時刻死亡病人的集合與存活病人的集合(危險集);k為觀測編號使得t0k＜t(j)≤tk;δij為0-1變量，表示病人死亡(=1)或截尾(=0)。在病例-對照研究中，上述“死亡”與“存活”的概念分別對應“病例”與“對照”。比照似然函數(shù)(2)和(4)可知，使用精確偏似然法調(diào)整結(jié)點的Cox比例風險模型與條件logistic回歸模型完全等價。

實例分析

1．資料概述

為研究嬰兒出現(xiàn)低出生體重(體重≤2500g)的危險因素，Hosmer和Lemeshow按照母親年齡采取配對的病例-對照研究〔7〕，調(diào)查共涉及189名女性，其中59名產(chǎn)有低出生體重嬰兒，130名所產(chǎn)的嬰兒體重正常。數(shù)據(jù)文件可在南方醫(yī)科大學生物統(tǒng)計學系網(wǎng)站(http://www．echobelt．org)下載。

表1 變量說明

2．條件logistic模型的軟件實現(xiàn)

STATA 9．0和R語言均提供了條件logistic分析模塊，即clogit模塊。

(1)STATA 9．0對實例資料的分析過程為:

clogit low lwt smoke ht ut，group(age)/* 輸出參數(shù)估計值及置信域*/clogit，or/*輸出OR值及置信域*/

(2)R2．11．1對實例資料的分析過程為:

library(splines)

library(survival)#載入自帶軟件包鏡像

clogit(Low～LWT+Smoke+HT+UT+strata(Age))

STATA和R軟件的logit模塊分析結(jié)果完全相同，整理結(jié)果于表2。表中最后兩列展示的是比數(shù)比(Odds Ratio，OR)的95%置信區(qū)間。

表2 clogit模塊結(jié)果整理

3．Cox比例風險模型的軟件實現(xiàn)

各軟件處理結(jié)點的方法大體相同，但需要注意區(qū)分各軟件命名與含義的差異:STATA中有breslow(默認)、efron、exactm(精確邊際似然法)與exactp(精確偏似然法);R中“method”選項里包括 efron(默認)、breslow與exact(精確偏似然法);SAS中“ties”選項里包括breslow(默認)、efron、exact(精確邊際似然法)與discrete(精確偏似然法);SPSS沒有結(jié)點處理的可選項，運算過程使用內(nèi)置的breslow法。

使用Cox比例風險模型前，需在原始數(shù)據(jù)基礎(chǔ)上增加生存時間變量Time=2－Low(注:對照組生存時間大于病例組即可)。

(1)STATA 9．0的Cox比例風險模型為stcox模塊，實例資料的分析過程為:

stset time，failure(low)/*定義時間變量與反應變量*/

stcox lwt smoke ht ut，strata(age)exactp/* 輸出HR值及置信域*/

stcox，nohr/*輸出參數(shù)估計值及置信域*/

(2)R2．11．1的 Cox比例風險模型為 coxph模塊，實例資料的分析過程為:

coxph(Surv(Time，Low)～ LWT+Smoke+HT+UT+strata(Age)，method="exact")

(3)SAS9．2對于1∶1資料的條件logistic回歸可以使用logistic過程(proc logistic)來完成，但對于1∶M和N∶M資料，就必須使用Cox比例風險模型(proc phreg)，實例資料的分析過程為:

data LBW;

input Age Low LWT Smoke HT UT;

*［在此輸入數(shù)據(jù)］

Time=2－Low;

run;

proc phreg data=LBW;

model Time*Low(0)=LWT Smoke HT UT/ties=discrete;

strata Age;

run;

上述三種方法的結(jié)果匯總于表3，可見與表2中條件logistic回歸分析結(jié)果完全一致，但輸出結(jié)果通常稱為風險比(Hazard Ratio，HR)，而非比數(shù)比。表中展示的是風險比的95%置信區(qū)間。

表3 各軟件Cox模塊結(jié)果整理

4．SPSS17．0 軟件的誤用

在SPSS17．0軟件的菜單操作中沒有提供處理條件logistic回歸的模塊，如果使用SPSS軟件的Cox模塊對實例資料進行分析，結(jié)果如表4:

表4 SPSS軟件Cox Regression模塊結(jié)果整理

可以看出上述估計值低估了真實結(jié)果，原因就在于結(jié)點的處理使用的是Breslow法。因此SPSS軟件目前在菜單操作中尚無法實現(xiàn)N∶M資料的條件logistic回歸分析。但對于1∶1與1∶M的資料是可行的。

討 論

STATA軟件與R軟件中均有條件logistic回歸的專用模塊(clogit)，非常簡便靈活;SAS軟件中沒有類似的專用模塊，而是使用與之等價的Cox風險比例回歸模塊(phreg)進行分析。使用前需要在原數(shù)據(jù)基礎(chǔ)上加入生存時間變量(time)，并且使用精確偏似然法處理結(jié)點，這在STATA中的stcox模塊與R軟件中coxph模塊也得到了印證。1∶1和1∶M配對設(shè)計是N∶M配對的特例，故本文介紹的N∶M條件logistic回歸分析的方法同樣適合于1∶1或1∶M條件logistic回歸。

1．方亞，胡海蘭．女性乳腺癌危險因素及其變化．中國衛(wèi)生統(tǒng)計，2009，26(3):241-246．

2．劉靜，王潔貞，薛付忠，等．腎綜合征出血熱發(fā)病率與氣象因素關(guān)系的研究．中國衛(wèi)生統(tǒng)計，2006，23(4):326-329．

3．陸云霞，施侶元，余紅平，等．武漢市居民飲食模式與食管癌發(fā)病關(guān)系的條件logistic回歸分析．中國衛(wèi)生統(tǒng)計，2005，22(3):146-148．

4．STATA Base Reference Manual Volume 1 Release 10．Texas，USA:STATA，2007，286．

5．張業(yè)武．Cox比例風險模型對條件logistic回歸參數(shù)估計原理和方法．中國衛(wèi)生統(tǒng)計，2002，19(1):23-25．

6．STATA Survival Analysis and Epidemiological Table Reference Manual Release 10．Texas，USA:STATA，2007，153．

7．Hosmer D，Lemeshow S．Applied Logistic Regression．2nd edit．New York，USA:John Wiley ＆ Sons，1989，319．