非線性混合效應模型和廣義線性模型擬合隨機效應logistic回歸的應用比較

楊志雄 袁岱菁

在臨床藥物試驗中藥物療效的評價經常遇到二分類資料，即反應變量有兩個水平如有效、無效;成功、失敗等。二分類變量服從二項分布，可采用logistic回歸模型。

運用logistic回歸模型對分類資料進行分析，能給實際研究帶來很多便利。與多元線性回歸相比，logistic回歸具有許多獨特的優點，如對正態性和方差齊性不做要求，系數的可解釋性等。正是這些優點，使得logistic回歸成為廣受歡迎的分析工具。但是需要指出的是，logistic回歸模型只能處理具有獨立性的資料，即觀測數據應來自完全獨立的隨機樣本。

在實際研究中經常會碰到樣本之間具有相互關系的觀測數據(correlated data)。比如，臨床研究中的交叉設計試驗，在不同試驗條件下對同一對象的重復觀測。

由于重復測量數據間存在自相關性，因而增加了傳統統計方法對該類數據分析的難度。適合此類數據的統計分析方法就必須考慮數據之間的這種相關性。這些方法大多是傳統logistic回歸的擴展，大致可歸為兩類:邊際模型(marginal model)和隨機效應模型(random effect model)〔1〕。

對于交叉設計試驗中，在不同試驗條件下對同一對象(受試者)的重復觀測的數據，考慮到同一受試者在不同處理，不同階段下測量數據之間的相關關系是由受試者的內在特性引起的，而這種特性在不同受試者之間是不相同的，且無法實際觀測到，所以適用隨機效應的logistic回歸模型。

我們既可以用非線性混合效應模型，也可以用廣義線性模型來擬合隨機效應的logistic回歸。SAS提供了不同的過程步來實現非線性混合效應模型和廣義線性模型，分別是PROC NLMIXED，PROC GLMMIX。

下面就用一個實例來介紹并比較用兩種模型擬合交叉設計試驗下隨機效應logistic回歸。

模型與原理

1.非線性混合效應模型

非線性混合效應模型亦稱為多水平非線性模型、非線性隨機效應模型或非線性分層模型。它可以直接擬合非線性模型。不僅能識別和估計個體間和個體內的變異，而且也考慮了解釋變量與反應變量參數間的非線性關系，允許固定效應和隨機效應進入模型的非線性部分，相對于線性模型的正態假定，非線性模型對資料的分布無特殊要求，資料可以是正態，也可以是二項分布、Poisson分布等。

非線性混合效應模型可作如下表述:

其中，yij為第i個體第j次測量預測值，或經過某種單調聯系函數(1ink function)轉換的值;f(·)為非線性函數，如果其為線性，則退化為線性的混合效應模型;xij為P維解釋變量向量;eij為獨立正態分布隨機誤差向量;β為P維固定效應參數;bi為隨機效應因子;Ai、Bi為已知的設計矩陣。其參數估計可以通過偽數據步(pseudo—data step)和線性混合效應步(1inear mixed effects step)兩步之間的迭代完成，可分別使用Gauss—Newton迭代法和 EM 算法解決〔2〕。

2.廣義線性混合效應模型

廣義線性混合模型(gneralized linear mixed models，GLMMs)是廣義線性模型(generalized linear models，fGLMs)與線性混合效應模型(1inear mixed model)的擴展．通過在模型中納入隨機效應來解釋數據間的相關，過度離散(overdispersion)、異質性(heterogeneity)等問題。

其模型表述為:

ui為隨機效應項，反應變量Yij的條件分布的期望;

條件均數uij(考慮了隨機效應)通過聯結函數g(．)與條件線性預測值ηij聯結;

上式為廣義線性混合模型的一般式，Yij:第i類的第 j個觀測的反應變量，i=1，…，m，j=l，…ni。考慮隨機效應ui的條件下獨立，服從指數分布族，可以是二項分布、Poisson分布、Gamma分布等。

Xij為解釋變量，β為固定效應參數矢量，ui隨機效應，服從均數為0，方差協方差矩陣為γ的多變量正態分布，ui解釋了由于不可測因子引起的類間的異質性和同一類內觀測到的相關，不同類間的ui是相互獨立的，Zij為與隨機效應相關的解釋變量。設計矩陣分固定效應X與隨機效應Z兩部分。分析的數據不同可以選擇不同的聯結函數g(·)，可以擬合含隨機效應的 logistic回歸等多種模型〔3〕。

應用實例

某公司開發一種新藥用于治療泌尿系統疾病。以安慰劑為對照，將所有患者隨機分為兩組，一組患者先服用新藥A，再服用對照藥B;另一組患者順序相反，即先服用對照藥B，再服用新藥A，每個階段用藥2周，期間洗脫期1周，進行2×2交叉設計試驗，共納入病例30例。主要結果指標療效為每階段結束后，病人填寫PRTI量表。該量表由三個問題組成:包括病人對治療的總的滿意度，是否優先選擇研究中使用的治療和是否愿意再使用研究藥物。每一個問題上的得分是1～5分。其中1分為最滿意，5分為最不滿意。為了便于分析，按照PRTI量表評分的情況將療效劃分的2個等級，即有效和無效，分別賦值為1，0。當評分為1分和2分時，認為有效，賦值為1。當評分為3分，4分和5分時，認為無效，賦值為0。這樣三個問題中每個問題都會產生一個二分類變量。因篇幅有限，本例中的結果僅針對于對第一個問題即病人對治療的總的滿意度的分析。

統計分析方法以病人對治療的總的滿意度即藥物是否有效為應變量，處理(treat)、階段(period)、受試者(subject)為解釋變量，其中處理和階段納入為固定效應，受試者(subject)作為隨機效應，建立非線性混合效應模型。其模型基本形式為:

隨機效應 ui～Normal(0，δ2)。

所有分析過程在SAS軟件中完成。

用非線性混合效應模型PROC NLMIXED分析數據，所建立的模型在經過10次迭代后收斂。采用常用的Dual quasi-Newton最優化技術與Adaptive Gaussian積分方法，－2Loglikelihood=77.8，參數估計結果表1。

表1 非線性混合效應模型參數估計表

本表列出了4個參數和它們的最大似然估計值，標準誤，以及統計推斷。beta0是截距，表示處理和階段效應為0時的對數優勢(log-odds)。接下來幾行分別是處理(beta1)，階段(beta2)，隨機效應(s2u)的估計。每個系數都可以通過取冪來轉換成優勢比。檢驗水準α=0.05處理效應，階段效應和受試者效應均無統計學意義，如果某項效應有統計學意義，解釋為，通過取冪可以計算在控制了有關混雜因素后該效應的優勢比。

用廣義線性混合效應模型PROC GLMMIX分析數據，所建立的模型在經過8次迭代后收斂。采用常用的Newton-Raphson with Ridging最優化技術，－2 Res Log Pseudo-Likelihood=259.04，參數估計結果見表2和表3。

表2 隨機效應估計表

表3為模型中固定效應的解決方案。結果顯示處理效應和階段效應研究藥物組和安慰劑組相比兩組間比較均無統計學意義差異。

表3 固定效應估計表

GLIMMX參數估計的解釋與NLMIXED類似。同時有注意到變量在GLIMMX中顯示了變量名，協方差參數1．0051與NLMIXED表中s2u等價。可以看出，雖然用NLMIXED過程步產生的參數估計值與GLIMMX結果不同，但是兩個方法得出的結論是一致的。

討 論

醫學研究中常會遇到非線性重復觀測的二項分布數據分析的問題。對于二項分布數據的分析，常常采用logistic回歸模型，但是標準logistic回歸模型要求數據滿足獨立性的要求。重復觀測的數據不獨立。如本例中交叉設計的臨床試驗數據，數據間存在自相關性且隨機誤差至少分為兩個層次，即個體間誤差和個體內重復測量誤差。分析此類數據不僅需要考慮不同層次的誤差，而且也需要考慮參數間的非線性關系〔4〕。

非線性混合效應模型考慮了不同層次的誤差和參數間的非線性關系。允許固定效應和隨機效應進入模型的非線性部分，可以擬合具有隨機效應的logistic回歸模型，在臨床研究領域有著越來越廣泛的應用。

而廣義線性混合模型中隨機效應也可以以非線性形式進入模型，反應變量既可是連續變量，也可是分類變量，常用來處理相關數據、縱向觀測數據、過度離散的數據。隨機效應變量之間既可以獨立也可以相關。可以根據需要擬合logistic隨機效應模型，在醫學領域中同樣應用廣泛。

實例中，采用SAS軟件中的PROC NLMIXED和PROC GLMMIX過程來分析數據。比較NLMIXED和GLIMMIX，二者有幾處重要的不同，在實際應用中選擇哪個應仔細考慮。

首先是二者使用的估計方法不同。二者都利用了參數估計的方法。NLMIXED利用高斯積分來積分近似，而GLIMMIX依賴線性混合模型進行參數估計(線性法)。每個方法都有其優缺點。NLMIXED的優點是它更加準確，產生真實的對數似然擬合模型。用戶可以有相當大的彈性去定義似然函數。與之相反的是，GLIMMIX對固定效應和協方差參數的估計可能不是無偏估計，尤其是當所分析的數據是二項分布時。

GLIMMIX進行Wald-type檢驗和生成相應可信區間，嵌套的模型不及真實似然比檢驗精確。估計不精確換來的是GLIMMIX可以擬合不同類型的模型，隨機效應的數量可以有多個，也可以設置擬合的選項。例如，GLIMMIX允許多個嵌套的，交叉的隨機效應，但是NLMIXED所能處理的隨機效應少于5個。

此外，GLIMMIX允許用戶使用REML(restricted maximum likelihood)的方法，NLMIXED過程則不能用REML，只能用 ML(maximum likelihood)。GLIMMIX也支持以模型為基礎的對標準誤的Sandwich estimation，而NLMIXED只支持基于模型的標準誤。即使方差函數錯誤，Sandwich estimation也能提供一致的結果。但是，正如前面所提及的，GLIMMIX要求反應變量服從指數族分布，而NLMIXED允許用戶寫出自己的函數。

這兩個過程步另外一個不同是初值的設定和應用。NLMIXED，用戶必須先產生參數的初值然后把初值代入SAS程序中。一般地，可以使用PROC MIXED或PROC GENMOD來產生NLMIXED的初值。與之對照的是，GLIMMIX使用的是一個雙迭代方案。參數的初值來自于線性混合模型，而此初值又用來擬合線性模型。對于用戶而言，GLIMMIX不需要費力去設定初值。因此在某些情況下，對模型的控制缺乏也帶來一些不足。最后的參數估計對初值非常敏感，如果設定不好，很可能會導致模型不收斂。NLMIXED過程可以允許用戶去定義初值，因而也給模型的收斂提供了更多的機會。積累初始值的選擇經驗、判斷是否滿足收斂條件及模型評價的標準是非線性分析的關鍵〔5，6〕。

總之，NLMIXED更適合相對簡單，反應變量只有兩個水平，隨機效應的數量少的情況。當二分類數據需要精確的協方差參數估計，要求用戶自定義響應分布，或需要使用似然比檢驗比較嵌套模型時，可以考慮使用NLMIXED。NLMIXED更適合二項分布數據分析。GLIMMIX則更適合隨機效應數量多于兩個以上的復雜模型。本例中用NLMIXED分析數據顯然更加適合。

1．王全眾．兩類分析相關數據的 logistic回歸模型．統計研究，2007，24(2):81-83．

2．陳衛中，楊曉虹，陳朝瓊，等．非線性混合效應模型在交叉設計等級資料分析中的應用．成都醫學院學報，2007，2(3-4):181-183．

3．李麗霞，郜艷暉，張丕德，等．廣義線性混合效應模型及其應用．現代預防醫學，2007，34(11):2103-2104．

4．羅天娥，劉桂芬．重復測量資料非線性混合效應模型應用與實現．中國衛生統計，2006，23(2):104-107．

5．Flom PL，McMahon JM，Pouget ER．Using PROC NLMIXED and PROC GLMMIX to analyze dyadic data with a dichotomous dependent variable．SASGlobal Forum 2007 Proceedings:Paper 179．

6．SAS Institute Inc．SAS/STAT ? 9.2 User’s Guide．Cary，NC，USA:SAS Institute Inc．，2008，4337-4340．

附錄

用非線性混合效應模型PROC NLMIXED分析數據，相應的SAS程序為:

proc nlmixed data=mpr3 alpha=0.05 corr gconv=1e-10;

parms beta0=0.3855 beta1=－0.1599 beta2=0.07278 s2u=0.1709;/*設定參數初值*/

eta=beta0+beta1*treat+beta2*period+u;

expeta=exp(eta);

p=expeta/(1+expeta);

model question1～binary(p);

random u～normal(0，s2u)subject=patient;

ods output parameterestimates=para;

run;

用廣義線性混合效應模型PROC GLMMIX分析數據，相應的SAS程序為:

proc glimmix data=mpr3 order=internal;

class patient treat period;

model question1(ref=first)=treat period/solution link=logit dist=binary;

random intercept/subject=patient gcorr;

lsmeans treat/diff=control(“－0.5”)oddsratio cl alpha=0.05;

lsmeans period/diff=control(“－0.5”)oddsratio cl alpha=0.05;

ods output parameterestimates=solf1;

run;