非配對設計多值名義資料多水平多重Logistic回歸模型

李長平，張甜甜，宋德勝，胡良平

（1.天津醫科大學公共衛生學院衛生統計學教研室，天津300070；2.世界中醫藥學會聯合會臨床科研統計學專業委員會，北京100029；3.軍事科學院研究生院，北京100850*通信作者：李長平，E-mail：1067181059@qq.com）

1 基本概念

當反應變量為無序離散型結局變量且具有3個或更多類別時，這樣的結局變量稱為多值名義結果變量。在眾多因素中，欲探討哪些因素對多值名義結果變量具有統計學意義的影響，可采用多值名義資料多重logistic回歸模型進行分析。該模型是二值資料多重logistic回歸模型的擴展。如本刊上一期科研方法專題已發表的文章［1］所述，多水平數據具有非獨立性，需要采用能處理多水平數據的多水平模型。

對于一個具有M個類別的名義結果變量的多水平資料，可采用廣義線性混合效應模型進行分析。該模型是二值結果變量多水平logistic回歸模型的擴展。該模型將構建（M-1）個logistic回歸模型，且估計（M-1）組參數［2］。

2 問題與數據結構

【例1】以美國國家毒品濫用研究所開展的一項以社區為基礎的艾滋病干預研究數據［1］為例。共收集20個調查點，包含9 824名靜脈注射吸毒者的基本情況。收集的信息包括受檢者年齡、性別、種族、受教育程度以及所在地區的HIV流行程度和調查前30天內毒品的注射次數等，試研究各因素對靜脈注射吸毒程度的影響。部分結果見表1。

表1 9 824名靜脈注射吸毒者基礎情況

3 非配對設計多值名義資料多水平多重logistic回歸模型的構建原理

3.1 多值名義資料logistic回歸模型

當名義結局變量y具有M個類別時，即m=1，2，…，M，可用式（1）表示多值名義資料logistic回歸模型：

假設結局變量有M個類別，則會擬合（M-1）個logistic回歸模型，會產生（M-1）個截距和（M-1）組斜率估計值。在這（M-1）個logistic回歸模型中，每個對數發生比都是多個結局中的一個類別與參考類別進行比較。表達式如下：

在此前的模型中，假定以第M個類別為參照類別。由于Pr(y=1)+Pr(y=2)+…+Pr(y=M)=1，則：

3.2 多值名義資料多水平logistic回歸模型

當結局變量為多值名義變量且具有M個類別時，反應變量取第m類值的多值名義多水平logistic回歸模型可用廣義線性混合效應模型表達［2-3］，其表達式為：

在式（5）中，m=1，2，…，M；設計矩陣X的固定效應向量為βm；設計矩陣Z的隨機效應向量為Um。該模型假定隨機效應服從正態分布，其均數為零，方差/協方差為矩陣G［即μ～N(0，G)］。當結局變量為具有3個類別的名義變量時，其多值名義資料多水平logistic回歸模型有兩個logistic回歸模型：

其中，我們將結局變量取類別3作為參照類別，且同時擬合兩個logistic回歸模型，估計兩組固定效應（β1和β2）和兩組隨機效應（U1和U2）。根據式（6）和式（7），結局變量為各類別的條件概率模型為：

4 基于SAS分析實例

4.1 分析與解答

例1中，若將響應變量靜脈注射吸毒程度作為多值名義變量，試研究各因素對靜脈注射吸毒程度的影響。該數據中個體可以看作1水平單位，調查點（site）看作2水平單位。可進行多值名義資料多水平多重logistic回歸模型的構建。基于此，結局變量吸毒程度作為具有三個類別的名義變量，且隨機系數為截距和水平1解釋變量race的斜率，而其他水平1協變量，如sex、cenage、educ等都為有固定效應的自變量。另外，水平2協變量只有一個，即region。其模型表達式為：

這里，同時對兩個logistic回歸模型進行分析。SAS程序如下：

【說明】程序共分為3步，第1步是導入數據集，后面是分析過程，分別使用的是GLIMMIX過程和NLMIXED過程。

在GLIMMIX過程中，method選項指定廣義線性混合模型的參數估計方法是虛擬殘差似然法。Class語句指定了site和inject為分類變量。Model語句等號前指定inject為因變量，等號后指定模型中的自變量。S指定顯示固定效應的解，dist指定響應變量的分布類型為多項式分布，link指定連接函數。由于結局變量為多值名義變量，故此處指定連接函數為廣義logit函數，ddfm選項指定分配自由度的方式，參數值bw表示如果固定效應在一個對象內發生變化，分配對象內自由度，否則分配對象間自由度。Random語句指定了隨機效應為截距和race。Subject選項指定了模型中的對象，group指定協方差參數的分組。Nloptions語句中的tech指定了使用Newton-Raphson嶺穩定優化法進行非線性參數估計的優化，以解決在某些分布中默認的二元準牛頓算法存在的收斂問題。

NLMIXED過程中parms語句設定了一些系數的初始值，這些初始值來自于GLIMMIX的參數估計值。Model語句中“～”之前是因變量，后面是因變量服從的分布以及相關參數。General（ll）指定model語句之前使用編程語句構造的廣義對數似然函數。Random語句指定隨機效應以及它的分布。Subject選項定義隨機效應的唯一識別變量。

4.2 主要輸出結果及解釋

以上SAS結果是GLIMMIX的輸出結果。SAS輸出的Dimensions顯示模型有6個隨機效應，而Convariance Parameter Estimates僅顯示了4個隨機效應，其原因是隨機效應不適用于參照組。由于GLIMMIX過程不提供隨機截距的統計顯著性檢驗，因此我們將采用NLMIXED過程確認其統計學顯著性。

SAS輸出的Solutions for Fixed Effects部分列出14個固定參數估計值，即每個logistic回歸模型有7個參數估計值。例如，Intercept1和Intercept2分別是第1個和第2個logit的截距系數估計值。控制協變量之后，類別1發生幾率小于類別3發生幾率：優勢比為exp(-0.9002)=0.41(P=0.0022)；發生類別2與類別3的幾率沒有統計學差異：優勢比為exp(-0.0773)=0.93(P=0.6700)。兩 個logit中，region均有統計學差異，類別1與類別3比較，發生比率為exp(1.4047)=4.07(P=0.0068)；類別2與類別3比較，發生比率為exp(0.9622)=2.67(P=0.0115)。表明高HIV流行區的靜脈注射吸毒者更易成為重、中度吸毒者。值得注意的是，多水平logistic回歸模型采取了多次logit變換，因此解釋協變量效應時具體的logit應與效應一一對應。region和race之間的跨層交互作用在兩個logit中均無統計學意義（P值分別為0.3696和0.7121）。

以上是NLMIXED過程的輸出結果。在NLMIXED過程的PARMS語句中，帶有前綴GA和GB的參數分別是第1個和第2個logistic回歸模型的系數，該回歸系數的設定來源于GLIMMIX過程的結果。Eta1和eta2分別是兩個logit函數的線性預測指標。根據公式（7）（8）（9），IF THEN語句設定結局變量的類比概率。

Random語句設定4個隨機效應的方差/協方差矩陣，即logit1的u0a、u1a和logit2的u0b、u1b。隨機效應的均數設定為零，方差/協方差是根據G矩陣的下三角定義的，其對角線上的元素為方差，對角線以外的元素為協方差。結果顯示，logit1和logit2的隨機截距方差均有統計學意義（v_u0a=0.8，P=0.0210；V_U0B=0.37，P=0.0143）；但變量race的隨機斜率方差無統計學意義（v_u1a=0.15，P=0.0827；V_U1B=0.06，P=0.2333）。結果表明，種族對成為重度注射吸毒類別的幾率的效應不隨各項目實施點顯著變化。以上兩個過程的輸出結果顯示，NLMIXED和GLIMMIX兩個過程對固定效應的估算結果相近。

【結論】對結局變量而言，HIV流行程度、種族、年齡和性別的影響都是不可忽視的。至于每個影響因素各水平對結果影響的差異情況，需結合每個因素的參照類別和回歸系數的正負號，方可給出具體的解釋。

5 討論與小結

對于響應變量為多值名義變量的資料，一般采用廣義logistic回歸模型。但是該模型要求觀測結局相互獨立。當資料為多層的多值名義資料時，考慮到數據之間的聚集性，應當采用多水平回歸模型進行分析，這樣可以使個體的隨機誤差更純［3-4］。

我們采用GLIMMIX和NLMIXED兩個過程來構建模型：前者構建模型的速度快且用法簡單，但在模型比較時通常不適用，且沒有提供隨機效應的假設檢驗，不能采用t檢驗計算相應的P值；而NLMIXED過程可以提供真實的對數似然值，并提供隨機效應假設檢驗的結果，也可以通過似然比檢驗對嵌套模型的擬合效果進行比較，但用法復雜，需設置模型參數的初始值。因此，一般以GLIMMIX過程得到的參數估計值作為NLMIXED過程的模型參數初始值，最后以NLMIXED過程的結果為準。