如何用SAS軟件正確分析生物醫(yī)學(xué)科研資料XXIV. 結(jié)果變量為多值有序變量的高維列聯(lián)表資料的統(tǒng)計(jì)分析與SAS軟件實(shí)現(xiàn)（二）

鮑曉蕾，王小利，胡良平

·講座·

如何用SAS軟件正確分析生物醫(yī)學(xué)科研資料XXIV. 結(jié)果變量為多值有序變量的高維列聯(lián)表資料的統(tǒng)計(jì)分析與SAS軟件實(shí)現(xiàn)（二）

鮑曉蕾，王小利，胡良平

編者按

生物統(tǒng)計(jì)學(xué)是生物學(xué)領(lǐng)域科學(xué)研究和實(shí)際工作中必不可少的工具，在分子生物學(xué)迅速發(fā)展的今天，生物統(tǒng)計(jì)學(xué)更顯示出了它的重要性。實(shí)驗(yàn)設(shè)計(jì)與數(shù)據(jù)統(tǒng)計(jì)分析是現(xiàn)代生物學(xué)的基石，是生物學(xué)研究者檢驗(yàn)假說(shuō)、尋找模式、建立生物學(xué)理論的有利工具，也是生物學(xué)研究者探索微觀和宏觀生物世界的必備基礎(chǔ)知識(shí)。對(duì)于每天甚至是每時(shí)每刻涌現(xiàn)的大量的、以天文數(shù)字計(jì)量的分子遺傳數(shù)據(jù)，必須借助統(tǒng)計(jì)學(xué)知識(shí)加以分析處理，才能從中獲得有意義的信息?！吧锒鄻有詳?shù)據(jù)分析”是開(kāi)展生物多樣性研究的一個(gè)重要方面，數(shù)據(jù)分析能力的高低極大地影響著我們對(duì)各種生態(tài)學(xué)現(xiàn)象認(rèn)識(shí)的深度和廣度。現(xiàn)在，電子計(jì)算機(jī)的普及使得生物統(tǒng)計(jì)分析過(guò)程大大簡(jiǎn)化，生物統(tǒng)計(jì)分析軟件包的普及將生物統(tǒng)計(jì)學(xué)從統(tǒng)計(jì)學(xué)家的書(shū)本里解放了出來(lái)，簡(jiǎn)化了生物統(tǒng)計(jì)分析過(guò)程，使之成為生物學(xué)研究者的常用工具。本刊特邀軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計(jì)學(xué)咨詢(xún)中心主任胡良平教授，以“如何用 SAS 軟件正確分析生物醫(yī)學(xué)科研資料”為題，撰寫(xiě)系列統(tǒng)計(jì)學(xué)講座，希望該系列講座能對(duì)生物醫(yī)學(xué)科研工作者有所幫助。

對(duì)于結(jié)果變量為多值有序變量的高維列聯(lián)表資料可使用 CMH 校正的秩和檢驗(yàn)或有序變量多重 logistic 回歸分析等統(tǒng)計(jì)分析方法。在上一期中，我們已經(jīng)詳細(xì)介紹了CMH 校正的秩和檢驗(yàn)。本期將詳細(xì)介紹多值有序變量的logistic 回歸分析。

1 原理

結(jié)果變量為多值有序變量的 logistic 回歸又稱(chēng)為累計(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)將各相鄰類(lèi)別分開(kāi)。即：

若 y* ≤ μ1，則 y = 1；

若 μ1< y* ≤ μ2，則 y = 2；

……

若 y* > μJ–1，則 y = J。

給定 χ 值的累積概率可以按如下形式表示：

與二值變量的 logit變換相似，累積 logit 變換定義如下：

其中 1 – P(y ≤ j | χ) 即為 P(y ≥ j + 1 | χ)，這樣就依次將 J 個(gè)可能的結(jié)局合并成兩個(gè)，從而進(jìn)行 logistic 回歸分析。

相應(yīng)地，累積概率可通過(guò)以下公式進(jìn)行預(yù)測(cè)：

統(tǒng)計(jì)軟件在實(shí)際運(yùn)行中，定義 β0j為各類(lèi)中截距 α 與分界點(diǎn) μj的綜合，所以上式就轉(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無(wú)統(tǒng)計(jì)學(xué)意義，則只能獲得J – m – 1 個(gè)累積 logit 函數(shù)）。累積 logistic 回歸模型對(duì)每一個(gè)累積 logit 函數(shù)各有一個(gè)不同的 β0j估計(jì)，然而對(duì)所有的累積 logit 函數(shù)，變量 χk卻有一個(gè)相同的 βk估計(jì)，因?yàn)槠浼僭O(shè)條件為自變量的作用與所有累積 logit 的截?cái)帱c(diǎn)無(wú)關(guān)。在此假設(shè)條件下，不同累積對(duì)數(shù)發(fā)生比的回歸線(xiàn)相互平行，只是截距參數(shù)有所不同。這被稱(chēng)為成比例發(fā)生比假設(shè)條件或平行線(xiàn)假設(shè)條件。

運(yùn)用累積 logistic 回歸模型首先需要對(duì)平行線(xiàn)假設(shè)條件進(jìn)行檢驗(yàn)。如果這一假設(shè)條件被拒絕，便說(shuō)明自變量 χk對(duì)不同的 logit 有不同的 βk，因而說(shuō)明累積 logistic 回歸模型不適合，需要采用其他模型來(lái)進(jìn)行資料的分析，如可在模型中引入二次項(xiàng)或交互項(xiàng)。

2 變量賦值

在 logistic 回歸模型中，原因變量可以是連續(xù)型變量，也可以是二值變量、多值有序變量或多值名義變量。

如果原因變量是連續(xù)變量，一般不需要進(jìn)行處理，可直接建立 logistic 回歸模型。但有時(shí)根據(jù)專(zhuān)業(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ù)專(zhuān)業(yè)知識(shí)和試驗(yàn)?zāi)康臎Q定）肺癌發(fā)病的危險(xiǎn)性是基礎(chǔ)狀態(tài)時(shí)的多少倍，而年齡每增加 1 歲肺癌發(fā)病的危險(xiǎn)性是基礎(chǔ)年齡時(shí)的多少倍往往沒(méi)有多大實(shí)際意義。

如果原因變量是二值變量，一般可按 0、1 賦值。

如果原因變量是多值名義變量，需引入啞變量（dummy variable），每個(gè)啞變量都是一個(gè)二值變量，所需啞變量的數(shù)目為多值名義變量的類(lèi)別數(shù)減 1。如“血型”是一個(gè)多值名義變量，有 A、B、AB、O 四種，若以 O 型血為基準(zhǔn)，需引入 3 個(gè)（4 – 1 = 3）啞變量來(lái)描述。令

SAS 程序如下，程序名為 example3。

則可得到下面的對(duì)應(yīng)關(guān)系（表 1）：

表1 用 3 個(gè)啞變量描述一個(gè)四值名義的血型變量

如果原因變量是多值有序變量，一般可按影響結(jié)果變量由小到大的順序賦值為 1、2 ……，并將它當(dāng)作連續(xù)型變量處理，但這樣做有時(shí)并不科學(xué)，因?yàn)樵撚行蜃兞棵可仙蛳陆狄粋€(gè)等級(jí)對(duì)結(jié)果變量的影響并非總是“線(xiàn)性”的。若遇到有序原因變量按連續(xù)變量處理所得回歸效果不理想時(shí)，建議引入啞變量的方式來(lái)處理。

結(jié)果變量為多值有序變量，可按其程度或一般默認(rèn)的順序由小到大地賦值為 1、2 ……。

【例 1】 在一項(xiàng)臨床試驗(yàn)中，研究病程與依沙酰胺療效的關(guān)系。試驗(yàn)在三所醫(yī)院中同時(shí)進(jìn)行，具體療效數(shù)據(jù)見(jiàn)表 2，試對(duì)病程與治療效果之間的關(guān)系進(jìn)行分析。

表2 病程與依沙酰胺療效的數(shù)據(jù)

程序說(shuō)明：首先建立數(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 表示無(wú)效；變量 f 表示頻數(shù)。調(diào)用 logistic 過(guò)程進(jìn)行有序變量多重 logistic 回歸分析。由于試驗(yàn)中心是多值名義變量，需對(duì)其賦啞變量；病程是多值有序變量，可將其當(dāng)成連續(xù)型變量直接賦值，但最好是對(duì)其賦啞變量。若自變量是二值變量，則只需直接賦值 0、1 即可。class 語(yǔ)句可實(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 賦啞變量，則只需將語(yǔ)句改成 class a (ref = ‘1’) b/param = ref 即可。model 語(yǔ)句表示建模，等號(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 level information

以上是通過(guò) class 語(yǔ)句對(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 因素的啞變量意義類(lèi)似。

Summary of stepwise selection

以上為逐步篩選法的篩選結(jié)果，最終只有因素 b 進(jìn)入了回歸方程（χ2= 13.4863，P = 0.0037）。

Score test for the proportional odds assumption

以上為平行線(xiàn)假設(shè)的檢驗(yàn)結(jié)果：χ2= 0.8909，P = 0.8276 > 0.05，說(shuō)明資料滿(mǎn)足平行線(xiàn)假設(shè)。

Type 3 analysis of effects

以上是將因素 b 作為一個(gè)整體的假設(shè)檢驗(yàn)結(jié)果：wald χ2= 13.5102，P = 0.0037 < 0.05，說(shuō)明病程對(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)和無(wú)效發(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ì)及專(zhuān)業(yè)結(jié)論：病程對(duì)療效的影響有統(tǒng)計(jì)學(xué)意義（P = 0.0037 < 0.05），而試驗(yàn)中心對(duì)療效的影響沒(méi)有統(tǒng)計(jì)學(xué)意義。OR1vs4= 2.576，其 95% 置信區(qū)間為（1.445，4.592），說(shuō)明病程小于 1 個(gè)月的療效是病程大于 5 年療效的2.576倍；OR2vs4= 2.473，其 95% 置信區(qū)間為（1.205，5.072），說(shuō)明病程為 1 ~ 3 個(gè)月的療效是病程大于 5 年療效的 2.473 倍；OR3vs4= 1.382，其 95% 置信區(qū)間為（0.830，2.302），說(shuō)明病程為 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)

胡良平. 面向問(wèn)題的統(tǒng)計(jì)學(xué)——(2)多因素設(shè)計(jì)與線(xiàn)性模型分析.北京: 人民衛(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.

10.3969/cmba.j.issn.1673-713X.2013.04.014

100850 北京，軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計(jì)學(xué)咨詢(xún)中心（鮑曉蕾、胡良平）；100850 北京，軍事醫(yī)學(xué)科學(xué)院門(mén)診部（王小利）

胡良平，Email：lphu812@sina.com