非配對設計多值有序資料多水平多重Logistic回歸分析

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

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

1 基本概念

1.1 多值有序資料

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

1.2 多水平的概念

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

1.3 多重logistic回歸模型

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

2 數據結構

【例1】研究者選擇8所醫(yī)院開展多中心臨床試驗，每所醫(yī)院均選取400名受試者，在各醫(yī)院內隨機等分成兩組，分別接受試驗藥物和對照藥物治療，治療結果為多值有序變量（好、一般、差），試比較兩種藥物的療效。基本信息見表1。

表1 多中心臨床試驗的基本信息

3 回歸模型的構建與求解

3.1 模型的構建

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

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

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

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

多水平累積logistic回歸模型是對固定效應和隨機效應做了更細致的考察，其模型可以表達如下：

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

3.2 模型的參數估計和假設檢驗

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

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

4 SAS程序及結果解釋

4.1 SAS程序

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

4.2 主要輸出結果及解釋

以下為GLIMMIX過程方差/協方差參數估計的結果，給出了隨機效應方差的估計值。其中隨機截距的方差（即）的估計值為0.4447，標準誤為0.1243。但此處未給出隨機截距方差是否為0的假設檢驗結果，故不能判斷與0之間的差異是否有統計學意義，尚不能說明是否存在隨機效應。

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

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

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

5 討論與小結

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

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

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