合理進(jìn)行多重Logistic回歸分析
——結(jié)合多水平模型分析

胡純嚴(yán) ，胡良平 ，2*

（1.軍事科學(xué)院研究生院，北京 100850；2.世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計學(xué)專業(yè)委員會，北京 100029*通信作者：胡良平，E-mail：lphu927@163.com）

統(tǒng)計資料可分為無層級結(jié)構(gòu)和有層級結(jié)構(gòu)兩種類型。通常情況下，研究者收集的統(tǒng)計資料為無層級結(jié)構(gòu)的資料。此時，可運用常規(guī)統(tǒng)計分析方法處理資料。然而，在調(diào)查研究中，受試對象常歸屬于不同層級，例如某省某市某縣某區(qū)。又例如，學(xué)生常歸屬于某市某中學(xué)某年級某班級。當(dāng)資料中的受試對象歸屬于多層級結(jié)構(gòu)時，一般來說，不適合選用處理無層級結(jié)構(gòu)資料的統(tǒng)計分析方法。本文將結(jié)合一個具有2個層級結(jié)構(gòu)的常見臨床試驗資料，介紹二水平多重Logistic回歸模型的構(gòu)建方法和SAS實現(xiàn)方法。

1 基本概念

1.1 分層變量

在統(tǒng)計資料中，有一類變量被稱為“分層變量”。例如，在人口普查資料中，會涉及省、市、縣、區(qū)、街道或村委會，它們都可以被稱為“分層變量”。從全國角度考量，可以把“省”作為最高級別的分層變量，把“市”作為第二級別的分層變量，以此類推，把“街道或村委會”作為最低級別的分層變量。

1.2 層級結(jié)構(gòu)

在研究人的身高與體重關(guān)系的資料中，若受試對象的來源涉及省、市、縣3個行政級別，在統(tǒng)計學(xué)上就稱此資料為具有3個層級結(jié)構(gòu)的資料。也就是說，省下面嵌套著市、市下面又嵌套著縣。同理，在一個多中心臨床試驗研究中，常把“試驗中心”和“受試對象”視為2個不同層級上的變量，因為“試驗中心”下面嵌套著“受試對象”。

1.3 k水平單位

在一個涉及省、市、縣3個行政級別的關(guān)于居民健康情況的調(diào)查研究中，受調(diào)查者為1水平單位，縣為2水平單位，市為3水平單位，省為4水平單位；另一種等價的表述為：受調(diào)查者為水平1單位，縣為水平2單位，市為水平3單位，省為水平4單位。

1.4 個體水平與組水平

所謂個體水平，實際上就是指多層級資料中的最底層，通常為“個體”。若試驗中需要對每個個體進(jìn)行多次重復(fù)觀測，則時間點就成為最底層了，此時，時間點就叫做“個體水平”，而每個個體就被稱為“組水平”。在沒有重復(fù)觀測的試驗研究資料中，“個體”就叫做“個體水平”，而它們的上一級就叫做“組水平”。在多層級結(jié)構(gòu)中，必須給“組水平”附上相應(yīng)的標(biāo)記，以示區(qū)別。例如，在一個涉及省、市、縣3個行政級別的關(guān)于居民健康情況的調(diào)查研究中，就有3種“組水平”，即“省級組水平”“市級組水平”和“縣級組水平”。

2 多水平模型的構(gòu)建步驟

多水平模型是基于層級結(jié)構(gòu)數(shù)據(jù)形成的一種統(tǒng)計模型［1］，對于具有層級結(jié)構(gòu)的二值結(jié)果變量的統(tǒng)計資料，需要采用多水平多重Logistic回歸模型分析，它是一般Logistic回歸模型的延伸，通過在模型中納入隨機(jī)效應(yīng)項來處理層級結(jié)構(gòu)數(shù)據(jù)的組內(nèi)相關(guān)［2］。其模型表達(dá)見式（1）。

式（1）中，X是固定效應(yīng)的解釋變量設(shè)計矩陣，Z是隨機(jī)效應(yīng)的解釋變量設(shè)計矩陣，β是水平1固定回歸系數(shù)向量，U是隨機(jī)回歸系數(shù)向量，服從均值為0、協(xié)方差為矩陣G的正態(tài)分布。

下面介紹2水平多重Logistic回歸模型的建模過程及分析步驟。假定在一個臨床試驗中，涉及試驗藥種類（drug）和臨床試驗中心（center）兩個主要的原因變量，還涉及一個代表療效的結(jié)果變量（y）。假設(shè)：y=0表示治療成功，y=1表示治療失敗，以pij表示個體y=0發(fā)生的概率。建模過程如下：

第一步，建立空模型，計算組內(nèi)相關(guān)系數(shù)（ICC）的值。空模型中僅有一個隨機(jī)截距而不包含任何自變量，其模型見式（2）。

式（2）中，j代表2水平單位的編號，C代表2水平單位的數(shù)目；i代表第j個2水平單位組中個體的編號，mj代表第j個2水平單位組中個體的數(shù)目，β0j代表第j個2水平單位組中的截距，其表達(dá)式見式（3）。

式（2）和式（3）兩個模型可合并成式（4）。

式（4）中，β0為y=0的總平均logit值，μ0j為組水平（本資料為中心）的平均logit變異值，它表示第j個組的平均logit值與總平均logit值之間的差異，且

多水平多重Logistic回歸模型的組間變異也可用ICC進(jìn)行評估，因Logistic回歸模型的殘差方差為π2/3，所以ICC的計算公式可由式（5）表示。

第二步，建立包含自變量的隨機(jī)截距模型，即在隨機(jī)截距的基礎(chǔ)上再考察變量drug的固定效應(yīng)。構(gòu)建的模型見式（6）。

式（6）中的β0j由式（7）給出。

式（6）和式（7）兩個模型可合并成式（8）的形式。

式（8）中，β0+ β1drugij為固定效應(yīng)，μ0j為隨機(jī)效應(yīng)，且

第三步，建立包含自變量的隨機(jī)截距-隨機(jī)斜率模型，即截距項和自變量drug的回歸系數(shù)中均含有隨機(jī)效應(yīng)項。構(gòu)建的模型見式（9）。

式（9）中等號右邊的β0j和β1j分別由式（10）和式（11）給出。

式（9）、式（10）、式（11）的3個模型可合并成式（12）的形式。

式（12）中，β0+β1drugij為 固 定 效 應(yīng) ，μ0j+μ1jdrugij為隨機(jī)效應(yīng)，且μ0j與μ1j之間的協(xié)方差可能有統(tǒng)計學(xué)意義；若無統(tǒng)計學(xué)意義，則將它們之間的協(xié)方差設(shè)定為0。

3 實例與SAS實現(xiàn)

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

【例1】某臨床研究中，研究者選擇16所醫(yī)院同時開展臨床試驗，每所醫(yī)院均選取受試者120人，在醫(yī)院內(nèi)隨機(jī)等分為兩組，分別接受試驗藥物和對照藥物的治療。治療結(jié)果見表1，試比較兩種藥物的療效［3］。

表1 多中心臨床試驗數(shù)據(jù)Table 1 Data from multicenter clinical trials

3.2 檢驗各中心優(yōu)勢比之間是否具有齊性

為了檢驗各中心優(yōu)勢比之間是否具有齊性，可利用 Breslow-Day檢驗［4］。設(shè)所需要的 SAS程序如下：

【SAS主要輸出結(jié)果及解釋】優(yōu)勢比齊性的Breslow-Day檢驗結(jié)果：χ2=37.994，df=15，P=0.000 9。此結(jié)果表明：16個中心的優(yōu)勢比不滿足齊性要求，說明不應(yīng)將16個中心的資料簡單合并后構(gòu)建多重Logistic回歸模型。

3.3 采用4種方法構(gòu)建多重Logistic回歸模型

3.3.1 對試驗中心產(chǎn)生啞變量后構(gòu)建多重Logistic回歸模型

設(shè)所需要的SAS程序如下：

【SAS主要輸出結(jié)果及解釋】對二重Logistic回歸模型中兩個參數(shù)進(jìn)行的“3型效應(yīng)分析”的結(jié)果見表2。

表2 模型中兩個參數(shù)的“3型效應(yīng)分析”的結(jié)果Table 2 Results of "Type 3 effect analysis" of two parameters in the model

由以上輸出結(jié)果可知，16個中心療效之間的差異有統(tǒng)計學(xué)意義，兩種藥物療效之間的差異也有統(tǒng)計學(xué)意義。因模型中參數(shù)估計部分的輸出結(jié)果很多，以下僅簡單呈現(xiàn)與第16中心比較差異有統(tǒng)計學(xué)意義的結(jié)果，見表3。

表3 與第16中心比較差異有統(tǒng)計學(xué)意義的結(jié)果Table 3 Statistically significant results compared with the 16th center

由以上輸出結(jié)果可知，僅5個中心與第16中心的療效之間的差異有統(tǒng)計學(xué)意義；試驗藥與對照藥療效之間的差異有統(tǒng)計學(xué)意義。由于程序中設(shè)定：擬合“y=0（治療成功）”的概率模型，并且，drug=0代表試驗藥，drug=1代表對照藥，drug的回歸系數(shù)為-0.826，即對照藥相對于試驗藥的優(yōu)勢比OR=exp（-0.826）=0.438。也就是說，試驗藥相對于對照藥的優(yōu)勢比為1/OR=1/0.438=2.283倍。

3.3.2 將試驗中心視為分層變量構(gòu)建多重Logistic回歸模型

相當(dāng)于在每個試驗中心內(nèi)部進(jìn)行配對設(shè)計，再進(jìn)行條件Logistic回歸分析。設(shè)所需要的SAS程序如下：

【SAS主要輸出結(jié)果及解釋】基于center分層的一重Logistic回歸模型的分析結(jié)果見表4。

表4 條件最大似然估計的分析結(jié)果Table 4 Analysis results of conditional maximum likelihood estimation

由以上輸出結(jié)果可知，參數(shù)drug的回歸系數(shù)為-0.819，即對照藥相對于試驗藥的優(yōu)勢比OR=exp（-0.819）=0.441。也就是說，試驗藥相對于對照藥的優(yōu)勢比為1/OR=1/0.441=2.268倍。

3.3.3 構(gòu)建隨機(jī)截距多水平多重Logistic回歸模型

從16個“不同中心”抽取受試對象，同一個中心的個體之間的差異被稱為1水平上的差異；而不同中心的個體之間的差異被稱為2水平上的差異，需要構(gòu)建隨機(jī)截距多水平多重Logistic回歸模型。設(shè)所需要的SAS程序如下：

【SAS主要輸出結(jié)果及解釋】與center對應(yīng)的截距項的計算結(jié)果：截距項中隨機(jī)部分V_u0=0.154（即隨機(jī)部分的方差）。模型中固定效應(yīng)的計算結(jié)果見表5。

表5 模型中固定效應(yīng)的計算結(jié)果Table 5 Calculation results of fixed effects in the model

由以上輸出結(jié)果可知，b0=0.530，b1=-0.814，它們與0之間的差異均有統(tǒng)計學(xué)意義。drug的回歸系數(shù)為-0.814，即對照藥相對于試驗藥的優(yōu)勢比OR=exp（-0.814）=0.443。也就是說，試驗藥相對于對照藥的優(yōu)勢比為1/OR=1/0.443=2.257倍。

為了檢驗截距項中隨機(jī)部分V_u0=0.154與0之間的差異是否有統(tǒng)計學(xué)意義，設(shè)所需要的SAS程序如下：

/*將以上模型估計的3個參數(shù)值代入下面的程序中*/

【SAS主要輸出結(jié)果及解釋】模型中各參數(shù)的估計結(jié)果見表6。

表6 模型中各參數(shù)的估計結(jié)果Table 6 Estimation results of parameters in the model

由以上輸出結(jié)果可知，對3個參數(shù)進(jìn)行了重新估計，并且，對它們都進(jìn)行了假設(shè)檢驗，P值均小于0.05。b0=0.533，b1=-0.819，V_u0=0.144。

3.3.4 構(gòu)建隨機(jī)截距-隨機(jī)斜率多水平多重Logistic回歸模型

還可以嘗試構(gòu)建隨機(jī)截距-隨機(jī)斜率多水平多重Logistic回歸模型。設(shè)所需要的SAS程序如下：

【SAS主要輸出結(jié)果及解釋】模型中參數(shù)的協(xié)方差參數(shù)估計結(jié)果見表7。

表7 模型中參數(shù)的協(xié)方差參數(shù)估計結(jié)果Table 7 Covariance parameter estimation results of parameters in the model

由以上輸出結(jié)果可知，截距中的隨機(jī)部分V_u0=0.432，斜率中的隨機(jī)部分V_u1=0.427，它們之間的協(xié)方差Cov_u01=-0.199。模型中參數(shù)的固定效應(yīng)的分析結(jié)果見表8。

表8 模型中參數(shù)的固定效應(yīng)的分析結(jié)果Table 8 Analysis results of fixed effects of parameters in the model

由以上輸出結(jié)果可知，截距和斜率的固定部分的估計值分別為b0=0.532，b1=-0.821，它們與0之間的差異均具有統(tǒng)計學(xué)意義。drug的回歸系數(shù)為-0.821，即對照藥相對于試驗藥的優(yōu)勢比OR=exp（-0.821）=0.440。也就是說，試驗藥相對于對照藥的優(yōu)勢比為1/OR=1/0.440=2.273倍。

為了檢驗3個隨機(jī)部分的估計值與0之間的差異是否有統(tǒng)計學(xué)意義，設(shè)所需要的SAS程序如下：

/*將以上模型估計的5個參數(shù)值代入下面的程序中*/

【SAS主要輸出結(jié)果及解釋】模型中參數(shù)的估計結(jié)果見表9。

表9 模型中參數(shù)的估計結(jié)果Table 9 Analysis results of fixed effects of parameters in the model

由以上輸出結(jié)果可知，后三行為隨機(jī)部分，它們的標(biāo)準(zhǔn)誤差無法計算出來，因此，無法檢驗它們與0之間的差異是否有統(tǒng)計學(xué)意義。從SAS日志可知，無法計算的原因是最終的海森矩陣不是正定矩陣，因此，估計的協(xié)方差矩陣不是滿秩的［4］。

3.4 結(jié)論

比較“第3.3.1節(jié)”到“第3.3.4節(jié)”的計算結(jié)果，可得出以下結(jié)論：第一，采用4種方法構(gòu)建多重Logistic回歸模型，所獲得的截距和斜率的固定部分的估計值比較接近；第二，采用多水平Logistic回歸模型分析，可以把不同中心資料優(yōu)勢比之間的不齊性以隨機(jī)效應(yīng)定量地呈現(xiàn)出來，即各中心試驗藥與對照藥的成功率之間存在明顯的跨中心效應(yīng)。也就是說，采用無層級結(jié)構(gòu)的多重Logistic回歸模型［5-7］分析層級之間存在明顯差異的統(tǒng)計資料是不合適的。

4 討論與小結(jié)

4.1 討論

由于本例資料中有一個多值名義變量“center（試驗中心）”，使統(tǒng)計分析者有多種可能的分析策略。除了本文已經(jīng)采用的4種建模方法之外，還有兩種誤用的做法：一是忽視變量“center”的存在；二是將變量“center”視為有序變量或定量變量（賦值為1，2，…，16）。

4.2 小結(jié)

本文介紹了與多水平模型分析有關(guān)的4個基本概念；介紹了構(gòu)建二水平Logistic回歸模型的三個步驟；針對一個多中心藥物臨床試驗的實例，介紹了用SAS實現(xiàn)分析的全過程，包括如何檢驗各試驗中心優(yōu)勢比的齊性，如何將試驗中心變換為啞變量后構(gòu)建多重Logistic回歸模型，如何將試驗中心視為分層變量構(gòu)建多重Logistic回歸模型，如何構(gòu)建隨機(jī)截距多重Logistic回歸模型，以及如何構(gòu)建隨機(jī)截距和隨機(jī)斜率多重Logistic回歸模型。對4種可行的建模方法所獲得的估計結(jié)果進(jìn)行了比較，作出了結(jié)論。最后，在討論中還指出了兩種可能被誤用的做法。