一類線性混合效應模型分位數估計的影響分析

李宏凱 肖松濤 歐陽應根 李志強*

(1.北京化工大學 數理學院， 北京 100029； 2.中國原子能科學研究院 放射化學研究所， 北京 102413)

引 言

線性混合效應模型是處理縱向數據最主要的模型之一，Diggle等[1]首次基于線性混合模型詳細地研究了針對縱向數據的統計分析方法。隨著縱向數據研究的不斷發(fā)展，線性混合模型在生物、醫(yī)學、衛(wèi)生、心理學、農學和經濟學等各個領域都有了非常廣泛的應用。

影響分析作為統計建模中非常重要的一部分，主要被應用于識別數據中的強影響點和異常值點。近年來，關于線性混合模型的影響分析問題已經成為了研究熱點之一。例如，當總體誤差服從正態(tài)分布時，Pan等[2]基于似然函數及其導出的Q函數研究了線性混合模型在不同的隨機效應協方差結構下的影響診斷問題；當模型誤差服從指數族分布時，Xu等[3]基于數據刪除模型對廣義線性混合效應模型進行影響分析，基于似然函數關于隨機效應的條件期望構造了Q函數，進而建立了廣義Cook距離和似然距離等；Pinho等[4]考慮了廣義線性混合模型不同參數的影響分析問題，并且基于固定效應和隨機效應估計的聯合影響構建了廣義Cook距離。此外，孫慧慧等[5]與Tang等[6]還分別研究了線性混合模型和部分線性混合模型的局部影響分析問題。然而，在很多情況下縱向數據模型的誤差可能服從非正態(tài)分布或者非對稱分布，此時一般的估計精度很容易受到影響[4]。分位數回歸估計作為一種穩(wěn)健的估計方法，不僅可以克服上述問題，并且能夠獲得數據的總體分布信息，尤其是當回歸誤差項服從重尾分布或其分布受到污染時，分位數回歸估計具有更高的穩(wěn)健性。

Koenker[7]最早對混合效應模型的分位數估計方法進行了研究，他利用加權分位數和對隨機效應進行L1懲罰的方式，給出了模型的損失函數，獲得了模型系數和隨機效應的分位數估計。Geraci等[8]利用Monte Carlo expectation maximization(MCEM)算法消除了隨機效應的影響，并求解了線性混合模型系數的分位數估計。盡管分位數回歸估計具有穩(wěn)健性，但是Narula等[9]的研究表明，在分位數τ較大或者較小的情況下，估計結果也會受到極端異常點的影響，此外他們的數值模擬結果和實際數據分析也說明了針對最小絕對偏差(least absolute deviation，LAD)回歸進行影響分析等的統計診斷技術對異常數據的檢測是非常重要的。盡管目前已經有很多關于影響分析的研究，但由于計算及統計推斷的復雜性，針對線性混合效應模型分位數估計的影響分析的文獻還很少見到，因此基于線性混合效應模型分位數估計對強影響點的診斷度量和異常值的檢測進行研究是十分必要的。

本文利用線性混合效應模型的分位數估計處理服從非對稱、非正態(tài)，尤其是重尾或輕尾分布的縱向數據模型的影響分析問題。為了檢驗強影響點并減少計算量，通過對分位數回歸的光滑逼近得到了數據刪除模型參數的一步近似估計，并構造了檢驗強影響點的損失函數距離和Cook距離。在識別異常值時，首先證明了數據刪除模型和均值漂移模型參數估計的等價性，然后基于參數等價性給出計算Wald統計量的方法，并利用Bootstrap方法得到異常值檢驗的拒絕域。

1 實驗數據描述及線性混合模型

1.1 數據描述

為了對氨基羥基脲(HSC)水溶液及硝酸水溶液在γ射線作用下的輻解產物進定量分析，考察不同輻解劑量、HSC初始濃度以及硝酸濃度條件下主要輻解產物的生成量，分別配制濃度為0.1、0.2、0.4、0.6、0.8、1.0 mol/L的HSC溶液和濃度為0.2、0.4、0.6、0.8、1.0、1.5 mol/L的硝酸溶液，兩兩組合得到36種不同濃度配比的溶液，委托中國原子能科學研究院輻照中心，采用實驗型鈷源裝置(60Co源裝置)，分別置于輻解劑量為1×103、2×103、4×103、6×103、8×103、10×103Gy的條件下，進行三因素下的無重復實驗，總共得到216個樣品。

為了考察不同條件下輻解劑量(103Gy) 、HSC初始濃度(mol/L)以及硝酸濃度(mol/L)對銨根離子濃度(mmol/L)變化情況的影響，首先對數據進行初步的相關分析和散點圖對比，結果表明輻解劑量與銨根離子濃度大體呈線性關系，而單獨的HSC初始濃度和硝酸濃度對銨根離子濃度的影響則無明顯的趨勢變化規(guī)律。因此，本文將輻解劑量作為重要的解釋變量，而將硝酸濃度和HSC濃度的聯合影響作為銨根離子濃度的隨機影響因素，建立線性混合效應模型，分別考察兩個因素的影響大小。

1.2 線性混合模型及其分位數估計

(1)

由于很多化學實驗數據通常服從非對稱、非正態(tài)的概率分布，為了克服誤差分布的影響，采用穩(wěn)健的分位數估計方法。根據Koenker[7]的定義，在τ分位數下的線性混合模型具有如下形式

(2)

為了解決估計隨機效應時可能會面臨的維數較高的問題，Koenker[7]提出一種基于隨機效應的L1懲罰方法，通過加權多個分位數下的信息來獲取參數與隨機效應的聯合估計，其損失函數定義如下

(3)

式中，ρτ(r)=rτI(r>0)+r(τ-1)I(r<0)，q是分位數個數，wk是第k個分位數的權重，λ是懲罰因子。Koenker[7]建議將懲罰因子λ設置為λ=σε/σb，其中σε和σb可通過計算軟件進行求解，具體選取過程詳見第3節(jié)。

由文獻[7]可知，由于隨機效應服從零均值的正態(tài)分布，為了保證固定效應的估計量具有較好的統計性質，要求固定效應中必須包含截距項。通過極小化損失函數(式(3))獲得的參數估計記為=(,)T。為求解損失函數的極小化解，以往的學者們提出了不同的求解思想和算法。其中，Koenker[7]將問題轉化為線性規(guī)劃問題，然后利用單純形法或內點法等優(yōu)化算法進行計算。針對式(3)的損失函數，R語言提供了regression quantiles for panel data (RQPD)程序包用于求解縱向數據的分位數估計。

2 線性混合效應模型分位數估計的影響分析

實際數據中常常會存在一些異常數據，通常將對回歸估計影響過大的數據點稱為強影響點，而將偏離回歸直線較遠的點稱為異常值點。有時強影響點也很有可能同時是異常值點。影響分析正是為了識別出這些數據點而發(fā)展起來的一個統計分支，影響分析中常用的診斷模型包括數據刪除模型(case deletion model，CDM)和均值漂移模型(mean shift outlier model，MSOM)。

2.1 數據刪除模型的影響分析

2.1.1數據刪除模型

數據刪除模型主要通過依次刪除每個數據點來求解參數，之后再對比各個參數估計，若數據中存在較強的影響點，那么當刪除該點時，模型的參數估計相對于原估計就會有較大的變化。

假定刪去第(i1,j1)個觀測值時得到的模型為

(4)

2.1.2參數一步(光滑)近似估計

模型(4)相較原模型(1)只是減少了1個數據，所以參數估計的方法及計算復雜度與原模型基本相同，但是當數據量很大時，由于需要逐一檢驗每個數據點，將會產生很大的計算負擔。針對線性模型，Cook等[10]提出了利用似然函數的Taylor展開式構造參數近似估計的思想。

為了能夠對不可導的損失函數(3)進行光滑逼近，以建立數據刪除模型參數的一步近似估計，本文結合一般分位數參數估計的minorize- maximization(MM)算法[11]和絕對值函數的光滑方法[12]，構造如下目標函數

(5)

(6)

對式(6)兩端關于η[i1j1]求導可得[10]

(7)

通常情況下，參數估計都是多維向量，不便于直接比較，所以需要通過構建診斷度量來衡量參數變化的大小。在普通的回歸模型中，常用的診斷度量有似然距離和Cook距離等，本文將似然距離推廣到基于損失函數的距離。

1)根據損失函數(3)，構建基于損失函數的距離

LD=2[L(,)-L[i1j1]([i1j1],[i1j1])]

(8)

式中，L[i1j1]([i1j1],[i1j1])是刪去第(i1,j1)個觀測值的數據刪除模型的損失函數。該距離類似于似然距離，但是使用的目標函數不同。另外針對損失函數(3)中的權重選取，可以利用損失函數距離(8)分別研究數據在不同分位數下的情況。

2)Cook距離定義為[10]

LCD=([i1j1]-)T[Var()]-1([i1j1]-)/p

(9)

式中p為參數η的維數。由于沒有對模型(2)的誤差項作任何的分布假定，所以在中小樣本的情況下很難求得Cook距離中的Var()。為了解決此問題，本文在中小樣本條件下采用Bootstrap方法來近似計算樣本方差Var()。文獻[7]中介紹了在τ分位數下構建Bootstrap樣本的抽樣方法，但是由于本文采用的加權分位數可以同時用于多個分位數下的參數估計，因此本文抽樣誤差需要從τ=0.5時的分位數下產生，具體操作步驟如下。

1)利用容量為N的原始樣本求解分位數參數估計=(,)T。假設共有q個分位數，則記(τq/2)為τ=0.5分位數下的參數估計，利用生成殘差集合{ij}。

最近在市場上多次聽到“兩塊錢，你買不了吃虧；兩塊錢，你買不了上當”這一句式，句中有可能補語標識“買不了(liǎo)”?！俺蕴?、上當”都為動賓結構的謂詞性詞語，正常來說“買”后應該帶名詞性的詞語，這種語言形式從事理上當不符合語法規(guī)則。如何來解釋其存在的原因呢？

6)計算τk分位數下參數的樣本協方差陣

2.2 均值漂移模型的影響分析

2.2.1均值漂移模型

對(i1,j1)個觀測值添加漂移項

(10)

式中γi1j1是一個漂移項。記模型(10)的參數估計為{i1j1}=({i1j1},{i1j1},i1j1)T。此時均值漂移模型在分位數τk(k=1,2，…,q)下的加權損失函數為

(11)

式中γi1j1k是分位數τk下的漂移項。均值漂移模型主要用來對每個數據點下的漂移項逐一進行假設檢驗：H0k∶γi1j1k=0或H1k∶γi1j1k≠0，其中k=1,2，…,q。若拒絕H0k，則說明在τk分位點下該數據點為異常值點。

2.2.2γi1j1k的估計及數據刪除模型和均值漂移模型系數估計的等價性

在分位數τk下，對式(11)關于γi1j1k進行極小化，可得

(12)

(13)

根據式(3)和式(11)可以看出數據刪除模型和均值漂移模型的系數β的估計是等價的，由此可得

(14)

2.2.3診斷統計量

為了得到統計量的抽樣分布，韋博成等[13]在小樣本中針對線性回歸模型的異常值檢驗問題，采用方差分析的方法構造出檢驗統計量；曾林蕊等[14]在誤差的正態(tài)分布假設下，針對大樣本的半參數廣義線性模型，利用懲罰對數似然函數構建了Score檢驗統計量，給出了上述問題的拒絕域。本文在中小樣本和非正態(tài)或非對稱的模型誤差假定下，構造了原假設H0k下的Wald統計量

(15)

式中，Var(γ)是漂移項的方差。統計量(15)在大樣本下近似服從卡方分布，但是在中小樣本下其分布難以確定。因此，本文采取Bootstrap法來構造上述統計量的拒絕域，過程可以分為兩步：首先對統計量中的Var(i1j1k)進行近似計算；其次利用統計量的Bootstrap樣本構造假設檢驗的拒絕域。具體步驟如下所述。

2)利用初始樣本(x,y)求得i1j1k，結合步驟1)中的根據式(15)計算檢驗統計量[15]，記為T0。

3 數值模擬

在模擬中yij代表第i個個體的第j個觀測值，根據如下模型生成數據。

(16)

圖1是0.5分位數下模型的系數估計隨懲罰因子λ的變化曲線，在λ大于7時系數不再發(fā)生變化，說明此時懲罰因子過大，導致隨機效應不再起作用。針對線性混合模型(16)，利用R語言lme4包對其隨機效應方差和誤差方差進行估計，然后根據Koenker[7]建議的方法，并結合圖1，在模擬中經過多次實驗，最終選取λ=2。

圖1 0.5分位數下系數估計與λ的關系Fig.1 The relationship between coefficient estimation and λ in the 0.5 quantile

表1 CDM參數近似估計與真實估計對比Table 1 Comparison between approximate estimationsand true values of CDM parameters

3.1 針對單個點的影響分析

首先對單個數據中是否是強影響點或異常值點進行驗證，實驗可以分為兩部分：針對單個數據是否為強影響點的診斷，和是否為異常值點的診斷。針對由模型(16)所生成的模擬數據，對第30個點和第50個點增加擾動，分別為y30=y30+6和y50=y50-6，然后采用本文提出的基于損失函數的距離和Cook距離來檢測強影響點。強影響點檢驗結果如圖2所示。圖2(a)從上到下分別為分位數等于0.1、0.3、0.5、0.7、0.9時的損失函數距離?？梢钥闯鰮p失函數距離可以很好地檢驗出第30個和第50個點為強影響點，并且在不同分位數下不同位置的強影響點的表現不同：低分位數利于檢驗位于數據下方的強影響點，高分位數則對位于數據上方的強影響點敏感。圖2(b)是Cook距離的診斷結果，盡管結果的波動幅度與圖2(a)不同，但是結論相同，可以相互驗證，避免出現方法不準確的問題。

圖2 單個強影響點的影響分析Fig.2 Influence analysis of a single influential observation

基于相同的模擬數據，對其進行異常值點檢驗，結果展示在圖3中。圖3中虛線上方是Wald統計量的拒絕域，可以明顯看出，圖3(a)中第30個和第50個點位于拒絕域，圖3(b)中第50個點是異常值點，同時結合圖2結果可知，異常值點與強影響點的檢驗結果相同，但是對比圖3(a)和(b)可以看出，低分位數下的結果對位于數據上方的異常值點更敏感。

圖3 單個異常值點的影響分析Fig.3 Influence analysis of a single outlier

3.2 針對同一個體一組點的影響分析

對屬于同一個體的一組數據是否為強影響點或異常值點的情形進行檢驗，類似于3.1節(jié)實驗，分為針對強影響點的診斷和針對異常值點的診斷。對初始模擬數據的第3組數據和第14組數據進行擾動：y3·=y3·+6和y14·=y14·-6。篇幅所限，圖4僅展示一組強影響點的損失函數距離的檢驗結果。從圖4可以看出第3組和第14組數據的損失函數距離明顯大于其他組，因此為強影響點，說明即使是在同一個體的一組數據都是強影響點的情況下，本文所提方法也具有良好的檢驗能力。

圖4 一組強影響點的損失函數距離Fig.4 Distance of loss function for a group of influence observations

4 實驗驗證與分析

實驗數據共包含216個觀測值，其中輻解劑量有1、2、4、6、8、10(103Gy)共6個水平，硝酸濃度也有0.2、0.4、0.6、0.8、1.0、1.5 mol/L共6個水平，HSC濃度為0.1、0.2、0.4、0.6、0.8、1.0 mol/L共6個水平。根據1.1節(jié)數據描述的結果，以輻解劑量作為解釋變量Xij，硝酸濃度和HSC濃度的聯合影響作為隨機效應μi，共有36個取值。響應變量cij表示為不同實驗條件下生成的銨根離子濃度(mmol/L)。為了對響應變量進行影響因素分析，可建立如下線性混合效應模型

cij=β0+β1Xij+μi+εij

(17)

式中，i=1,…,36，j=1,…,6。基于模型(17)在0.5分位數下的參數估計作誤差散點圖，如圖5所示，可以看出誤差在零點上下波動，但是在上側波動小而密集，在下側波動大而疏松，可以明顯得出數據總體上既不是對稱分布，也不是正態(tài)分布，因此選擇采用本文方法對其進行影響分析。

圖5 0.5分位數下的誤差散點圖Fig.5 Scatter plot of error in the 0.5 quantile

由于在低分位數下檢驗結果不存在明顯的強影響點，所以圖6僅展示了τ=0.9時的檢驗結果。從圖6可以看出，在0.9的分位數下，第125個數據點和第35組數據相較于其他數據明顯突出，結合圖5可知二者均是位于數據上方的強影響點。

圖6 實證數據影響分析結果(τ=0.9)Fig.6 Influence analysis results for real data (τ=0.9)

5 結束語

本文針對服從非對稱、非正態(tài)，尤其是重尾或輕尾分布的縱向數據的影響分析問題，在中等樣本下利用線性混合模型的分位數估計，分別構建檢測強影響點的基于損失函數的距離和Cook距離，以及識別異常值點的Wald統計量。此外，為了減少計算量，給出數據刪除模型參數的一步近似估計，并且利用Bootstrap方法求得統計量的拒絕域。模擬結果表明，本文所提的診斷度量和診斷統計量都可以很好地檢驗出數據中的強影響點和異常值點，實證分析結果也進一步證明了本文方法的實用性。