SPSS學習筆記：因變數二分類資料的logistic回歸分析

03-04

轉自個人微信公眾號【Memo_Cleon】的統計學習筆記

http://weixin.qq.com/r/wjv27irEQHc-rW6N927Q (二維碼自動識別)

兩個概念：RR和OR
二分類資料的logistic回歸SPSS操作示例
幾個需要注意的問題：樣本量、啞變數、模型擬合效果和擬合優度檢驗、多重共線

【1】兩個概念

RR（Relative Risk）：相對危險度，也稱危險比（Risk Ratio）或率比（Rate Ratio），在前瞻性研究中用以表示暴露與疾病發生的關聯強度，說明暴露組發病危險是非暴露組發病危險的多少倍，是兩組發病率之比，計算公式為：

OR（Odds Ration）：比值比，優勢比，比數比。在病例對照研究中說明暴露與疾病的關聯強度，它是暴露在某危險因子下的發病率P1與不發病率之比（1-P1）與未暴露在某危險因子下的發病率P2與不發病率（1-P2）之比的比值。

比值（odd）表示發生的可能性與不發行的可能性之比，odds1=P1/（1-P1），odds2=P2/（1-P2）。OR=odds1/odds2=ad/bc

在logistic回歸中，各自變數的回歸係數βi表示自變數Xi每改變一個單位，比值比的自然對數值該變數，而exp(βi)即OR，表示自變數Xi每改變一個單位，陽性結果出現概率與不出現概率之比是變化前相應比值的倍數，即優勢比。當陽性結果出現概率較小時（一般小於0.05）或者較大時（大於0.95），OR=（P1/（1-P1））/（P2/(1-P2)）≈P1/P2=RR

【2】二分類資料的logistic回歸SPSS操作示例

適用條件：①因變數為二分類變數，自變數可以是連續變數也可以是分類變數；②各觀測間相互獨立；③自變數與因變數logit（P）之間存在線性關係；④自變數間不存在多重共線；⑤盡量避免異常值，殘差服從二項分布合計為零；⑥研究對象無論是病例組還是對照組，樣本量至少是需要分析的自變數個數的10倍，多分類自變數的啞變數的參照水平頻數至少為30。

示例：探討冠心病發生的有關危險因素。分析數據來自孫振球主編的《醫學統計學》第三版。

（1）數據錄入賦值。二分類變數取值只有「0」和「1」兩個值，本例年齡和體重指數嚴格來講是要設成啞變數的（可參見本文後面關於啞變數的介紹），但由於樣本量的問題，為更好的演示spss操作，本例按等級資料直接賦值。

（2）分析（Analyze）>>回歸（Regression）>>二分類logistic（Binary Logistic…）

因變數（Dependent）：選入冠心病[Y]

協變數（Covariates）：選入自變數[X1]-[X8]

篩選變數方法（Method）：

進入（Enter）：強行引入，全部變數一次進入模型，不進行篩選
向前：條件（Forward:Conditional）：向前逐步選擇法，剔除變數的標準是條件參數估計的似然比檢驗
向前：似然比（Forward:LR）：向前逐步選擇法，剔除變數的標準是最大偏似然估計的似然比檢驗
向前：Wald（Forward:Wald）：向前逐步選擇法，剔除變數的標準Wald卡方檢驗
不同的向前法選入自變數時均採用比分檢驗，只是採用的剔除變數的檢驗方法不同
向後法也三種，剔除標準分別對應條件、似然比和Wald

本例選中Forward:LR。

分類變數（Categorial…）：可將多分類自變數設置為啞變數。

保存（Save…）：可將中間的計算結果存儲起來供以後分析，共有預測值、殘差和影響強度因子三類

選項（Options…）：可以選中相應選項對模型進行描述、預測和診斷，設置逐步法的進入和剔除標準等。本例選中CI for exp(B)選項

結果解釋：

（1）數據處理情況匯總，本例共有54條記錄納入分析，無缺失值。本例樣本量時偏少，要分析8個自變數至少需要160例記錄，本例僅用於操作演示。

（2）因變數分類賦值

（3）預測分類表：此步開始進行擬合。Block 0擬合的是不包含任何自變數只有常數的無效模型，本例有28例對照的記錄被預測為對照，有26例病例也被預測為對照，總預測準確率為51.9%，這是不納入任何解釋變數時的預測準確率，相當於比較基線。

Block 0: Beginning Block

（4）納入方程的變數及檢驗情況。因此步擬合僅有常數項，β0=-0.074；檢驗方法Wald檢驗，檢驗統計量Wald χ2=0.074，P=0.786>0.05，H0假設是回歸係數為0。Exp(B)為自然對數e的β0次方，實際意義是總體研究對象中（54例）中患病率（26/54）與未患病率（28/54）的比值0.929。

（5）尚未納入模型方程的變數及其比分檢驗結果。所作的檢驗是分別將他們納入方程，方程的改變是否有統計學意義。本例如果將X1-X8分別納入方程，則方程的改變都是有統計學意義的，總的統計量也有統計學意義。逐步回歸法（Stepwise）是一個一個的選入變數，下一步將會先納入P值最小的變數X6重新計算選擇。

（6）模型係數的綜合檢驗，此步開始Block 1的擬合，本例採用的方法為向前逐步選擇法（似然比），本例依次引入了4個變數，結果顯示每一步新引入的變數和最終的模型均有統計學意義（新引入變數係數不為零，所有引入變數係數不全為零）

Block 1: Method = Forward Stepwise (Likelihood Ratio)

此處輸出的即為每一步的似然比檢驗結果，模型新引入一個或幾個變數後形成的新模型（當前模型）與「上一步（step）的模型/上一個區塊（block）模型/初始模型（model）」進行比較，卡方值是上一個步模型/上一個區組模型/初始模型的-2logLR（負2倍的似然比值）與當前的-2logLR的差值，可通過下一張表格（Model Summary）中的值驗證，step的卡方值結果上一步的模型與納入新變數後的模型-2logLR之差，model的卡方值是初始模型（step0，各自變數係數均為0）與納入新變數後的模型-2logLR之差，本例只設置了一個block，因此block與model相同。step檢驗的是每個step被納入的變數引起的似然比變化，model則提示納入新變數後整個模型的似然比變化。

SPSS軟體中自帶的幫助對幾個卡方的描述如下：

Model Chi-Square

2(log-likelihood function for current model ? log-likelihood function for initial model)。 The initial model contains a constant if it is in the model; otherwise, the model has no terms. The degrees of freedom for the model chi-square statistic is equal to the difference between the numbers of parameters estimated in each of the two models. If the degrees of freedom is zero, the model chi-square is not computed.

Block Chi-Square

2(log-likelihood function for current model ? log-likelihood function for the final model from the previous method)。The degrees of freedom for the block chi-square statistic is equal to the difference between the numbers of parameters estimated in each of the two models.

Improvement Chi-Square

2(log-likelihood function for current model ? log-likelihood function for the model from the last step)。The degrees of freedom for the improvement chi-square statistic is equal to the difference between the numbers of parameters estimated in each of the two models.

（7）模型概況，輸出每一步的-2logLR（-2倍的對數似然比值）以及兩個廣義決定係數（也稱偽決定係數），這兩個指標是模型擬合效果的的判斷指標。logistic模型估計一般採用最大似然法，即是得到模型的似然函數值LR達到最大值（LR取值在0-1之間）。模型預測效果越好，LR越大，-2logLR越小。Cox&Snell R2和Nagelkerke R2表示當前模型中的自變數導致的因變數變異占因變數總變異的比例。本例在擬合過的4個模型中-2logLR逐步減小，廣義決定係數在增大，模型擬合效果還是不錯的。

（8）每一步的預測分類表，預測準確率由Block 0的51.9%上升至77.8%。逐步回歸得出的結果是保證模型的最大似然函數值最大，但不能保證此時的預測精度最高。

每一步第一行最後的百分比為真陰性率（特異度），第二行為真陽性率（靈敏度），以第四步結果為例，特異度75.0%，靈敏度80.8%。

（9）每一步擬合入選方程的變數檢驗情況，輸出每一步模型中參數估計值，包括常數項、各自變數的回歸係數及標準誤、Wald卡方值、自由度、P值及OR值（Exp(B)）。最後一步的入選變數作為最終結果，最終篩選出的危險因素有4個，分別是年齡(X1)、高血脂史(X5)、動物脂肪攝入(X6)、A型性格(X8)，說明年齡、高血脂史、動物脂肪攝入及A型性格與冠心病呈正相關，優勢比分別是2.519、4.464、23、7.008。以A型性格為例，其結果可做如下解釋：不考慮其他因素的印象，A型性格的對非A型性格的人發生冠心病的優勢比是7.008（一般文獻中會近似的表達為：不考慮其他因素的印象，A型性格的人發生冠心病的可能是非A型性格的人發生冠心病可能的7.008倍）。對於本例年齡和體重指數，按統計結果直接表達年齡每增加10歲/體重指數每增加一級引起的優勢比，但此解釋應當慎重，因為我們不能確定年齡和體重指數對冠心病的影響和我們的賦值水平是一致的，嚴格來說此處兩變數應該設置成啞變數進行分析，但由於樣本量的原因本例僅做演示。啞變數的設置可參見本文後面關於啞變數的設置。

logit(P)=-4.705+0.924X1+1.496X5+3.135X6+1.947X8

需要特別指出的是，本例X1的Wald檢驗P=0.053>0.05依然保留在了最終的模型中，是因為未達到變數剔除的標準，剔除檢驗見下一張表格：Model if Term Removed。SPSS默認的進入標準是0.05，剔除標準是0.10（可在選項（Options…）中設置），每一步P值小於等於0.1的變數均保留在方程中。

在對尚未進入模型的自變數進行分析時（結果見表(11):Variables not in the Equation），在進行到step3引入X6後，X1的比分檢驗結果,P=0.042<0.05，符合納入標準，因此將其納入方程。在進行變數剔除時（Model if Term Removed）採用的似然比檢驗，P=0.041<0.05，尚未達到剔除標準，因此保留在模型中。從檢驗方法上看，似然比檢驗是基於整個模型的擬合情況，比分檢驗一般與似然比檢驗類似，而Wald檢驗沒有考慮各因素的綜合作用，結果最不可靠，因此當結果出現衝突時應以似然比檢驗結果為準。

（10）輸出每一步被納入的自變數是否需要被剔除，採用的方法是似然比檢驗。結果顯示每一步單獨移除納入的自變數後模型的改變均有統計學意義，因此需要全部保留在方程中。採用不同的變數篩選方法（剔除方法），該處的統計量會有差異。

（11）尚未進入模型的自變數是否可能被納入的比分檢驗結果：尚未進入方程的自變數如果再進入現有模型，模型的改變有無統計學意義。可見在Step 1時，還有多個變數可以引入，其中X5的P值最小，Step2時首先引入；在Step 2時首先考慮引入X8，Step3時首選考慮引入X1，而Step4時在現有的進入標準下已經沒有需要引入的變數了。

需要特別說明的是，變數是否最終入選，最重要的準則是專業判斷，其次為多變數分析，最差的是單變數分析。如在專業上認為某自變數對結果有影響，不論其檢驗水準如何，都可強行納入方程。多變數模型綜合考慮了各個自變數的影響，結果更為客觀，如多變數模型與單變數模型出現矛盾時應以多變數分析結果為準。

強行納入方法：第一種方法是將欲強行納入的變數和已篩選出的變數重新進行logistic回歸，變數帥選辦法為Enter。第二種方法是利用區組（Block）來進行。將欲強行納入的變數歸為一個區組（block），在篩選方法中選擇「ENTER」來確保這些變數被納入，其他變數則歸入另一區組(block)，篩選方法可以是其他方法，比如Forward:LR。

【3】幾個需要注意的問題

（1）樣本量

越多越好。理論上講，每個自變數所需要的事件數至少為10個，也就是研究對象中較少組的記錄數與自變數個數之商大於等於10。本例有8個自變數需要分析，病例組和對照組至少各需要80例以上，總樣本160例以上，如設置啞變數，自變數增多樣本量要求也就更大，而且參照水平至少30例以。本例對照組只有28例，病例組只有26例，最多只能分析的自變數個數是26/10=2.6個，若要分析所有的自變數，樣本量就需要增加，否則就需要減少自變數個數。本例僅演示操作。

（2）啞變數

二分類變數取值只有「0」和「1」兩個值，它們分別代表兩種成對出現的邏輯概念，如「是」和「否」、「有」和「無」、「高」和「低」等，spss默認取值水平高的為陽性結果。當自變數為無序多分類時，SPPS中對變數的賦值大小並不代表自變數間存在次序或者程度的差異，比如研究血型對性格的影響，對A、B、AB、O血型分別賦值1、2、3、4，並不代表性格的logit（P）會按此順序進行線性遞增或遞減，此時需要設置啞變數，自變數有n個水平，需要n-1個啞變數，如上面的4個血型，則需要3個啞變數來設置。

無序多分類變數可直接設置成啞變數，對於有序多分類，則視情況而定。張文彤老師給出的解決方案是專業上不能判斷不同等級的自變數對因變數的影響程度是一致的時候，需要將有序變數分別以啞變數和連續性變數引入模型，對兩個模型進行似然比檢驗（似然比卡方值等於兩模型-2logLR之差，自由度為兩模型自變數個數之差），如果似然比檢驗無統計學意義，且各啞變數回歸係數間存在相應的等級關係，則可將該自變數作為連續性變數引入模型，否則還是採用啞變數方式引入。本例年齡和體重指數均為有序多分類變數，我們不能確定冠心病的發病率logit轉換值會與設定的年齡等級遞增相同，體重指數不同分類之間劃分切點是不等距的，直接按等距賦值為1、2、3也不太符合實際情況。上面的示例將年齡和體重指數直接引入模型進行篩選是不恰當的，嚴格來說這兩個變數是應該採用啞變數來分析的，但本來例數太少，引入啞變數後所需樣本量應該更多，如果引入啞變數分析可能會使結果誤差很大。實際上本例直接引入後，就出現了年齡不再是冠心病發病的危險因素。為更好的演示啞變數的設置及解讀，將原數據複製一倍，變數只保留X2、X4、X6、X7，重新進行回歸分析。

步驟如下：

分析>>回歸>>二分類logistic…

因變數：選入冠心病[Y]

協變數：選入自變數[X2]、[X4]、[X6]、[X7]

篩選變數方法（Method）：Forward:LR

分類變數（Categorial…）：可將多分類自變數設置為啞變數，將「體重指數(BMI)[X7]」選入分類協變數（Categorial Covariats）框中，低BMI水平（BMI＜24）設為參照水平，將Contrast-Indicator，Reference Category-First，Change，Continue

大部分結果解讀同示例結果，以下僅結合啞變數對應的結果進行解讀。

分類變數的啞變數編碼見下表

擬合納入方程的變數及其檢驗情況見下表。最後一步的入選變數作為最終結果，最終篩選出的危險因素有3個，分別是X2、X6、X7。對於X2的解讀：不考慮其他因素的影響，有高血壓病史的患者發生冠心病的可能是無高血壓史患者發生冠心病可能的4.646倍（近似，實際是OR=4.646），X6的解釋與此類同，高動物脂肪攝入對低脂肪攝入發生冠心病的優勢比是22.042。對啞變數的解釋則是與參照水平相比的結果：不考慮其他因素的影響，BMI24-26水平對BMI＜24水平的優勢比是3.279，BMI＞26水平對BMI＜24水平的優勢比是4.039。

要特別注意，啞變數需要遵循同進同出的原則，即在一個模型中同一個多分類變數的所有啞變數要麼全部納入模型，要麼全部不納入模型。目前在各很多統計論壇上，普遍認為結果的第一行是該自變數的總體檢驗，總體檢驗有統計學意義，所有啞變數都應保留，本例X7總P值=0.04<0.05，X7(1):P=0.04<0.05，X7(2):P=0.53>0.05。

一般來說第一行的總體檢驗有統計學意義，其後的啞變數至少有一個是有統計學意義的，有時可能出現後面的啞變數都沒有統計學意義的情況，因為啞變數選取不同的參照水平，其他啞變數與之相比的P值會不同，出現這種情況可能是其他啞變數與選取的參照水平相比剛好沒有統計學意義。

還有一種情況是此表啞變數總檢驗P>0.05，結果仍保留在方程中，原因可能是不同的檢驗方法造成的，跟正文結果（9）後面的解釋一樣，納入標準、剔除標準以及方程中變數的檢驗方法是不一樣的，以本例為例，納入變數方法是比分檢驗（默認納入標準P=0.05），剔除標準是最大似然比檢驗（默認標準P=0.10），而本表方程中變數的檢驗採用的是Wald卡方檢驗。Wald與似然比檢驗出現矛盾時以似然比檢驗為準。

實際上，自變數設為啞變數後，Variables in the Equations表中直接顯示的都是最終入選的變數，具體的納入和剔除檢驗則需要從Variables not in the Equation和Model if Term Removed表中查看。如本例，結合上表X7在step3時被納入，從模型係數的綜合檢驗表格（Omnibus Tests of Model Coefficient）中可以看出，Step3與上一步的似然比檢驗有統計學意義，納入該變數時模型的改變有統計學意義（Variables not in the Equation表中step2比分檢驗結果：P=0.033<0.05），剔除該變數時模型的改變也有統計學意義（Model if Term Removed表中step3的似然比檢驗結果：P=0.034<0.05）。

（3）模型效果檢驗

模型擬合效果可以通過對似然比值的變化和廣義決定係數（Model Summary）、預測分類表（Classification Table），解讀如正文。另外對模型效果的檢驗還有ROC曲線。

製作ROC曲線步驟如下：

分析>>回歸>>二分類logistic…

因變數：選入冠心病[Y]
協變數：選入自變數[X2]、[X4]、[X6]、[X7]
篩選變數方法（Method）：Forward:LR
保存（Save…）：選中預測值中的概率(Probabilities)
Continue
OK

分析>>ROC曲線（ROC Curve…）

檢驗變數（Test Variable）：選入新生產的預測概率（Predicted probability[PRE_1]）
狀態變數（State Variable）：選入冠心病[Y]，將Y=1的研究對象指定為病人（狀態變數值(Value of State Variable)：1）。
顯示（Dispaly）選中所有複選框：ROC曲線、帶對焦參考線、標準誤和置信區間、ROC曲線的坐標點
OK

結果顯示曲線下面積為0.883，P<0.001。P值的H0假設是曲線下面積為0.5。

ROC曲線下面積取值範圍0.5-1，0.5-0.7表示診斷價值較低，0.7-0.9表示診斷價值中等，0.9以上表示診斷價值較高。如果ROC曲線沿著對角線方向分布. 表示分類是機遇造成的，正確和錯誤的概率各位50%，此時該診斷方法無效，較好的診斷方法的ROC曲線應該是從左下角垂直上升至頂線.然後水平向右至右上角，更外面離對角線更遠的曲線其靈敏度和特異度均高於裡面的離對角線更近的曲線。

（4）模型的擬合優度的檢驗

可採用Pearson擬合優度檢驗和Deviance擬合優度檢驗，此兩法在SPSS的二分類的logistic回歸沒有輸出

似然比檢驗是計算經篩選出的變數的飽和模型的-2logLR值和只引入篩選出的變數的主效應模型（簡單模型）的-2logLR值，兩者之差服從卡方分布，自由度為兩者參數個數之差，若P>0.05則說明擬合篩選的自變數主效應的logistic回歸模型已經足夠。註：飽和模型可以理解為納入各自變數主效應及其交互項的模型，簡約模型為非飽和模型，比如缺少交互項。

除了似然比檢驗外，SPSS中還有Hosmer-lemeshow檢驗，可在選項（Options…）中選入Hosmer-lemeshow goodness-of-fit，結果出現在「模型概況（Model Summary）」表格之後，輸出「Hosmer和Lemeshow檢驗」以及「Hosmer和Lemeshow檢驗列聯表」。其中Hosmer和Lemeshow檢驗列聯表是十分位組的觀測值和期望值，按模型求得某個體的得病概率，概率遞增排序後分為10等份，據此期望值與相應的觀測值求得Pearson卡方（自由度=組數-2=8，但有時自變數組合和樣本量的原因組數可能少於10）。

本例最終結果χ2=8.727，P=0.273>0.05，表明簡約模型與飽和模型間無差異（H0假設是簡約模型與飽和模型無差異。如果結果顯著則說明現有簡約模型仍然需要加入新的變數以提升模型的解釋力度；相反若結果不顯著說明簡約模型中包括的自變數已足夠，即解釋力度已與飽和模型無差異）。

（5）模型診斷

可採用殘差分析、迭代記錄、分類圖等。

通過殘差分析可以查找異常值，可通過保存（Save…）中的殘差（ Residuals）部分來進行，一般來說殘差絕對值大於2，該條記錄可能是異常點。殘差圖可參見正態分布與方差齊性的檢驗方法與SPSS操作部分。異常點的查找也可通過選項（Option…）中選中Casewise listing of residuals來進行，本例結果如下：第5和40條記錄的的標準化殘差均大於2，應根據實際情況分析是否為異常值。如確定異常值，需要刪除記錄後重新進行回歸分析。

代記錄在選項（Option…）中選中Iteration history來獲得，健康迭代過程的迭代記錄似然值和自變數係數從迭代開始就向著一個方向發展，如中間出現波折，尤其是當引入新變數後變化方向改變了，則提示要好好進一步分析。

分類圖可以直觀的觀測模型的預測狀況，可在選項（Option…）中選中Classification plots來獲得

（6）多重共線的檢驗

運用相同的因變數和自變數擬合線性回歸模型進行共線性診斷

分析>>回歸>>線性（Linear…）

選入相應的因變數和自變數

統計（Statistics…）：選中Collinearity diagnostics，Continue，OK

查看結果中的回歸係數表（Coefficients），如容忍度（Tolerance）小於0.1，方差膨脹因子（VIF）大於10，則提示共線性的存在。存在共線性最簡單的處理方法就是剔除引起共線性的因素。

END

轉自個人微信公眾號【Memo_Cleon】的統計學習筆記

http://weixin.qq.com/r/wjv27irEQHc-rW6N927Q (二維碼自動識別)