楊志雄 袁岱菁

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

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

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

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

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

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

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

模型與原理

1.非線性混合效應(yīng)模型

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

非線性混合效應(yīng)模型可作如下表述:

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

2.廣義線性混合效應(yīng)模型

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

其模型表述為:

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

條件均數(shù)uij(考慮了隨機效應(yīng))通過聯(lián)結(jié)函數(shù)g(．)與條件線性預(yù)測值ηij聯(lián)結(jié);

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

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

應(yīng)用實例

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

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

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

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

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

表1 非線性混合效應(yīng)模型參數(shù)估計表

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

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

表2 隨機效應(yīng)估計表

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

表3 固定效應(yīng)估計表

GLIMMX參數(shù)估計的解釋與NLMIXED類似。同時有注意到變量在GLIMMX中顯示了變量名，協(xié)方差參數(shù)1．0051與NLMIXED表中s2u等價?？梢钥闯?，雖然用NLMIXED過程步產(chǎn)生的參數(shù)估計值與GLIMMX結(jié)果不同，但是兩個方法得出的結(jié)論是一致的。

討 論

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

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

而廣義線性混合模型中隨機效應(yīng)也可以以非線性形式進入模型，反應(yīng)變量既可是連續(xù)變量，也可是分類變量，常用來處理相關(guān)數(shù)據(jù)、縱向觀測數(shù)據(jù)、過度離散的數(shù)據(jù)。隨機效應(yīng)變量之間既可以獨立也可以相關(guān)?？梢愿鶕?jù)需要擬合logistic隨機效應(yīng)模型，在醫(yī)學(xué)領(lǐng)域中同樣應(yīng)用廣泛。

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

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

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

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

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

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

1．王全眾．兩類分析相關(guān)數(shù)據(jù)的 logistic回歸模型．統(tǒng)計研究，2007，24(2):81-83．

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

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

4．羅天娥，劉桂芬．重復(fù)測量資料非線性混合效應(yīng)模型應(yīng)用與實現(xiàn)．中國衛(wèi)生統(tǒng)計，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．

附錄

用非線性混合效應(yīng)模型PROC NLMIXED分析數(shù)據(jù)，相應(yīng)的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;/*設(shè)定參數(shù)初值*/

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;

用廣義線性混合效應(yīng)模型PROC GLMMIX分析數(shù)據(jù)，相應(yīng)的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;