纵向数据下广义经验似然方法的有效稳健估计

（上海理工大学 理学院 上海 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