劉紅偉，張?zhí)鹛穑铋L平，2*，胡良平

（1. 天津醫(yī)科大學(xué)公共衛(wèi)生學(xué)院衛(wèi)生統(tǒng)計(jì)學(xué)教研室，天津 300070；2. 世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計(jì)學(xué)專業(yè)委員會，北京 100029；3. 軍事科學(xué)院研究生院，北京 100850

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

1 基本概念

1.1 二值資料

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

1.2 多水平數(shù)據(jù)

多水平數(shù)據(jù)或具有多水平層次結(jié)構(gòu)的數(shù)據(jù)是多水平統(tǒng)計(jì)模型發(fā)展和應(yīng)用的基礎(chǔ)。多水平數(shù)據(jù)也就是自然形成的層次數(shù)據(jù)，例如，在多中心臨床試驗(yàn)中，每個中心是水平2 單位，受試者是水平1 單位；在動物試驗(yàn)中，小鼠是水平1單位，窩別是水平2單位。多水平數(shù)據(jù)具有非獨(dú)立性，故無法采用廣義線性模型進(jìn)行分析，因此提出了能夠處理多水平數(shù)據(jù)的多水平模型［1-2］。

與廣義線性模型相比，多水平模型稍顯復(fù)雜，因?yàn)樗瑫r包含了多個水平的數(shù)據(jù)，從而在多個水平上都存在殘差?？傮w來說，其建模的思想就是把高水平上的差異估計(jì)出來（傳統(tǒng)的線性模型不考慮這一差異，將其放到了一個統(tǒng)一的殘差中），就使得殘差變小，估計(jì)的結(jié)果更可靠［3-4］。

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

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

【例1】某地區(qū)開展多中心臨床試驗(yàn)，擬比較兩種藥物治療某疾病的效果。數(shù)據(jù)見表1。

表1 多中心臨床試驗(yàn)數(shù)據(jù)

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

3.1 回歸模型的表達(dá)式

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

對于響應(yīng)變量為二值變量的層級數(shù)據(jù)結(jié)構(gòu)，可采用多水平logistic回歸模型分析，其模型可表達(dá)為：

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

3.2 回歸模型的構(gòu)建

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

第一步，建立空模型，計(jì)算組內(nèi)相關(guān)系數(shù)ICC的值?？漳Ｐ椭袃H有一個隨機(jī)截距而不包含任何解釋變量，其模型為：，上述模型可合并為：

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

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

其中，β0j+β1drugij為固定效應(yīng)，μ0j為隨機(jī)效應(yīng)，且。

第三步，建立包含解釋變量的隨機(jī)截距-斜率模型，即截距項(xiàng)和解釋變量drug 的系數(shù)均為隨機(jī)系

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

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

4 基于SAS分析實(shí)例

4.1 分析與解答

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

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

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

【主要輸出結(jié)果及解釋】以下是第一個過程步輸出結(jié)果，即調(diào)用GLIMMIX過程運(yùn)行空模型。模型構(gòu)建是以”y=0”為基礎(chǔ)的，即計(jì)算”y=0”發(fā)生的概率模型。

以上是模型擬合的有關(guān)信息，第一行即為-2 倍的限制性∕殘差虛擬對數(shù)似然值，此統(tǒng)計(jì)量不能用于不同模型的比較。

以上是協(xié)方差參數(shù)估計(jì)的結(jié)果，給出了隨機(jī)效應(yīng)方差的估計(jì)值及相關(guān)假設(shè)檢驗(yàn)的結(jié)果?？梢婋S機(jī)截距方差（即σ2μ0）的估計(jì)值為0.5988，標(biāo)準(zhǔn)誤為0.2686。但此處未給出隨機(jī)截距方差是否為0的假設(shè)檢驗(yàn)結(jié)果。故沒有客觀依據(jù)判定σ2μ0與0 之間的差異是否有統(tǒng)計(jì)學(xué)意義。

以上是固定效應(yīng)的解。因?yàn)榇诉^程步運(yùn)行的是空模型，所以這里只有一個固定效應(yīng)，即截距，其值為-0.3018，表示y=0的總平均logit值為-0.3018。

以下是模型擬合的有關(guān)信息，包括三種信息標(biāo)準(zhǔn)的估計(jì)值和-2 倍的對數(shù)似然值。這些統(tǒng)計(jì)量本身不能說明模型擬合的優(yōu)劣，但可用于含相同自變量數(shù)目的不同模型的比較。

以上是模型中參數(shù)的估計(jì)結(jié)果，包括固定效應(yīng)和隨機(jī)效應(yīng)的參數(shù)估計(jì)值及相應(yīng)的假設(shè)檢驗(yàn)結(jié)果。注意，隨機(jī)效應(yīng)假設(shè)檢驗(yàn)給出的是雙側(cè)檢驗(yàn)的結(jié)果，而實(shí)際上檢驗(yàn)方差是否為0應(yīng)采用單側(cè)檢驗(yàn)，故此處所得的P 值應(yīng)除以2 才是正確的P 值，后同?！皏_u0”對應(yīng)的P值為0.0399∕2＜0.05，說明σ2μ0與0之間差異有統(tǒng)計(jì)學(xué)意義，分析時應(yīng)采用多水平logistic回歸模型分析。

以上是ICC的計(jì)算結(jié)果，其值為0.1487，對應(yīng)的P值為0.0187＜0.05，說明數(shù)據(jù)存在一定的組內(nèi)同質(zhì)性，需采用多水平logistic模型分析該資料。

以下是第三個過程步的輸出結(jié)果，即調(diào)用GLIMMIX 過程運(yùn)行含解釋變量drug的隨機(jī)截距模型的結(jié)果，其截距項(xiàng)設(shè)定為0?？梢?，σ2μ0值為0.5984，解釋變量系數(shù)為-0.4343。

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

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

5 討 論

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

采用PROC GLIMMIX 和PROC NLMIXED 過程來構(gòu)建模型：首先建立不包含任何解釋變量的空模型，以計(jì)算ICC值。若存在組內(nèi)相關(guān)，則構(gòu)建截距項(xiàng)不為0 的模型。若經(jīng)檢驗(yàn)得到此截距項(xiàng)與0 的差異沒有統(tǒng)計(jì)學(xué)意義，則構(gòu)建截距項(xiàng)為0的模型。