李宏凱 肖松濤 歐陽應(yīng)根 李志強(qiáng)*

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

引 言

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

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

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

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

1 實(shí)驗(yàn)數(shù)據(jù)描述及線性混合模型

1.1 數(shù)據(jù)描述

為了對(duì)氨基羥基脲(HSC)水溶液及硝酸水溶液在γ射線作用下的輻解產(chǎn)物進(jìn)定量分析，考察不同輻解劑量、HSC初始濃度以及硝酸濃度條件下主要輻解產(chǎn)物的生成量，分別配制濃度為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種不同濃度配比的溶液，委托中國(guó)原子能科學(xué)研究院輻照中心，采用實(shí)驗(yàn)型鈷源裝置(60Co源裝置)，分別置于輻解劑量為1×103、2×103、4×103、6×103、8×103、10×103Gy的條件下，進(jìn)行三因素下的無重復(fù)實(shí)驗(yàn)，總共得到216個(gè)樣品。

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

1.2 線性混合模型及其分位數(shù)估計(jì)

(1)

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

(2)

為了解決估計(jì)隨機(jī)效應(yīng)時(shí)可能會(huì)面臨的維數(shù)較高的問題，Koenker[7]提出一種基于隨機(jī)效應(yīng)的L1懲罰方法，通過加權(quán)多個(gè)分位數(shù)下的信息來獲取參數(shù)與隨機(jī)效應(yīng)的聯(lián)合估計(jì)，其損失函數(shù)定義如下

(3)

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

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

2 線性混合效應(yīng)模型分位數(shù)估計(jì)的影響分析

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

2.1 數(shù)據(jù)刪除模型的影響分析

2.1.1數(shù)據(jù)刪除模型

數(shù)據(jù)刪除模型主要通過依次刪除每個(gè)數(shù)據(jù)點(diǎn)來求解參數(shù)，之后再對(duì)比各個(gè)參數(shù)估計(jì)，若數(shù)據(jù)中存在較強(qiáng)的影響點(diǎn)，那么當(dāng)刪除該點(diǎn)時(shí)，模型的參數(shù)估計(jì)相對(duì)于原估計(jì)就會(huì)有較大的變化。

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

(4)

2.1.2參數(shù)一步(光滑)近似估計(jì)

模型(4)相較原模型(1)只是減少了1個(gè)數(shù)據(jù)，所以參數(shù)估計(jì)的方法及計(jì)算復(fù)雜度與原模型基本相同，但是當(dāng)數(shù)據(jù)量很大時(shí)，由于需要逐一檢驗(yàn)每個(gè)數(shù)據(jù)點(diǎn)，將會(huì)產(chǎn)生很大的計(jì)算負(fù)擔(dān)。針對(duì)線性模型，Cook等[10]提出了利用似然函數(shù)的Taylor展開式構(gòu)造參數(shù)近似估計(jì)的思想。

為了能夠?qū)Σ豢蓪?dǎo)的損失函數(shù)(3)進(jìn)行光滑逼近，以建立數(shù)據(jù)刪除模型參數(shù)的一步近似估計(jì)，本文結(jié)合一般分位數(shù)參數(shù)估計(jì)的minorize- maximization(MM)算法[11]和絕對(duì)值函數(shù)的光滑方法[12]，構(gòu)造如下目標(biāo)函數(shù)

(5)

(6)

對(duì)式(6)兩端關(guān)于η[i1j1]求導(dǎo)可得[10]

(7)

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

1)根據(jù)損失函數(shù)(3)，構(gòu)建基于損失函數(shù)的距離

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

(8)

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

2)Cook距離定義為[10]

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

(9)

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

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

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

6)計(jì)算τk分位數(shù)下參數(shù)的樣本協(xié)方差陣

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

2.2.1均值漂移模型

對(duì)(i1,j1)個(gè)觀測(cè)值添加漂移項(xiàng)

(10)

式中γi1j1是一個(gè)漂移項(xiàng)。記模型(10)的參數(shù)估計(jì)為{i1j1}=({i1j1},{i1j1},i1j1)T。此時(shí)均值漂移模型在分位數(shù)τk(k=1,2，…,q)下的加權(quán)損失函數(shù)為

(11)

式中γi1j1k是分位數(shù)τk下的漂移項(xiàng)。均值漂移模型主要用來對(duì)每個(gè)數(shù)據(jù)點(diǎn)下的漂移項(xiàng)逐一進(jìn)行假設(shè)檢驗(yàn)：H0k∶γi1j1k=0或H1k∶γi1j1k≠0，其中k=1,2，…,q。若拒絕H0k，則說明在τk分位點(diǎn)下該數(shù)據(jù)點(diǎn)為異常值點(diǎn)。

2.2.2γi1j1k的估計(jì)及數(shù)據(jù)刪除模型和均值漂移模型系數(shù)估計(jì)的等價(jià)性

在分位數(shù)τk下，對(duì)式(11)關(guān)于γi1j1k進(jìn)行極小化，可得

(12)

(13)

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

(14)

2.2.3診斷統(tǒng)計(jì)量

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

(15)

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

2)利用初始樣本(x,y)求得i1j1k，結(jié)合步驟1)中的根據(jù)式(15)計(jì)算檢驗(yàn)統(tǒng)計(jì)量[15]，記為T0。

3 數(shù)值模擬

在模擬中yij代表第i個(gè)個(gè)體的第j個(gè)觀測(cè)值，根據(jù)如下模型生成數(shù)據(jù)。

(16)

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

圖1 0.5分位數(shù)下系數(shù)估計(jì)與λ的關(guān)系Fig.1 The relationship between coefficient estimation and λ in the 0.5 quantile

表1 CDM參數(shù)近似估計(jì)與真實(shí)估計(jì)對(duì)比Table 1 Comparison between approximate estimationsand true values of CDM parameters

3.1 針對(duì)單個(gè)點(diǎn)的影響分析

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

圖2 單個(gè)強(qiáng)影響點(diǎn)的影響分析Fig.2 Influence analysis of a single influential observation

基于相同的模擬數(shù)據(jù)，對(duì)其進(jìn)行異常值點(diǎn)檢驗(yàn)，結(jié)果展示在圖3中。圖3中虛線上方是Wald統(tǒng)計(jì)量的拒絕域，可以明顯看出，圖3(a)中第30個(gè)和第50個(gè)點(diǎn)位于拒絕域，圖3(b)中第50個(gè)點(diǎn)是異常值點(diǎn)，同時(shí)結(jié)合圖2結(jié)果可知，異常值點(diǎn)與強(qiáng)影響點(diǎn)的檢驗(yàn)結(jié)果相同，但是對(duì)比圖3(a)和(b)可以看出，低分位數(shù)下的結(jié)果對(duì)位于數(shù)據(jù)上方的異常值點(diǎn)更敏感。

圖3 單個(gè)異常值點(diǎn)的影響分析Fig.3 Influence analysis of a single outlier

3.2 針對(duì)同一個(gè)體一組點(diǎn)的影響分析

對(duì)屬于同一個(gè)體的一組數(shù)據(jù)是否為強(qiáng)影響點(diǎn)或異常值點(diǎn)的情形進(jìn)行檢驗(yàn)，類似于3.1節(jié)實(shí)驗(yàn)，分為針對(duì)強(qiáng)影響點(diǎn)的診斷和針對(duì)異常值點(diǎn)的診斷。對(duì)初始模擬數(shù)據(jù)的第3組數(shù)據(jù)和第14組數(shù)據(jù)進(jìn)行擾動(dòng)：y3·=y3·+6和y14·=y14·-6。篇幅所限，圖4僅展示一組強(qiáng)影響點(diǎn)的損失函數(shù)距離的檢驗(yàn)結(jié)果。從圖4可以看出第3組和第14組數(shù)據(jù)的損失函數(shù)距離明顯大于其他組，因此為強(qiáng)影響點(diǎn)，說明即使是在同一個(gè)體的一組數(shù)據(jù)都是強(qiáng)影響點(diǎn)的情況下，本文所提方法也具有良好的檢驗(yàn)?zāi)芰Α?/p>

圖4 一組強(qiáng)影響點(diǎn)的損失函數(shù)距離Fig.4 Distance of loss function for a group of influence observations

4 實(shí)驗(yàn)驗(yàn)證與分析

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

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

(17)

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

圖5 0.5分位數(shù)下的誤差散點(diǎn)圖Fig.5 Scatter plot of error in the 0.5 quantile

由于在低分位數(shù)下檢驗(yàn)結(jié)果不存在明顯的強(qiáng)影響點(diǎn)，所以圖6僅展示了τ=0.9時(shí)的檢驗(yàn)結(jié)果。從圖6可以看出，在0.9的分位數(shù)下，第125個(gè)數(shù)據(jù)點(diǎn)和第35組數(shù)據(jù)相較于其他數(shù)據(jù)明顯突出，結(jié)合圖5可知二者均是位于數(shù)據(jù)上方的強(qiáng)影響點(diǎn)。

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

5 結(jié)束語

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