用SAS軟件實現(xiàn)高維列聯(lián)表資料的統(tǒng)計學(xué)分析：因變量為二值變量的多重logistic回歸分析

胡良平，沈 寧，柳偉偉

(軍事醫(yī)學(xué)科學(xué)院研究生部生物醫(yī)學(xué)統(tǒng)計學(xué)咨詢中心，北京 100850)

1 方法原理

Logistic回歸分析屬于概率型回歸分析，其應(yīng)用范圍很廣，不僅適用于流行病學(xué)中病因?qū)W的分析，也可用于臨床療效評價和衛(wèi)生服務(wù)研究等。它適用于因變量是定性變量的案例，包括二值變量、多值有序變量和多值名義變量。Logistic回歸分析對自變量的數(shù)目有一定要求，一般來說，要求樣本含量為自變量個數(shù)的10倍以上，結(jié)果穩(wěn)定性好；對自變量的性質(zhì)沒有特殊要求，自變量可以是連續(xù)性變量、二值變量、多值有序變量或多值名義變量。

程序編寫中，用Y代表因變量，其取值為0和1，分別表示發(fā)生陽性與陰性結(jié)果，X1、X2、…、Xm分別代表m個自變量。設(shè)P(Y=1|X1,X2,…,Xm)表示在自變量X1、X2、…、Xm存在的條件下出現(xiàn)陽性結(jié)果的概率，簡記為P，則多重logistic回歸分析模型可以表示為：

式中β0為常數(shù)項，β1、β2、…、βm分別為各個自變量所對應(yīng)的回歸系數(shù)。陽性結(jié)果發(fā)生的概率為P，則陰性結(jié)果發(fā)生的概率為(1-P)。P與(1-P)之比值被叫做優(yōu)勢比(OR)。對OR取自然對數(shù)，稱為對P作logit變換，用logitP來表示。

由此logistic模型又可以表示為如下的形式：

logitP=β0+β1X1+β2X2+…+βmXm

注意，概率P與自變量X1、X2、…、Xm之間的關(guān)系是非線性的，但是logitP和自變量之間呈線性關(guān)系。Logistic回歸分析模型中各參數(shù)都有明確的實際意義，回歸系數(shù)βi表示在其他影響因素不變的情況下，自變量Xi每變化一個單位時所引起的logitP的改變量，Xi對應(yīng)的ORi=expβi。當某種疾病的發(fā)病率或病死率很低時，可以用OR近似地估計相對風險度(relative risk,RR)。

從整體上理解logistic回歸分析，可以根據(jù)操作過程依次總結(jié)為以下幾個方面：自變量篩選，建立回歸模型，以及檢驗部分(包括回歸系數(shù)的檢驗、整體模型檢驗以及模型擬合優(yōu)度檢驗)。(1)自變量篩選：篩選自變量的方法主要有前進法、后退法、逐步法和最優(yōu)子集法。與多重線性回歸不同的是，logistic回歸分析篩選自變量時所用的檢驗不是F檢驗，而是似然比檢驗、計分檢驗和Wald檢驗。(2)檢驗部分：在logistic回歸分析中，參數(shù)的估計通常采用最大似然法。在估計出回歸系數(shù)后，對全部回歸系數(shù)是否均為0做出的檢驗和對單個回歸系數(shù)的檢驗，常用的檢驗方法有似然比檢驗、計分檢驗和Wald檢驗。(3)模型擬合優(yōu)度檢驗：在建立回歸方程后需要對模型做出評價，考察模型與實際數(shù)據(jù)的符合情況，稱為擬合優(yōu)度檢驗。在logistic回歸分析中，用于擬合優(yōu)度評價的統(tǒng)計量主要包括Pearson χ2值、偏差、Hosmer-Lemeshow統(tǒng)計量和一些信息測量指標。信息測量指標包括Akaike信息準則(AIC)和貝葉斯信息準則(BIC)。這兩個指標在其他眾多模型的評價中都可以看到，其取值越小，表明模型擬合越好。

2 實 例

以晚期非小細胞肺癌患者作為研究對象進行的回顧性研究中， 分析多西他賽聯(lián)合奧沙利鉑與多西他賽聯(lián)合順鉑的療效差異是否具有統(tǒng)計學(xué)意義。 評價指標為接受不同方案治療有效的患者數(shù)， 疾病完全緩解、 部分緩解和穩(wěn)定被定義為有效， 疾病進展和死亡定義為無效。 在考察過程中， 研究者發(fā)現(xiàn)對療效可能有影響的因素除藥物外， 還有年齡 (>55歲或≤55歲)以及患者是否復(fù)治。 表1為原始調(diào)查表格[1]。

表1 兩組患者化療后有效的人數(shù)a

文獻[2]應(yīng)用CMHχ2檢驗對表1中的數(shù)據(jù)進行了分析，結(jié)果顯示，在以年齡和復(fù)治情況為分層因素的基礎(chǔ)上，兩組患者的療效差異有統(tǒng)計學(xué)意義。但在對分層因素(年齡和復(fù)治)進行各層優(yōu)勢比齊性Breslow-Day檢驗時發(fā)現(xiàn)，各分層因素在效應(yīng)方向上不一致，表明運用CMHχ2檢驗獲得的結(jié)果并不可靠。本文在文獻[2]的基礎(chǔ)上，繼續(xù)分析本案例，闡述多重logistic回歸分析的正確應(yīng)用、SAS軟件實現(xiàn)及結(jié)果解釋。

3 分析思路

本研究的分析思路為：首先應(yīng)弄清數(shù)據(jù)結(jié)構(gòu)，資料中涉及的影響因素包括治療方案、年齡、是否為復(fù)治；因變量為治療是否有效，故稱為因變量為二值變量的四維列聯(lián)表資料。處理此種資料可以選擇的統(tǒng)計分析方法有加權(quán)χ2檢驗(限三維列聯(lián)表)、CMHχ2檢驗、logistic回歸分析和對數(shù)線性模型分析。由文獻[2]的分析可見，本案例數(shù)據(jù)采用CMHχ2檢驗進行分析無法得出可靠結(jié)論，故采用二值變量的多重logistic回歸分析為宜。

3.1 原始數(shù)據(jù)記錄表向標準型表格的轉(zhuǎn)換 根據(jù)資料類型以及要選擇的統(tǒng)計分析方法，將原始表1轉(zhuǎn)換為用于統(tǒng)計分析的標準型表格(表2)[2]。在轉(zhuǎn)換表格的過程中，應(yīng)該把同一因素不同水平組合下因變量在各水平下的頻數(shù)全部呈現(xiàn)出來。

3.2 多重logistic回歸分析的SAS實現(xiàn)

3.2.1 程序編寫 SAS程序如下，程序名為LOG1.SAS。

程序說明：自變量a、b、c均為二值變量，因變量Y也是二值變量，故程序中變量賦值均以0和1來代表每個變量的兩種狀態(tài)。a表示年齡分層，a=0表示>55歲，a=1表示≤55歲；b表示復(fù)治與否，b=0表示復(fù)治，b=1表示初治；c表示用藥方案，c=0表示多西他賽聯(lián)合奧沙利鉑，c=1表示多西他賽聯(lián)合順鉑；Y=0表示有效，Y=1表示無效；變量f表示頻數(shù)。數(shù)據(jù)的分析采用logistic過程，在Model語句中依次列出年齡、復(fù)治情況和用藥方案。Model語句后面的“selection=stepwise”表示用逐步法篩選變量。Ods html語句則要求將數(shù)據(jù)處理結(jié)果以網(wǎng)頁格式輸出。SAS在擬合非條件logistic回歸模型時，若不特別指定，計算的結(jié)果是與“Y=0”對應(yīng)的概率表達式Z中的參數(shù)估計值，即應(yīng)當寫成P(Y=0)=[exp(Z)]/[1+exp(Z)]。此處需要注意的是，若自變量中有多值名義變量，則需要引入啞變量來取代多值名義變量，啞變量個數(shù)為該多值名義變量的水平數(shù)減1。

表2因變量為二值變量的四維列聯(lián)表資料的標準型表格

Table2Standardtableofthefour-dimensionalcontingencytablewithabinarydependentvariable

a()>5576 468 2028 1220 20≤5568 1248 328 2012 16

a:對照組治療方案為多西他賽聯(lián)合順鉑，治療組治療方案為多西他賽聯(lián)合奧沙利鉑

SAS程序編寫

3.2.2 結(jié)果輸出及解釋 第一部分：logistic過程產(chǎn)生的第一部分結(jié)果，主要是模型信息和關(guān)于因變量的說明，其中數(shù)據(jù)集名稱為WORK.LOG1，因變量名稱為Y，因變量有兩個水平，權(quán)重變量為w，使用的模型是二值logit模型，參數(shù)估計時的優(yōu)化方法是Fisher’s scoring法。讀入觀測數(shù)和使用觀測數(shù)都是16，讀入和使用的樣本例數(shù)是464例(見表3)。因變量的取值順序為0和1，各自分別有328例和136例(見表4)。最后一行文字說明該模型是以Y=0為基礎(chǔ)，也就是以有效的概率為基礎(chǔ)建模。這里需要注意的是， 在logistic過程中，默認狀態(tài)下是以因變量取值較小的那個水平的發(fā)生概率為基礎(chǔ)建模。本例中使用Y=0代表治療有效，所以是以治療有效為基礎(chǔ)建立概率模型。

表3 Logistic過程中的模型信息

表4 響應(yīng)變量統(tǒng)計結(jié)果

第二部分：對模型擬合優(yōu)度進行檢驗。模型擬合統(tǒng)計量包括AIC、SC和-2倍的對數(shù)似然值，這三個統(tǒng)計量取值越小，表明模型擬合優(yōu)度越好(見表5)。檢驗結(jié)果可見，在包含自變量的模型中，這三種方法計算的統(tǒng)計量取值都小于不包含自變量的模型，表明模型擬合得較好。因此，對回歸模型整體是否成立進行假設(shè)：Beta=0，其原假設(shè)是所有的回歸系數(shù)都為0，分別使用似然比檢驗、計分檢驗和Wald檢驗三種方法。檢驗結(jié)果中依次給出了χ2值、自由度(df值)和P值，3種方法的P值都<0.001，可以認為該模型是成立的(見表6)。

第三部分：在對擬合多重logistic回歸分析方程過程中的變量進行篩選時，一般只保留對因變量有統(tǒng)計學(xué)意義的變量。這部分結(jié)果是回歸系數(shù)以及優(yōu)勢比的估計結(jié)果。逐步回歸法篩選變量的結(jié)果顯示，a、b、c都是有統(tǒng)計學(xué)意義的變量(P<0.000 1)，均被保留在回歸方程中(見表7)。由此可以寫出治療有效的概率表達式為：

表5 模型擬合優(yōu)度檢驗結(jié)果

表6 對回歸模型整體是否成立進行假設(shè)檢驗的結(jié)果

表7 最大似然比檢驗結(jié)果

計算優(yōu)勢比估計值發(fā)現(xiàn)，>55歲的患者接受治療后有效的可能性是≤55歲者的0.391倍(95% CI為0.250～0.611)；復(fù)治的患者接受治療后有效的可能性是初治患者的0.212倍(95% CI 為0.134～0.337)；接受多西他賽聯(lián)合奧沙利鉑治療有效的可能性是多西他賽聯(lián)合順鉑者的0.386倍(95% CI為0.247～0.605)，即多西他賽聯(lián)合順鉑的療效優(yōu)于多西他賽聯(lián)合奧沙利鉑。以下所有OR的區(qū)間均不包含1，提示OR與1之間的差異有統(tǒng)計學(xué)意義。

最后輸出預(yù)測概率和響應(yīng)變量之間的關(guān)聯(lián)性，其中統(tǒng)計量c的取值等于ROC曲線的曲線下面積(AUC)，是反映模型預(yù)測準確性的另外一個指標。計算得到c=0.747，表明該模型的預(yù)測效果較好(見表8)。

4 討 論

近年來，logistic回歸分析在眾多的臨床醫(yī)學(xué)研究中得到了應(yīng)用[3-6]，多數(shù)研究采用logistic回歸分析模型進行預(yù)測分析，呈現(xiàn)出以下幾個問題。(1)缺乏實用性：多數(shù)研究僅篩選出危險因素，卻沒有進一步列出預(yù)測模型，從而使得篩選出的危險因素缺乏利用度和推廣性[3，4]。(2)缺乏模型擬合優(yōu)度檢驗：與馮國雙等[7]分析111篇應(yīng)用logistic回歸分析的文章提出的問題相一致，許多l(xiāng)ogistic回歸分析均缺乏模型擬合優(yōu)度檢驗。評價模型的基本要點是穩(wěn)定性和準確性。在建立回歸方程以后，需要對模型做出評價，即考察模型與實際數(shù)據(jù)的符合情況。除了采用評價擬合優(yōu)度的主要統(tǒng)計量，如Pearson χ2值、偏差、Hosmer-Lemeshow統(tǒng)計量和一些信息測量指標以外，在模型的應(yīng)用后期，可以建立驗證集(validation set)對模型的準確度進行前瞻性的考察[6]。(3)用ROC曲線評價模型擬合優(yōu)度：在SAS 9.2統(tǒng)計軟件中，可以采用“ODS Statistical Graphics”功能繪制ROC曲線。根據(jù)ROC曲線與45度線的偏離程度反映模型預(yù)測效果，其優(yōu)點在于直觀，不足之處是不夠精確。ROC曲線中，45度線下面積為0.5，ROC曲線與其偏離越大，就越向左上方靠攏，它下面的面積(AUC)也應(yīng)該越大。SAS的logistic回歸分析結(jié)果報告的“Association Statistics”中指標c是Wilcoxon-Mann-Whitney檢驗的統(tǒng)計量，等價于ROC曲線的AUC，是反映模型預(yù)測準確度的另外一種方式。

表8 預(yù)測概率與響應(yīng)變量之間的關(guān)聯(lián)性

本文采用logistic回歸分析對患者(分層因素：年齡、是否復(fù)治)接受多西他賽聯(lián)合奧沙利鉑治療和多西他賽聯(lián)合順鉑治療的療效進行了分析，結(jié)果顯示年齡、復(fù)治情況和用藥方案都是對療效有統(tǒng)計學(xué)意義的影響因素，根據(jù)三個因素建立的模型預(yù)測準確度較好。

【參考文獻】

[1] 裘光賢，周兆春，王 婕.多西他賽聯(lián)合奧沙利鉑治療晚期非小細胞肺癌的療效和安全性[J].藥學(xué)服務(wù)與研究,2011,11(3):210-212.

Qiu GuangXian,Zhou ZhaoChun,Wang Jie.Efficacy and safety of docetaxel combined with oxaliplatin in treatment of patients with advanced non-small-cell lung cancer[J].Pharm Care Res,2011,11(3):210-212.In Chinese with English abstract.

[2] 胡良平,沈 寧.用SAS軟件實現(xiàn)高維列聯(lián)表資料的統(tǒng)計學(xué)分析：CMH χ2檢驗[J].藥學(xué)服務(wù)與研究,2013,13(5):332-335.

Hu LiangPing,Shen Ning.Statistical analysis and SAS solutions for multi-dimensional contingency table:CMH χ2test[J].Pharm Care Res,2013,13(5):332-335.In Chinese with English abstract.

[3] 胡曉斌，楊軼男，白亞娜.蘭州市6825例住院孕產(chǎn)婦不良妊娠結(jié)局危險因素的logistic回歸分析[J].中國循證醫(yī)學(xué)雜志,2011,11(5):505-507.

Hu XiaoBin,Yang YiNan,Bai YaNa.Risk factor logistic regression on adverse pregnancy outcomes in 6825 hospitalized pregnant women in Lanzhou city[J].Chin J Evid-based Med,2011,11(5):505-507.In Chinese with English abstract.

[4] 李慶祥,張 瑩,朱小玲,等.心肌梗死患者非酒精性脂肪性肝病的發(fā)生率及危險因素分析[J].中國循證醫(yī)學(xué)雜志,2013,13(4):404-408.

Li QingXiang,Zhang Ying,Zhu XiaoLing,etal.Incidence and risk factors of myocardial infarction complicated by non-alcoholic fatty liver disease[J].Chin J Evid-based Med,2013,13(4):404-408.In Chinese with English abstract.

[5] Kikuchi M,Mikami T,Sato T,etal.High Ki67,Bax,and thymidylate synthase expression well correlates with response to chemoradiation therapy in locally advanced rectal cancers:proposal of a logistic model for prediction[J].Br J Cancer,2009,101(1):116-123.

[6] Roberts E B,Grayson A D,Alahmar A E,etal.Predicting angiographic outcome in contemporary percutaneous coronary intervention:a lesion-specific logistic model[J].J Interv Cardiol,2010,23(4):394-400.

[7] 馮國雙，陳景武，周春蓮.Logistic回歸應(yīng)用中容易忽視的幾個問題[J].中華流行病學(xué)雜志,2004,25(6):544-545.

Feng GuoShuang,Chen JingWu,Zhou ChunLian.Several problems easily ignored in applications of logistic regression [J].Chin J Epidemiol,2004,25(6):544-545.In Chinese.