N∶M條件logistic回歸分析在統計軟件上的實現*

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

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

問題的提出

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

原理及方法

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

2．Cox比例風險模型

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

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

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

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

有:

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

實例分析

1．資料概述

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

表1 變量說明

2．條件logistic模型的軟件實現

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

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

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

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

library(splines)

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

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

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

表2 clogit模塊結果整理

3．Cox比例風險模型的軟件實現

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

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

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

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

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

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

(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;

*［在此輸入數據］

Time=2－Low;

run;

proc phreg data=LBW;

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

strata Age;

run;

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

表3 各軟件Cox模塊結果整理

4．SPSS17．0 軟件的誤用

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

表4 SPSS軟件Cox Regression模塊結果整理

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

討 論

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

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

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

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

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

5．張業武．Cox比例風險模型對條件logistic回歸參數估計原理和方法．中國衛生統計，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．