鮑曉蕾,王小利,胡良平
對(duì)于結(jié)果變量為多值有序變量的高維列聯(lián)表資料可使用cmH 校正的秩和檢驗(yàn)或有序變量多重 logistic 回歸分析等統(tǒng)計(jì)分析方法。在上一期中,我們已經(jīng)詳細(xì)介紹了CMH 校正的秩和檢驗(yàn)。本期將詳細(xì)介紹多值有序變量的logistic 回歸分析。
結(jié)果變量為多值有序變量的 logistic 回歸又稱為累計(jì)logistic 回歸。累積 logistic 回歸模型可視為二值變量logistic 回歸的擴(kuò)展,其回歸模型可定義如下:
其中y* 表示觀測(cè)現(xiàn)象的內(nèi)在趨勢(shì),不能被直接測(cè)量;ε為誤差項(xiàng)。當(dāng)結(jié)果變量有J個(gè)可能的結(jié)局,相應(yīng)的取值為y= 1、y= 2 ……y=J時(shí),共有J– 1 個(gè)分界點(diǎn)將各相鄰類別分開。即:
若y* ≤μ1,則y= 1;
若μ1 …… 若y* >μJ–1,則y=J。 給定x值的累積概率可以按如下形式表示: 與二值變量的 logit變換相似,累積 logit 變換定義如下: 其中 1 –P(y≤j|x) 即為P(y≥j+ 1 |x),這樣就依次將J個(gè)可能的結(jié)局合并成兩個(gè),從而進(jìn)行 logistic 回歸分析。 相應(yīng)地,累積概率可通過以下公式進(jìn)行預(yù)測(cè): 統(tǒng)計(jì)軟件在實(shí)際運(yùn)行中,定義β0j為各類中截距α與分界點(diǎn)μj的綜合,所以上式就轉(zhuǎn)化為: 值得注意的是,SPSS 和 SAS 在對(duì)累積 logistic 回歸模型進(jìn)行參數(shù)化時(shí)采用的形式是不同的,SPSS 軟件中采用的線性形式是,這與式 ⑸ 相同,而 SAS 中采用的是,所以式 ⑸ 就轉(zhuǎn)化為: 在使用不同軟件時(shí),應(yīng)該注意回歸系數(shù)符號(hào)的差別。 由上面的討論可以看出,若結(jié)果變量有J個(gè)可能的結(jié)局,則可獲得J– 1 個(gè)累積 logit 函數(shù)(當(dāng)進(jìn)行統(tǒng)計(jì)分析時(shí),若有m個(gè)截距項(xiàng)β0j無統(tǒng)計(jì)學(xué)意義,則只能獲得J–m– 1 個(gè)累積 logit 函數(shù))。累積 logistic 回歸模型對(duì)每一個(gè)累積 logit 函數(shù)各有一個(gè)不同的β0j估計(jì),然而對(duì)所有的累積 logit 函數(shù),變量xk卻有一個(gè)相同的βk估計(jì),因?yàn)槠浼僭O(shè)條件為自變量的作用與所有累積 logit 的截?cái)帱c(diǎn)無關(guān)。在此假設(shè)條件下,不同累積對(duì)數(shù)發(fā)生比的回歸線相互平行,只是截距參數(shù)有所不同。這被稱為成比例發(fā)生比假設(shè)條件或平行線假設(shè)條件。 運(yùn)用累積 logistic 回歸模型首先需要對(duì)平行線假設(shè)條件進(jìn)行檢驗(yàn)。如果這一假設(shè)條件被拒絕,便說明自變量xk對(duì)不同的 logit 有不同的βk,因而說明累積 logistic 回歸模型不適合,需要采用其他模型來進(jìn)行資料的分析,如可在模型中引入二次項(xiàng)或交互項(xiàng)。 在 logistic 回歸模型中,原因變量可以是連續(xù)型變量,也可以是二值變量、多值有序變量或多值名義變量。 如果原因變量是連續(xù)變量,一般不需要進(jìn)行處理,可直接建立 logistic 回歸模型。但有時(shí)根據(jù)專業(yè)知識(shí)需對(duì)其進(jìn)行分級(jí)以獲得更有實(shí)際意義的結(jié)果時(shí),連續(xù)變量就轉(zhuǎn)換成了有序變量,此時(shí)可按影響結(jié)果變量由小到大的順序賦值為1、2 ……,并將它當(dāng)作連續(xù)型變量處理或直接引入啞變量,建立 logistic 回歸模型。比如在肺癌危險(xiǎn)因素的病例-對(duì)照研究中,研究者往往感興趣的是年齡每增加 5 歲(根據(jù)專業(yè)知識(shí)和試驗(yàn)?zāi)康臎Q定)肺癌發(fā)病的危險(xiǎn)性是基礎(chǔ)狀態(tài)時(shí)的多少倍,而年齡每增加 1 歲肺癌發(fā)病的危險(xiǎn)性是基礎(chǔ)年齡時(shí)的多少倍往往沒有多大實(shí)際意義。 如果原因變量是二值變量,一般可按 0、1 賦值。 如果原因變量是多值名義變量,需引入啞變量(dummy variable),每個(gè)啞變量都是一個(gè)二值變量,所需啞變量的數(shù)目為多值名義變量的類別數(shù)減 1。如“血型”是一個(gè)多值名義變量,有 A、B、AB、O 四種,若以 O 型血為基準(zhǔn),需引入 3 個(gè)(4 – 1 = 3)啞變量來描述。令 則可得到下面的對(duì)應(yīng)關(guān)系(表1): 表1 用 3 個(gè)啞變量描述一個(gè)四值名義的血型變量 如果原因變量是多值有序變量,一般可按影響結(jié)果變量由小到大的順序賦值為1、2 ……,并將它當(dāng)作連續(xù)型變量處理,但這樣做有時(shí)并不科學(xué),因?yàn)樵撚行蜃兞棵可仙蛳陆狄粋€(gè)等級(jí)對(duì)結(jié)果變量的影響并非總是“線性”的。若遇到有序原因變量按連續(xù)變量處理所得回歸效果不理想時(shí),建議引入啞變量的方式來處理。 結(jié)果變量為多值有序變量,可按其程度或一般默認(rèn)的順序由小到大地賦值為1、2 ……。 【例 1】在一項(xiàng)臨床試驗(yàn)中,研究病程與依沙酰胺療效的關(guān)系。試驗(yàn)在三所醫(yī)院中同時(shí)進(jìn)行,具體療效數(shù)據(jù)見表2,試對(duì)病程與治療效果之間的關(guān)系進(jìn)行分析。 表2 病程與依沙酰胺療效的數(shù)據(jù) SAS 程序如下,程序名為example3。 data example3;do a=1 to 3;do b=1 to 4;do c=1 to 3;input f@@;output;end; end; end;cards;24 10 3 10 4 2 35 26 10 10 8 4 25 12 3 11 4 2 30 22 8 10 10 5 20 8 2 10 4 1 20 17 5 12 10 6;run;ods html;proc logistic data=example3;class a b/param=ref;model c=a b/selection=stepwise;freq f;run;ods html close; 程序說明:首先建立數(shù)據(jù)集,程序中的 a 表示試驗(yàn)中心,a = 1 表示中心 1,a = 2 表示中心 2,a = 3 表示中心 3;b 表示病程,b = 1 表示 < 1 個(gè)月組,b = 2 表示 1~3 個(gè)月組,b = 3 表示 3 個(gè)月~5年組,b = 4表示 > 5年組;c 表示療效,c = 1 表示治愈,c = 2 表示好轉(zhuǎn),c = 3 表示無效;變量 f 表示頻數(shù)。調(diào)用 logistic 過程進(jìn)行有序變量多重 logistic 回歸分析。由于試驗(yàn)中心是多值名義變量,需對(duì)其賦啞變量;病程是多值有序變量,可將其當(dāng)成連續(xù)型變量直接賦值,但最好是對(duì)其賦啞變量。若自變量是二值變量,則只需直接賦值 0、1 即可。class 語句可實(shí)現(xiàn)對(duì)自變量自動(dòng)賦啞變量,同時(shí)還能保證啞變量在回歸方程中同進(jìn)同出。class a b 表示對(duì) a、b 兩個(gè)因素均自動(dòng)賦啞變量。選項(xiàng)param = ref 指定將其中的一個(gè)水平作為基準(zhǔn)實(shí)現(xiàn)啞變量賦值,默認(rèn)以輸入的該定性變量最后一個(gè)水平為基準(zhǔn)。若想以其他水平為基準(zhǔn),比如希望以 a 的第一個(gè)水平為基準(zhǔn)對(duì)因素 a 賦啞變量,則只需將語句改成 class a (ref = ‘1’)b/param = ref 即可。model 語句表示建模,等號(hào)前表示因變量,等號(hào)后表示自變量。選項(xiàng) selection = stepwise 表示用逐步回歸法進(jìn)行變量篩選,其默認(rèn)的進(jìn)入和剔除方程的顯著性水準(zhǔn)為0.05;若希望改變這一標(biāo)準(zhǔn),可在選項(xiàng)后加入 sle =xx 和 sls = xx 選項(xiàng),前者指定進(jìn)入方程的標(biāo)準(zhǔn),后者指定剔除方程的標(biāo)準(zhǔn),xx 表示具體的數(shù)值(0~1 之間)。freq f指定 f 變量為頻數(shù)變量。 SAS 程序運(yùn)行結(jié)果: 以上是通過 class 語句對(duì)因素 a 和因素 b 自動(dòng)賦啞變量的結(jié)果??梢钥吹?,兩因素均以最后一個(gè)水平為基準(zhǔn)。以 a 因素為例,啞變量 a1 表示 a 因素的第一個(gè)水平相對(duì)于第三個(gè)水平進(jìn)行分析;啞變量 a2 表示 a 因素的第二個(gè)水平相對(duì)于第三個(gè)水平進(jìn)行分析。b 因素的啞變量意義類似。 Summary of stepwise selection 以上為逐步篩選法的篩選結(jié)果,最終只有因素 b 進(jìn)入了回歸方程(χ2= 13.4863,P= 0.0037)。 Score test for the proportional odds assumption 以上為平行線假設(shè)的檢驗(yàn)結(jié)果:χ2= 0.8909,P= 0.8276> 0.05,說明資料滿足平行線假設(shè)。 Type 3 analysis of effects 以上是將因素 b 作為一個(gè)整體的假設(shè)檢驗(yàn)結(jié)果:waldχ2= 13.5102,P= 0.0037 < 0.05,說明病程對(duì)療效的影響有統(tǒng)計(jì)學(xué)意義。 Analysis of maximum likelihood estimates 以上為參數(shù)估計(jì)及假設(shè)檢驗(yàn)的結(jié)果。本例結(jié)果變量有 3個(gè)水平,故模型包含 2 個(gè)截距項(xiàng)。若P1、P2和P3分別表示治愈、好轉(zhuǎn)和無效發(fā)生的概率,則回歸方程為: Odds ratio estimates 以上是對(duì)優(yōu)勢(shì)比的估計(jì)結(jié)果:OR1vs4= 2.576,其 95%置信區(qū)間為(1.445,4.592);OR2vs4= 2.473,其 95%置信區(qū)間為(1.205,5.072);OR3vs4= 1.382,其 95%置信區(qū)間為(0.830,2.302)。 統(tǒng)計(jì)及專業(yè)結(jié)論:病程對(duì)療效的影響有統(tǒng)計(jì)學(xué)意義(P= 0.0037 < 0.05),而試驗(yàn)中心對(duì)療效的影響沒有統(tǒng)計(jì)學(xué)意義。OR1vs4= 2.576,其 95%置信區(qū)間為(1.445,4.592),說明病程小于 1 個(gè)月的療效是病程大于 5年療效的2.576倍;OR2vs4= 2.473,其 95%置信區(qū)間為(1.205,5.072),說明病程為1~3 個(gè)月的療效是病程大于 5年療效的 2.473 倍;OR3vs4= 1.382,其 95%置信區(qū)間為(0.830,2.302),說明病程為3 個(gè)月~5年的療效是病程大于 5年療效的 1.382 倍。 [1]Hu LP.Statistics facing practical scientific issues -- (2) multi-factor designs and linear model analysis.Beijing: People’s Medical Publishing House, 2012:508-517.(in Chinese)胡良平.面向問題的統(tǒng)計(jì)學(xué)——(2)多因素設(shè)計(jì)與線性模型分析.北京: 人民衛(wèi)生出版社, 2012:508-517. [2]Hu LP.Medical statistics-analysis of quantitative and qualitative data applying the triple-type theory.Beijing: People’s Military Medical Press, 2009:363-375.(in Chinese)胡良平.醫(yī)學(xué)統(tǒng)計(jì)學(xué)-運(yùn)用三型理論分析定量與定性資料.北京:人民軍醫(yī)出版社, 2009:363-375.2 變量賦值