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

劉紅偉，張甜甜，李長平，2*，胡良平

（1. 天津醫科大學公共衛生學院衛生統計學教研室，天津 300070；2. 世界中醫藥學會聯合會臨床科研統計學專業委員會，北京 100029；3. 軍事科學院研究生院，北京 100850

*通信作者：李長平，E-mail：1067181059@qq.com）

1 基本概念

1.1 二值資料

結局變量只有兩個取值的資料稱為“二值資料”，例如，在表1多中心臨床試驗數據中，治療結局取值為“成功”或“失敗”。

1.2 多水平數據

多水平數據或具有多水平層次結構的數據是多水平統計模型發展和應用的基礎。多水平數據也就是自然形成的層次數據，例如，在多中心臨床試驗中，每個中心是水平2 單位，受試者是水平1 單位；在動物試驗中，小鼠是水平1單位，窩別是水平2單位。多水平數據具有非獨立性，故無法采用廣義線性模型進行分析，因此提出了能夠處理多水平數據的多水平模型［1-2］。

與廣義線性模型相比，多水平模型稍顯復雜，因為它同時包含了多個水平的數據，從而在多個水平上都存在殘差。總體來說，其建模的思想就是把高水平上的差異估計出來（傳統的線性模型不考慮這一差異，將其放到了一個統一的殘差中），就使得殘差變小，估計的結果更可靠［3-4］。

雖然理論上多水平模型可以有多個層次，但實際中最常用的是二水平模型，因此這里主要通過一份二水平數據簡要介紹多水平模型構建與求解的思路。

2 數據結構

【例1】某地區開展多中心臨床試驗，擬比較兩種藥物治療某疾病的效果。數據見表1。

表1 多中心臨床試驗數據

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

3.1 回歸模型的表達式

對于響應變量為二值變量的非層級結構數據，一般采用普通logistic 回歸模型分析，又稱為固定效應logistic回歸模型分析。設P（y=1|X）（簡記為P）表示暴露因素為X 時個體發生陽性事件（以y=1 表示發生陽性事件）的概率，而陽性事件發生的概率P與陰性事件發生的概率（1-P）之比稱為優勢比。對優勢比進行自然對數變換即為對P的logit變換，得：

對于響應變量為二值變量的層級數據結構，可采用多水平logistic回歸模型分析，其模型可表達為：

其中，U是隨機回歸系數向量，服從N(0，G)，G為協方差矩陣，β是水平1固定回歸系數向量，X和Z分別是固定效應和隨機效應的解釋變量設計矩陣。

3.2 回歸模型的構建

以例1 為例，以變量zhongxin 表示“中心”，以變量drug 表示“藥物種類”，以變量y 表示療效（y=0 表示治療成功，y=1表示治療失敗），以Pij表示個體y=0發生的概率。建模過程如下：

第一步，建立空模型，計算組內相關系數ICC的值。空模型中僅有一個隨機截距而不包含任何解釋變量，其模型為：，上述模型可合并為：

其中，βo為y=0 的總平均logit 值，μ0j為組水平（本資料為中心）的平均logit 值的變異量，表示第j個組的平均logit值與總平均logit值之間的差異，且。多水平logistic 回歸模型的組間變異也可用組內相關系數進行評估，因logistic 回歸模型的殘差方差為π2∕3，所以：

第二步，建立包含解釋變量的隨機截距模型，即在隨機截距的基礎上再考察變量drug 的固定效應。模型如下，上述模型可合并為：

其中，β0j+β1drugij為固定效應，μ0j為隨機效應，且。

第三步，建立包含解釋變量的隨機截距-斜率模型，即截距項和解釋變量drug 的系數均為隨機系

其中，β0+β1drugij為固定效應，μ0j+μ1jdrugij為隨機效應，且，μ0j與μ1j之間的協方差可能有統計學意義；若無統計學意義，則將它們之間的協方差設為0。

例1較為特殊，其最終模型是包含解釋變量drug的隨機截距模型，但又與第二步所建模型略有不同，區別在于本資料截距項為0，即βo=0。例1 的最終模型中包含一個固定效應和一個隨機效應，模型如下：

4 基于SAS分析實例

4.1 分析與解答

在例1中，研究者欲考察試驗藥與對照藥治療某種感染的效果。資料中涉及兩個原因變量——中心和藥物種類，響應變量為二值變量，由于不同醫院對同一疾病的治療效果可能有差異，而在同一醫院中，相同疾病的治療效果也并不完全獨立。故可考慮采用多水平logistic回歸模型分析，SAS程序如下：

【說明】程序共5步，包括1個數據步和4個過程步。程序第2 步、第3 步是建立不包含任何解釋變量的空模型，以計算ICC 值。程序第4、5 步均是建立包含解釋變量drug 的隨機截距模型。GLIMMIX過程的計算結果與NLMIXED 過程的結果會略有差異，前者運算速度快、用法簡單，但在評估模型擬合效果時使用虛擬的對數似然值，而非真實值，不能用于模型的比較，且SAS 9.4 的中GLIMMIX 過程沒有提供隨機效應的假設檢驗，其結果雖有隨機系數方差的參數估計值及標準誤，但兩者的比值只能作為參考，不能采用t 檢驗計算相應的P 值。NLMIXED過程可以提供真實的對數似然值，并為隨機效應提供假設檢驗的結果，也可以通過似然比檢驗對嵌套模型的擬合效果進行比較，但其用法較復雜，需設置模型和參數的初始值，不便于使用。因此一般以GLIMMIX 過程得到的參數估計值作為NLMIXED 過程的模型參數初始值，最后以NLMIXED 過程的結果為準。對于相對簡單的模型而言，NLMIXED 過程對參數初始值并不敏感，此時采用其默認的初始值1 即可。調用NLMIXED 過程運行包含解釋變量的隨機截距-斜率模型，所用程序與本節程序第5 步有較大修改。參考程序如下：

其中，“b0”“b1””v_u0”“cov_u01”“v_u1”分別相當于公式（7）中β0、β1、μ0j的方差，μ0j與μ1j之間的協方差，μ1j的方差。

【主要輸出結果及解釋】以下是第一個過程步輸出結果，即調用GLIMMIX過程運行空模型。模型構建是以”y=0”為基礎的，即計算”y=0”發生的概率模型。

以上是模型擬合的有關信息，第一行即為-2 倍的限制性∕殘差虛擬對數似然值，此統計量不能用于不同模型的比較。

以上是協方差參數估計的結果，給出了隨機效應方差的估計值及相關假設檢驗的結果。可見隨機截距方差（即σ2μ0）的估計值為0.5988，標準誤為0.2686。但此處未給出隨機截距方差是否為0的假設檢驗結果。故沒有客觀依據判定σ2μ0與0 之間的差異是否有統計學意義。

以上是固定效應的解。因為此過程步運行的是空模型，所以這里只有一個固定效應，即截距，其值為-0.3018，表示y=0的總平均logit值為-0.3018。

以下是模型擬合的有關信息，包括三種信息標準的估計值和-2 倍的對數似然值。這些統計量本身不能說明模型擬合的優劣，但可用于含相同自變量數目的不同模型的比較。

以上是模型中參數的估計結果，包括固定效應和隨機效應的參數估計值及相應的假設檢驗結果。注意，隨機效應假設檢驗給出的是雙側檢驗的結果，而實際上檢驗方差是否為0應采用單側檢驗，故此處所得的P 值應除以2 才是正確的P 值，后同。“v_u0”對應的P值為0.0399∕2＜0.05，說明σ2μ0與0之間差異有統計學意義，分析時應采用多水平logistic回歸模型分析。

以上是ICC的計算結果，其值為0.1487，對應的P值為0.0187＜0.05，說明數據存在一定的組內同質性，需采用多水平logistic模型分析該資料。

以下是第三個過程步的輸出結果，即調用GLIMMIX 過程運行含解釋變量drug的隨機截距模型的結果，其截距項設定為0。可見，σ2μ0值為0.5984，解釋變量系數為-0.4343。

以下是第四個過程步的輸出結果，即調用NLMIXED 過程運行含解釋變量drug 的隨機截距模型的結果，其截距項設定為0。可見，解釋變量 系 數 為-0.4409（P=0.0057＜0.05），σ2μ0值 為0.6221（P=0.0390＜0.05），二者與0 的差異均有統計學意義。另外，由“Fit Statistics”部分結果可知，-2 倍對數似然值為1225.1，略大于第四步構造的模型；但由“Parameter Estimates”部分結果可知，模型中共包含兩個參數，參數個數較之前少了一個，且兩個模型的擬合效果并無統計學差異（χ2=1225.1-1224.9=0.2，P＞0.05）。所以，使用此模型更合適。

解釋變量drug的系數為-0.4409，且與0的差異有統計學意義（P=0.0057），說明試驗藥組與對照藥組的療效之差具有統計學意義。因exp（-0.4409）=0.6435，所以對照藥組治療成功率是試驗藥組治療成功率的0.6435 倍。隨機截距的方差“v_u0”估計值為0.6221，與0的差異有統計學意義（P=0.0390∕2＜0.05），說明水平1截距跨中心變異顯著，即不同中心μ0j值存在差異。

5 討 論

一般而言，響應變量為二值變量的高維列聯表資料采用一般logistic 回歸分析，但此法要求所有觀測結局相互獨立。對于研究個體存在聚集性特征時，應采用多水平模型。這樣可將傳統模型中的隨機誤差分解到數據層級結構相應的水平上，使得個體的隨機誤差更純［5］。

采用PROC GLIMMIX 和PROC NLMIXED 過程來構建模型：首先建立不包含任何解釋變量的空模型，以計算ICC值。若存在組內相關，則構建截距項不為0 的模型。若經檢驗得到此截距項與0 的差異沒有統計學意義，則構建截距項為0的模型。