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

李長平，張?zhí)鹛?，宋德勝，胡良?/p>

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

1 基本概念

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

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

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

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

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

3 非配對設(shè)計多值名義資料多水平多重logistic回歸模型的構(gòu)建原理

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

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

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

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

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

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

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

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

4 基于SAS分析實例

4.1 分析與解答

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

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

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

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

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

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

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

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

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

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

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

5 討論與小結(jié)

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

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