縱向數據下廣義經驗似然方法的有效穩健估計

（上海理工大學 理學院 上海 200093）

縱向數據是指對同一個個體或受試單位在不同時間進行重復觀測得到的數據。觀測不同個體所得數據是獨立的，但對同一個個體觀測所得的數據往往具有相關性。因此，縱向數據具有組間獨立、組內相關的特點，是數理統計研究中的復雜數據類型之一，在生物、醫學及經濟學等領域都有廣泛的應用。在實際中經常需要對縱向數據進行統計分析和建模。

對于縱向數據的研究，學者們已經提出了各種參數模型和統計方法。參數模型易于解釋和統計分析，能很好地反映協變量與響應變量的關系，但存在模型誤判風險，一旦模型假設錯誤，直接導致得到錯誤的結論；非參數方法不需要對模型進行假定，具有非常高的靈活性，廣義經驗似然估計是非參數方法中常用的一種。Thomas等[1]最早闡述了經驗似然方法的思想；Owen[2]進一步系統地研究了經驗似然方法，他指出經驗似然是一種構造未知參數的置信區間和置信域的非參數推斷方法，其本質是在約束條件下求非參數似然比的極大值，而總體參數由約束條件帶入似然比中。經驗似然方法自提出后即引起了很多統計學家的興趣，他們將這一方法應用到統計的諸多領域；1991年Owen[3]將經驗似然方法應用于獨立同分布樣本總體均值的統計推斷中，構造了非參數經驗似然比統計量，并進一步對線性模型和廣義線性模型進行了深入的研究；Bondell等[4]提出了一種基于最小化C-R距離的兩階段廣義經驗似然方法,該方法對殘差大的樣本賦予小的權重，從而達到對異常值穩健的目的。

經驗似然方法與一些經典的統計方法相比，其優點在于它不需要對模型作任何參數形式的假設，也不需要估計方差，只要矩條件正確，則估計一定是漸近有效的。

本文研究的內容是縱向數據的均值和協方差的聯合估計，對此前人已作了許多研究，He等[5]提出了一種基于半參數廣義線性模型的穩健估計方程，用來估計縱向數據的均值和協方差，通過樣條回歸和得分函數從而達到對異常值穩健的目的；Qin等[6]基于廣義估計方程，提出了逆概率加權法來估計縱向數據的均值和協方差，估計效果較穩健；Qin等[7]建立了穩健估計方程，利用Cholesky分解，實現了縱向數據在缺失機制下，均值和協方差的聯合穩健估計，并且對數據異常值具有一定的抵抗性。Bondell等[4]通過廣義經驗似然方法中約束條件的限制，在保證穩健性的同時又提高了估計的有效性。但Bondell等[4]的方法只適用于橫截面數據，沒有考慮縱向數據。本文在Bondell等[4]的啟發下，結合Cholesky分解將縱向數據的線性模型重參數化，利用廣義經驗似然估計方法,實現了縱向數據的均值和協方差的穩健聯合估計，在保證估計的穩健性的同時，還提高了估計的有效性。

1 模型與估計方法

1.1 初始模型

1.2 模型重參數化

1.3 廣義經驗似然估計

2 交替優化算法

為了解決帶有約束條件的目標優化函數式（4），采用拉格朗日算法。

3 模擬研究

β0=(β1,β2,β3)T=(1,1,1)T，模擬重復次數為100次。 θ的初始估計值選取數據在最小二乘估計下的參數估計值OLS，記Xu為本文提出的方法，He是He等[5]提出的穩健廣義估計方程的方法，Qin是Qin等[7]提出的穩健方程的估計方法，由于Qin等[7]研究的縱向數據帶有缺失，而本文只考慮了污染，故在模擬時將缺失示性函數全部改為1，代表數據不缺失。此外，數據可能存在污染等異常情況，本文β0是最終估計式（7）中的前p個分量，協方差 Σ通過返還函數對反向分解和還原可得，在返還函數中，R表示采用了絕對中位數離差MAD方法，NR表示一般處理方法，以此來驗證本文提出的方法是否對異常值具備抵抗性。另外，兩種方法中R和NR均代表各自的穩健方法和非穩健方法，表1比較了3種方法得到的0的偏差(BIAS)、標準差(SE)和均方誤差(MSE)。

表1 數據不污染時參數估計對比Tab.1 Parameter estimation and comparison when data are not polluted

為了進一步說明本文提出方法的穩健性，對數據進行3%污染，污染方式為隨機挑選樣本3%的xi將其替換成xi-2，并將xi所對應的yi替換成yi+2，所得β0估計結果如表2所示。

通過比較可以看出：在數據不污染和3%污染的情況下，本文提出的穩健方法均要優于其他2種方法，參數的偏差和均方誤差都比其他2種方法要小，顯示估計更準確；當數據不污染時，本文提出的方法與其他方法比較，主要是參數的偏差比較小，從而使得本文方法的有效性比其他2種方法高；當數據被污染時，本文提出的方法的有效性更為顯著。相同條件下，本文提出方法的偏差遠遠小于其他2種方法，以參數β3為例，固定協方差為I單位結構，在數據不污染時，相同穩健的方法，本文β3偏差分別比Qin和He小0.0002，0.006；在數據3%污染時，相同穩健的方法，本文β3偏差分別比Qin和He小0.001，0.0111。

此外，本文還對數據進行了5%的污染，結果顯示，數據污染程度越大，本文的方法越有效，為了節省篇幅，相關表格在本文沒有給出。

為了進一步說明對縱向數據協方差矩陣的估計好壞程度，本文定義指標QL和EL來衡量與Σ的接近程度，QL和EL越小，表明協方差矩陣的估計值越接近真實值。

表2 數據 3% 污染時參數估計對比Tab.2 Parameter estimation and comparison when data are 3% polluted

在數據不污染和3%污染的情況下，將本文方法和Qin等[7]提出的穩健方程進行比較，對應的QL和EL結果如表3所示，其中，ELi,QLi(i=1,2,3,4)代表誤差項服從第i種協方差結構下的QL和EL指標，對應模型中設定的4種協方差結構。

由表3可以看出：在數據不污染時，本文提出的非穩健方法所估計出的協方差的EL指標，一直優于Qin提出的方法；在數據樣本受到3%污染時，本文提出的穩健方法所估計的協方差的EL和QL指標，都幾乎一直優于Qin提出的方法，表明本文方法估計出的協方差結構與真實值最為接近，說明本文提出的方法具有顯著穩健性和有效性，不僅能保證參數β0和 Σ估計的有效性，而且在樣本受到污染時仍能保證其準確性；此處還對數據進行了5%的污染，得到結果與上面類似，并且污染程度越大，本文方法的優越性越顯著。

為了進一步比較經驗似然方法的優越性，選取誤差項服從多元t分布的縱向數據進行模擬研究，研究方法同上，所得結果如表4所示。

由表4可知，在誤差服從t分布時，經驗似然方法的估計結果仍然非常可觀。3種方法的穩健結果顯示，在4種協方差結構下，本文方法的估計方差要遠遠小于其他2種方法，說明本文提出的經驗似然估計方法比其他2種方法更有效，參數估計更趨于一致。此外，對比誤差項服從多元t分布下本文方法和Qin等[7]的4種協方差估計情況如表5所示。

由表5可知，在誤差服從多元t分布時，經驗似然方法的有效性更為顯著。根據QL指標和EL指標顯示，本文估計的縱向數據協方差矩陣與真實協方差矩陣相似度更高，與表3對比可看出，當模擬數據的誤差服從多元t分布時，其協方差矩陣的估計效果更優于誤差服從多元正態分布時的情形，進一步說明了本文方法對誤差具有抵抗性。

表3 EL和 QL對比表Tab.3 EL and QL comparison table

表4 t分布誤差下的參數估計對比表Tab.4 Comparison of parameter estimates under t distribution error

表5 t分布誤差下 QL和 ELTab.5 EL and QL under t distribution error

4 實證分析

為了進一步驗證本文方法的實用性，分析了CD4細胞計數縱向數據。關于這個數據集的完整描述見 Diggle的主頁:http://www.lancs.ac.uk/diggle/.這里分析所用的數據是完整數據的截取部分,包含340個樣本的1020次觀測值。響應變量y是CD4計數的算數平方根,協變量包括血清轉換的時間x2，相對于一個起點的年齡x3，由吸煙的包數刻畫的吸煙狀況x4，娛樂性藥物使用是/否x5，性伴侶的個數x6，以及流行病中心給出抑郁狀態和抑郁程度x7。有許多學者研究過這個數據集。Zeger等[10]和Wang等[11]分別對這個數據擬合了半參數模型和部分線性半參數模型,其中，協變量x2作為非線性形式進入模型。

本文分析的主要目的是尋找CD4數量和協變量之間的關系，考慮到可能存在非線性關系，引入了除x5以外的平方項，加上截距項x1，記

由于實際數據無法知道真實參數，所以，現利用交互驗證的方法來比較各種方法的優劣，交互驗證過程中的均方誤差

這里的n=339，(-i)是去掉第i個個體外，由其余339個個體擬合所得的預測值，數據預處理后，代入算法中，本文方法與Qin和He方法的比較結果如表6所示。

表6說明，本文方法在實際應用中十分有效，CV值要明顯小于另外2種方法的CV值，值得推廣。

表6 CV值比較Tab.6 Comparison of CV value