非配對設(shè)計(jì)多值有序資料多水平多重Logistic回歸分析

鳳思苑，李長平，2，胡良平

（1.天津醫(yī)科大學(xué)公共衛(wèi)生學(xué)院衛(wèi)生統(tǒng)計(jì)學(xué)教研室，天津300070；2.世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計(jì)學(xué)專業(yè)委員會，北京100029；3.軍事科學(xué)院研究生院，北京100850*通信作者：胡良平，E-mail：lphu812@sina.com）

1 基本概念

1.1 多值有序資料

“多值有序”資料特指因變量或結(jié)局變量為多值有序變量（例如在描述藥物或手術(shù)療效時(shí)經(jīng)常用“治愈、顯效、好轉(zhuǎn)、無效和死亡”作為一個主要療效指標(biāo)的不同取值），而自變量沒有任何的限制，可以是定量的或定性的（包括二值的、多值有序的、多值名義的）變量。

1.2 多水平的概念

在社會科學(xué)研究中，“社會”的基本概念是一個具有分級結(jié)構(gòu)的整體。所謂的分級結(jié)構(gòu)就是指較低層次的單位嵌套在較高層次的單位之下，而這種社會分級結(jié)構(gòu)自然而然的使其產(chǎn)生的數(shù)據(jù)呈現(xiàn)多層次（多水平）結(jié)構(gòu)［1］。例如，在對學(xué)生成績的研究中，認(rèn)為學(xué)生的學(xué)習(xí)成績或狀態(tài)不僅與個人的內(nèi)在因素有關(guān)，還與所處的環(huán)境（學(xué)校、班級）有關(guān)，因此，在研究學(xué)生成績與個體水平變量的數(shù)量關(guān)系時(shí)，還需將其嵌套到相應(yīng)的學(xué)校和班級中去。由此形成3個層次的結(jié)構(gòu)數(shù)據(jù)：第一個層次的觀察單位是學(xué)生，第二個層次的觀察單位是班級，第三個層次的觀察單位是學(xué)校。這里的“多水平”是指層次結(jié)構(gòu)數(shù)據(jù)中的多個層次，其中學(xué)生為低水平即水平1單位，班級為中水平即水平2單位，而學(xué)校則為高水平即水平3單位；而在通常的回歸分析中，只有一種觀察單位，那就是“個體”或“受試對象”。此時(shí)，若資料中出現(xiàn)了“學(xué)校”“班級”等變量，則它們就被視為定性的“影響因素”（即自變量），通常需要將它們產(chǎn)生啞變量后引入回歸模型中去［2］。

1.3 多重logistic回歸模型

多重logistic回歸模型是一種廣義線性回歸模型，適用于研究一個定性因變量與多個自變量之間的依賴關(guān)系，其因變量y可以是二值變量、多值名義變量或多值有序變量。它不同于一般的多重線性回歸模型，其本質(zhì)屬于非線性概率回歸模型，在這種回歸模型中，真正的因變量是y取某特定值時(shí)所對應(yīng)的概率［如P(y=0)或P(y=1)］。

2 數(shù)據(jù)結(jié)構(gòu)

【例1】研究者選擇8所醫(yī)院開展多中心臨床試驗(yàn)，每所醫(yī)院均選取400名受試者，在各醫(yī)院內(nèi)隨機(jī)等分成兩組，分別接受試驗(yàn)藥物和對照藥物治療，治療結(jié)果為多值有序變量（好、一般、差），試比較兩種藥物的療效?；拘畔⒁姳?。

表1 多中心臨床試驗(yàn)的基本信息

3 回歸模型的構(gòu)建與求解

3.1 模型的構(gòu)建

分析結(jié)局變量為多值有序變量時(shí)，一般構(gòu)建累積logistic回歸模型，也稱為比例優(yōu)勢模型。累積logistic回歸模型其實(shí)就是結(jié)局變量為二值變量的logistic回歸模型的擴(kuò)展，從潛在變量的概念出發(fā)，模型可定義如下：

其中y*表示觀察現(xiàn)象的內(nèi)在趨勢，不能被直接測量；e為誤差項(xiàng)。當(dāng)實(shí)際的觀測結(jié)果變量有J個不同的類別時(shí)（j=1，2，…，m，…，J），相應(yīng)的取值即為y=1，y=2，…，y=J。于是，（J-1）個分界點(diǎn)將相鄰各類別分開［3-4］。

與結(jié)局變量為二值變量的logit變換類似，logit變換后的累積logistic回歸模型表達(dá)如下：

在該模型中，P（y≤m|x）實(shí)際是結(jié)局變量取值≤m的累積概率，即為P（y=1|x）+P（y=2|x）+…+P（y=m|x）的概率之和。該模型是將結(jié)局變量的J個等級人為分成兩類｛1，2，…，m｝和｛m+1，…，J｝，在這兩類基礎(chǔ)上定義的logit函數(shù)，實(shí)為前m個等級的累積概率與后（J-m）個等級累積概率比值的對數(shù)。該模型中共有（J-1）個累積的logits，β0m是第m個logit的截距，βk是協(xié)變量xk的斜率。模型的一個重要特征就是（J-1）個截距互不相同，但每個logit中相同自變量的系數(shù)相同，故而又稱比例優(yōu)勢模型［1，4］。

多水平累積logistic回歸模型是對固定效應(yīng)和隨機(jī)效應(yīng)做了更細(xì)致的考察，其模型可以表達(dá)如下：

該公式與普通的（單水平）累積logistic回歸模型相似，對應(yīng)了（J-1）個logit，但不同的是：此處的每個logits的截距可能是隨機(jī)系數(shù)，因而可體現(xiàn)宏觀水平（本例為2水平）單位間的差異。公式中的X是含有固定斜率的協(xié)變量的設(shè)計(jì)矩陣，β代表固定效應(yīng)，而Z是含有隨機(jī)斜率的協(xié)變量的設(shè)計(jì)矩陣，U代表隨機(jī)效應(yīng)［1，3-4］。

3.2 模型的參數(shù)估計(jì)和假設(shè)檢驗(yàn)

多水平累積logistic回歸模型由于存在水平1和水平2殘差組成的復(fù)合殘差結(jié)構(gòu)，模型的參數(shù)估計(jì)較為復(fù)雜，需同時(shí)估計(jì)固定回歸系數(shù)、隨機(jī)回歸系數(shù)以及矩陣G和R的方差/協(xié)方差矩陣（矩陣G為水平2殘差的方差/協(xié)方差矩陣、矩陣R為水平1殘差的方差/協(xié)方差矩陣）。目前SAS的GLIMMIX、NLMIXED過程進(jìn)行參數(shù)估計(jì)的方法主要有RSPL、MSPL、RMPL、MMPL，其本質(zhì)都是基于最大似然的估計(jì)方法。

多水平累積logistic回歸模型的假設(shè)檢驗(yàn)包括固定效應(yīng)的假設(shè)檢驗(yàn)、隨機(jī)效應(yīng)的假設(shè)檢驗(yàn)以及模型比較的檢驗(yàn)。固定效應(yīng)即模型中的固定參數(shù)包括總體的截距、協(xié)變量的斜率。隨機(jī)效應(yīng)是指模型中的隨機(jī)部分，主要指宏觀水平（本例為2水平）殘差的方差/協(xié)方差。當(dāng)采用不同的模型擬合相同的數(shù)據(jù)時(shí)，可以用似然比檢驗(yàn)，有關(guān)的統(tǒng)計(jì)量有-2倍的對數(shù)似然值。當(dāng)模型中包含的參數(shù)數(shù)目相同時(shí)，-2倍的對數(shù)似然值越小，模型對數(shù)據(jù)的擬合效果越好。

4 SAS程序及結(jié)果解釋

4.1 SAS程序

【程序說明】程序共3步，包括1個數(shù)據(jù)步和2個過程步。首先建立例1的數(shù)據(jù)集MLMO，利用do循環(huán)語句輸入變量Hospital（醫(yī)院編號）、Drug（藥物類型）、Gender（性別）和結(jié)局變量y（療效類型）。程序第2步調(diào)用GLIMMIX過程運(yùn)行多水平累積logistic回歸模型，其中Class語句創(chuàng)建分類變量Hospital，model語句中設(shè)置y為響應(yīng)變量，“dist=multi”和“l(fā)ink=clogit”選項(xiàng)分別設(shè)定分布為多項(xiàng)式分布，連接函數(shù)為累積logit函數(shù)。Random語句用來設(shè)定隨機(jī)效應(yīng)，“type=chol”選項(xiàng)采用chol-esky分解法來設(shè)定G矩陣，目的是保證G矩陣具有正特征根，以保證模型參數(shù)估計(jì)的穩(wěn)定。程序第三步利用NLMIXED過程實(shí)現(xiàn)多水平累積logistic回歸模型，parms語句給出模型中有關(guān)參數(shù)的初始值，此處初始值為由GLIMMIX過程計(jì)算所得。z為定義的線性預(yù)測值，由固定效應(yīng)部分和隨機(jī)效應(yīng)u組成。

4.2 主要輸出結(jié)果及解釋

以下為GLIMMIX過程方差/協(xié)方差參數(shù)估計(jì)的結(jié)果，給出了隨機(jī)效應(yīng)方差的估計(jì)值。其中隨機(jī)截距的方差（即）的估計(jì)值為0.4447，標(biāo)準(zhǔn)誤為0.1243。但此處未給出隨機(jī)截距方差是否為0的假設(shè)檢驗(yàn)結(jié)果，故不能判斷與0之間的差異是否有統(tǒng)計(jì)學(xué)意義，尚不能說明是否存在隨機(jī)效應(yīng)。

以下為GLIMMIX過程輸出的固定效應(yīng)檢驗(yàn)結(jié)果。模型有兩個截距，這是因?yàn)轫憫?yīng)變量療效有三個水平。在響應(yīng)變量為J個水平的多水平累積logistic回歸模型中，有（J-1）個logits函數(shù)式，這些函數(shù)式中有（J-1）個不同的截距，但會有一組相同的協(xié)變量系數(shù)的估計(jì)值。因模型是以“y=1”為基礎(chǔ)，故截距值-0.4714表示協(xié)變量均取0值時(shí)治療結(jié)果為“好”的對數(shù)發(fā)生比；截距值為0.7312表示協(xié)變量均取0值時(shí)治療結(jié)果為“好”和“一般”的對數(shù)發(fā)生比［注意：療效單獨(dú)為“一般”的截距應(yīng)為“0.7312-（-0.4714）=1.2026”］。正（負(fù)）斜率表示治療效果為“好”的可能性高（低）。例如，Drug的斜率為0.3627（P＜0.0001），表示試驗(yàn)組藥物的治療效果為“好”的概率比對照組藥物治療效果為“好”的概率高［1，5］。此外，還可以在程序中model語句的“/”之后添加選項(xiàng)oddsratio獲得各個協(xié)變量的OR估計(jì)值及95%CI。

NLMIXED過程輸出了與GLIMMIX過程類似的結(jié)果，即模型的總體信息、優(yōu)化信息以及迭代史，其中重要的是模型各參數(shù)的初始值信息：b0為模型的總體截距，b1為性別的效應(yīng)，b2為藥物的效應(yīng)，V_u0為隨機(jī)效應(yīng)的方差，這些參數(shù)的設(shè)定來源于GLIMMIX過程計(jì)算結(jié)果。NLMIXED過程模型的初始參數(shù)如下：

以下為NLMIXED過程輸出的模型的擬合信息和參數(shù)估計(jì)，包括固定效應(yīng)和隨機(jī)效應(yīng)方差的參數(shù)估計(jì)以及相應(yīng)的假設(shè)檢驗(yàn)結(jié)果。其中b0、b+b0、b1和b2分別表示截距1、截距2、Gender和Drug的系數(shù)值。對于隨機(jī)效應(yīng)的假設(shè)檢驗(yàn)，這里進(jìn)行的是雙側(cè)檢驗(yàn)。實(shí)際上，由于方差不可能為負(fù)值，所以檢驗(yàn)殘差方差應(yīng)選用單側(cè)檢驗(yàn)，故此處的V_u0對應(yīng)P值除以2后才是正確的P值，實(shí)際小于0.05，說明確實(shí)存在隨機(jī)效應(yīng)。有關(guān)其他固定效應(yīng)參數(shù)的解釋參考GLIMMIX過程輸出結(jié)果的解釋。當(dāng)然，由于NLMIXED過程所得的結(jié)果提供了隨機(jī)效應(yīng)的假設(shè)檢驗(yàn)，更為精確，最終結(jié)果應(yīng)以NLMIXED過程的輸出結(jié)果為準(zhǔn)。

5 討論與小結(jié)

對非配對的多值有序資料建立logistic回歸模型時(shí)，除了要考慮有充足的樣本量，以保證參數(shù)估計(jì)的穩(wěn)定性，還必須考慮研究個體是否存在聚集性特征。目前醫(yī)學(xué)研究試驗(yàn)設(shè)計(jì)大多數(shù)會產(chǎn)生多層次（即多水平）數(shù)據(jù)，而此類數(shù)據(jù)常存在組內(nèi)相關(guān)的問題，即組內(nèi)觀察值相互間是非獨(dú)立的。這種現(xiàn)象的存在會導(dǎo)致自變量和結(jié)局變量的關(guān)系隨著宏觀水平單位的不同而變化，此時(shí)若依然采用一般的累積logistic回歸模型，會導(dǎo)致錯誤的參數(shù)估計(jì)結(jié)果，而多水平累積logistic回歸模型可以很好地解決組內(nèi)同質(zhì)、組間異質(zhì)數(shù)據(jù)的回歸建模問題。

本文就多水平多值有序數(shù)據(jù)分別利用SAS的GLIMMIX過程和NLMIXED過程來擬合多水平累積logistic回歸模型，結(jié)果發(fā)現(xiàn)兩個過程參數(shù)估計(jì)的結(jié)果極為相似，但仍存在一些區(qū)別：NLMIXED過程的參數(shù)估計(jì)結(jié)果中直接提供了隨機(jī)效應(yīng)的假設(shè)檢驗(yàn)結(jié)果，有利于模型對于隨機(jī)效應(yīng)的取舍，若隨機(jī)效應(yīng)檢驗(yàn)的結(jié)果沒有統(tǒng)計(jì)學(xué)意義，可以直接采用普通的累積logistic回歸模型直接擬合數(shù)據(jù)。GLIMMIX過程并不提供該檢驗(yàn)，但卻為NLMIXED過程的初始參數(shù)設(shè)置提供了參考，極大地縮短了模型擬合的速度。建議二者同時(shí)使用，但以NLMIXED過程的輸出結(jié)果為準(zhǔn)。

采用多水平累積logistic回歸模型分析數(shù)據(jù)時(shí)還需要注意以下問題。①測量中心化：在多水平累積logistic回歸分析中，要注意同時(shí)關(guān)注水平1截距和斜率的變化。因?yàn)榧俣ㄒ粋€水平1截距為1.30的回歸模型，我們可以說當(dāng)模型中所有自變量都為0時(shí)，某種結(jié)局的對數(shù)優(yōu)勢比為1.30。但是所觀察的某些解釋變量若沒有實(shí)際的零值，則上述解釋便無任何實(shí)際意義。此種情況下要使截距變得有意義，必須通過中心化重新定義或轉(zhuǎn)化自變量的測量值［1］。②隨機(jī)效應(yīng)檢驗(yàn)：模型隨機(jī)部分的檢驗(yàn)主要指對宏觀水平殘差的方差/協(xié)方差檢驗(yàn)，根據(jù)定義，方差不能為負(fù)數(shù)，所以檢驗(yàn)殘差方差應(yīng)選用單側(cè)檢驗(yàn)，其統(tǒng)計(jì)量相應(yīng)的P值應(yīng)除以2；其次，用于模型比較的似然比檢驗(yàn)也可以用于隨機(jī)效應(yīng)的檢驗(yàn)。即先將特定的水平1回歸系數(shù)設(shè)定為固定系數(shù)，然后再將其設(shè)定為隨機(jī)系數(shù)，分別擬合并比較以篩選出適宜的模型［1］。