上海理工大学学报  2019, Vol. 41 Issue (4): 331-338   PDF    
纵向数据下广义经验似然方法的有效稳健估计
徐孝琳, 樊亚莉, 苏依官     
上海理工大学 理学院 上海 200093
摘要: 基于广义经验似然估计方法,提出了一种有效且稳健的估计,实现对纵向数据在线性模型下均值和协方差矩阵的联合估计。利用Cholysky分解将模型重参数化,利用拉格朗日乘子法求出估计值,再还原出均值和协方差矩阵的估计。在模拟研究中将所提方法同文献中其他稳健估计进行比较,结果显示所提方法效率更高。最后将所提方法用于分析CD4细胞数据,交叉验证结果显示所提方法更加可靠。
关键词: 纵向数据     稳健性     有效性     经验似然     交替优化    
Effective Robust Estimation Based on the Generalized Empirical Likelihood Method for Longitudinal Data
XU Xiaolin, FAN Yali, SU Yiguan     
College of Science, University of Shanghai for Science and Technology, Shanghai 200093, China
Abstract: Based on the generalized empirical likelihood estimation method, a robust and effective estimation was proposed to realize the joint estimation of the mean and covariance matrix of longitudinal data in linear models. Using the Cholysky decomposition, the model was re-parameterized, and using the Lagrangian multiplier method, estimates were obtained, then, an estimate of the mean and covariance matrix was restored. Comparing the proposed method with other robust estimates in the literature in a simulation study, the results show that the proposed method is more efficient. Finally, the proposed method was used to analyze CD4 cell data. Cross-validation results show that the proposed method is more reliable.
Key words: longitudinal data     robustness     effectivenes     empirical likelihood     alternating optimization    

纵向数据是指对同一个个体或受试单位在不同时间进行重复观测得到的数据。观测不同个体所得数据是独立的,但对同一个个体观测所得的数据往往具有相关性。因此,纵向数据具有组间独立、组内相关的特点,是数理统计研究中的复杂数据类型之一,在生物、医学及经济学等领域都有广泛的应用。在实际中经常需要对纵向数据进行统计分析和建模。

对于纵向数据的研究,学者们已经提出了各种参数模型和统计方法。参数模型易于解释和统计分析,能很好地反映协变量与响应变量的关系,但存在模型误判风险,一旦模型假设错误,直接导致得到错误的结论;非参数方法不需要对模型进行假定,具有非常高的灵活性,广义经验似然估计是非参数方法中常用的一种。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 初始模型

假设纵向数据有 $n$ 个对象,对于每一个对象有 $m$ 次观察值,考虑线性模型

${y_{ij}} = {{x}}_{ij}^{\rm T}{{{\beta}} _0} + {\varepsilon _{ij}},\;\;\;\;\;\;i = 1,2, \cdots, n,\;\;\;j = 1,2, \cdots, m$ (1)

式中: $\,{{{\beta }}_0}$ $p$ 维未知的回归系数向量; ${{{\varepsilon }}_i} = $ $ {({\varepsilon _{i1}},{\varepsilon _{i2}}, \cdots ,{\varepsilon _{im}})^{\rm T}}\left( {i = 1,2, \cdots, n} \right)$ 独立同分布,且服从于均值为0m×1,协方差 ${{{\varSigma }}_0}$ 未知的分布。

${{{Y}}_i} = {({y_{i1}}, \cdots, {y_{im}})^{\rm T}}$ ${{{X}}_i} = {({{{x}}_{i1}}, \cdots ,{{{x}}_{im}})^{\rm T}}$ ,在本文中主要研究模型(1)中 $\,{{{\beta }}_0}$ ${{{\varSigma }}_0}$ 的有效稳健估计。

1.2 模型重参数化

首先对 ${{{\varSigma }}_0}$ 进行Cholesky分解,实现 $\;{{{\beta }}_0}$ ${{{\varSigma }}_0}$ 的联合估计。即存在对角线元素都是1的下三角矩阵 ${{C}}$ ,使得 ${\rm{cov}} ({{C}}{{{\varepsilon }}_i}) = {{C}}{{{\varSigma }}_0}{{{C}}^{\rm T}} = {{B}}$ ,其中, ${{B}}$ 是对角矩阵,记 ${{B}} = {\rm {diag}}\;(b_1^2,b_2^2, \cdots ,b_m^2)$ 。令

${{C}} = \left( {\begin{array}{*{20}{c}} 1 & 0 & \cdots & 0\\ { - {c_{21}}}& 1& \cdots & 0\\ \vdots & \vdots & {}& \vdots \\ { - {c_{m1}}}& { - {c_{m2}}}& \cdots &1 \end{array}} \right)\;\;$

${{{e}}_i} \!\!=\!\! {({e_{i1}}, \!\cdots\!, {e_{im}})^{\rm T}}\! \!=\!\! {{C}}{{{\varepsilon }}_i}$ ${\rm{var}} ({e_{ij}}) \!\!=\!\! b_j^2\;,j \!\!=\!\! 1,2, \!\cdots\! ,m$

则有 ${\varepsilon _{i1}} = {e_{i1}}$ ${\varepsilon _{ij}} \!=\! {c_{j1}}{\varepsilon _{i1}} \!+\! {c_{j2}}{\varepsilon _{i2}} \!+\! \cdots +{c_{j,j - 1}}{\varepsilon _{i,j - 1}} + {e_{ij}}\;,$ $i = 1,2, \cdots n,j = 1,2, \cdots, m$

根据上述关系,模型(1)可以改写为

$\left\{ \begin{array}{l} {y_{i1}} = {{x}}_{i1}^{\rm T}{{{\beta }}_0} + {e_{i1}} \\ {y_{ij}} = {{x}}_{ij}^{\rm T}{{{\beta }}_0} + {c_{j1}}{\varepsilon _{i1}} + \cdots+ {c_{jj - 1}}{\varepsilon _{ij - 1}} + {e_{ij}},\\ \;\;\; i = 1,2, \cdots, n,\quad j = 2,3 \cdots, m \\ \end{array} \right.$ (2)

此处由于 ${\varepsilon _{ij}}$ 未知,故先对模型 $y{}_{ij} = {{x}}_{ij}^{\rm T}{{{\beta }}_0} + {{{\varepsilon }}_{ij}},\;i = $ $ 1,2, \cdots ,n,\;\;j = 1,2, \cdots ,m$ 中的数据用最小二乘估计出初始 $\;{{\hat{ \beta }}_0}$ ,则 ${\hat \varepsilon _{ij}} = {y_{ij}} - {{x}}_{ij}^{\rm T}{{\hat{ \beta }}_0}$ ,令 ${\varepsilon _{ij}} = {\hat \varepsilon _{ij}}$ ,代入模型(2)中。

当得到 ${c_{ij}}(i = 1,2 \cdots, m,j = 1,2\cdots ,\left( {m - 1} \right) )$ ${b_{j}}(j = $ $1,2 ,\cdots, m)$ $\;{{{\beta }}_0}$ 的相合估计值后,则可根据Cholesky分解得 ${{\hat{ \varSigma }}_0} = {{\hat{ C}}^{ - 1}}{\hat{ B}}{({{\hat{ C}}^{\rm T}})^{ - 1}}$ ,从而实现均值和协方差矩阵的联合估计。

基于Cholesky分解,为了将模型(1)重参数化后改写成一般的模型,需作如下处理:

${{{s}}_0} \!=\! {({c_{21}},{c_{31}},{c_{32}}, \cdots, {c_{m,m \!-\! 1}})^{\rm T}}$ ${{\hat{ Z}}_{ij}} = ({{\textit{0}}}_{(j - 2)(j - 1)/2}^{\rm T}, $ ${\hat \varepsilon _{i1}}, \cdots ,{\hat \varepsilon _{i,j - 1}},{{\textit{0}}}_{(m - 1)m/2 - (j - 1)j/2}^{\rm T})^{\rm T}$ ${{\hat{ Z}}_i}^{( \!-\! 1)} \!\!=\!\! {({{\hat{ Z}}_{i2}}, \cdots ,{{\hat{ Z}}_{im}})^{\rm T}}$ ${{{\theta }}_0} \!\!=\!\! {({{\beta }}_0^{\rm T},{{s}}_0^{\rm T})^{\rm T}}$ ${{X}}_i^{( \!-\! 1)}\! \!=\!\! ({{{x}}_{i2}}, \cdots ,$ ${{{x}}_{im}})^{\rm T} $ ${{{D}}_1} = ({{x}}_{i1}^{\rm T},{{\textit{0}}}_{m(m - 1)/2}^{\rm T})$

${\hat{ D}}_i^{( - 1)} = ({{X}}_i^{( - 1)},{\hat{ Z}}_i^{( - 1)})$ ${\hat{ D}}_i^{} = {({{{D}}_1},{\hat{ D}}_i^{( - 1)})^{\rm T}}$

在此变换后,原始的纵向数据的模型(1)被改写成 ${{{Y}}_i} = {{\hat{ D}}_i}{{\theta }} + {{{e}}_i}\;,\;\;\;i = 1,2, \cdots, n$ ,其中, ${{{Y}}_i}$ ${{{e}}_i}\;$ 都是 $m$ 维列向量,进一步记

   ${{Y}} ={({{{Y}}_1},{{{Y}}_2}, \cdots, {{{Y}}_n})^{\rm T}} = {({y_1},{y_2}, \cdots, {y_{nm}})^{\rm T}}\qquad$      ${\hat{ D}} = {({{\hat{ D}}_1},{{\hat{ D}}_2}, \cdots, {{\hat{ D}}_n})^{\rm T}} = {({{d}}_1^T,{{d}}_2^{\rm T}, \cdots, {{d}}_{nm}^{\rm T})^{\rm T}}$      ${{e}} = {({{{e}}_1},{{{e}}_2}, \cdots ,{{{e}}_n})^{\rm T}} = {({\omega _1},{\omega _2}, \cdots, {\omega _{nm}})^{\rm T}}\quad$

则模型最终化为

${y_k} = {{d}}_k^{\rm T}{{\theta }} + {\omega _k},\quad k = 1,2, \cdots ,N$ (3)

式中: $N = nm$ ${y_k}$ 是一个数; ${\omega _k}$ $(k=1,2,\cdots,N) $ 相互独立; ${{{d}}_k}$ $[p + (m - 1)m/2]$ 维列向量; ${{\theta }}$ $[p + (m - 1)m/2]$ 维列向量。

1.3 广义经验似然估计

定义广义经验似然目标函数

$L({{\theta }},{{w}}) = \min \left\{ {\sum\limits_{k = 1}^N {{w_k}\log (N{w_k})} } \right\}$ (4)

约束条件

$ \sum\limits_{k = 1}^N {{w_k}} = 1,\;\;\;\sum\limits_{k = 1}^N {{w_k}} ({y_k} - {{d}}_k^{\rm T}{{\theta }}){{{d}}_k} = 0\;\; $ (5)
$\sum\limits_{k = 1}^N {{w_k}} {({y_k} - {{d}}_{k}^{\rm T}{{\theta }})^2} = \tilde \sigma _T^2\;\;\;\;\;$ (6)

在上述条件中, ${w_k}$ 代表权重, $N = n m$ ${{w}} = \left({w_1},\right. $ $\left.{w_2}, \cdots ,{w_N}\right)$ 。本文提出的估计量是一个广义经验似然估计量,它是经验似然方法的延伸,是基于最小化C-R距离的统计量,C-R距离最小化也就是意味着其经验分布差异最小化。当权重无限接近于等权重时,才能使得估计的有效性达到最好,具体可参考文献[8-9]。

本文所提出的估计是两阶段的估计,第一阶段,在约束条件式(5)和式(6)下,目标函数式(4)达到最小值时,可得对应 ${{\theta }}$ ${w_k}$ 的估计值,即

$ {{\tilde \theta }} = \arg \min\; L({{\theta }},{{w}}), \quad{\hat w_k} = \arg \min \;L({{\theta }},{{w}})。$

第二阶段, ${{\theta }}$ 的最终估计值是基于模型(3)和 ${\hat w_k}$ 的加权估计量,记 ${{W}} = {\rm{diag}}({\hat w_1},{\hat w_2} ,\cdots, {\hat w_N})$ ${{D}} = ({{{d}}_1},{{{d}}_2}, \cdots, {{{d}}_N})$ ${{Y}} = {({y_1},{y_2}, \cdots, {y_N})^{\rm T}}$ 。则 ${{\theta }}$ 的最终估计值为

${\hat{ \theta }} = ({{{D}}^{\rm T}}{{WD}}){}^{ - 1}{{{D}}^{\rm T}}{{WY}}$ (7)

在上述条件中, ${w_k}$ 代表权重,在仅考虑条件式(5)时,目标函数将在 ${w_k} = 1/N$ $k = 1,2, \cdots ,N$ 达到全局最小。约束条件式(5)保证了该经验似然方法的有效性,加入条件式(6)则进一步保证本文提出方法的稳健性,其中, $\;\tilde \sigma _T^2\;$ 是一个需要预估的残差尺度参数,此处需要求 $\;\tilde \sigma _T^2\; \leqslant \hat \sigma _{\rm{OLS}}^2$ $\hat \sigma _{\rm{OLS}}^2 = $ $(1/N)\displaystyle\sum\limits_{k = 1}^N {({y_k}} - {{d}}_k^{\rm T}{{\hat{ \theta }}_{\rm {OLS}}}{)^2}$ $\;\tilde \sigma _T^2\;$ 的引入,使得本文提出的方法不仅有效,而且稳健。现进一步说明估计量式(7)稳健的原因。

$r_k^2({\hat{ \theta }}) = {y_k} - {{d}}_k^{\rm T}{\hat{ \theta }}$ ,因为, $\displaystyle\sum\limits_{k = 1}^N {{{\hat w}_k}\left\{ {r_k^2({{\theta }}) - \tilde \sigma _T^2} \right\} = } 0$ ,且 $\displaystyle\sum\limits_{k = 1}^N {\hat w{}_k} = 1$ 。在假设 $\;\tilde \sigma _T^2\; \leqslant \;\hat \sigma _{{\rm{OLS}}}^2$ 的前提下,可得

$\begin{split}& \quad \sum\limits_{k = 1}^N {{{\hat w}_k}} r_k^2({\hat{ \theta }}) =\sum\limits_{k = 1}^N {{{\hat w}_k}} \tilde \sigma _T^2 = \tilde \sigma _T^2 \leqslant \hat \sigma _{\rm{OLS}}^2= \\ &\qquad(1/N)\sum\limits_{k = 1}^N {r_k^2({{{\hat{ \theta }}}_{\rm{OLS}}})} \leqslant (1/N)\sum\limits_{k = 1}^N {r_k^2({\hat{ \theta }})}\qquad\qquad\qquad\qquad\\[-20pt]\end{split} $ (8)

因为, $\;(1/N)\displaystyle\sum\limits_{k = 1}^N {{w_k}} = 1/N\;$ ,所以,式(8)可变为

$\;\sum\limits_{k = 1}^N {{{\hat w}_k}r_k^2({\hat{ \theta }})} \leqslant (1/N)\sum\limits_{k = 1}^N {r_k^2({\hat{ \theta }})} \Rightarrow \sum\limits_{k = 1}^N {\left\{ {{{\hat w}_k} - (1/N)} \right\}r_k^2({\hat{ \theta }})} \leqslant 0$

因此, $r_k^2({\hat{ \theta }})$ ${\hat w_k}$ 是负相关的。当残差 $r_k^2({\hat{ \theta }})$ 比较大时, ${\hat w_k}$ 比较小;当残差 $r_k^2({\hat{ \theta }})$ 比较小时, ${\hat w_k}$ 比较大。这也是本文方法稳健的根本原因。

总结本文估计方法的步骤如下:

a.根据每一组对象的第一个观测值,利用最小二乘方法,得到 $\,{{{\beta }}_0}$ 的初始估计值,记为 $\,{{\hat{ \beta }}_b}$ ,从而有 ${\hat \varepsilon _{ij}} = {y_{ij}} - {{x}}_{ij}^{\rm T}{{\hat{ \beta }}_b}$ ,并且令 ${\varepsilon _{ij}} = {\hat \varepsilon _{ij}}$

b.对于模型(3), ${{\theta }} = {({{\beta }}_0^{\rm T},{{s}}_0^{\rm T})^{\rm T}}$ ,采用Bondell等[4]提出的两阶段方法,得到 ${\hat{ \theta }}$ ,提取 ${\hat{ \theta }}$ 的前 $p$ 个分量即为 $\,{{\hat{ \beta }}_0}$ ,再根据 $\,{\hat \omega _k} = y{}_k - {{d}}_k^{\rm T}{\hat{ \theta }}$ ,利用绝对中位数离差MAD方法得到 $b_j^2$ 的稳健估计。

c.当求得 ${{\hat{ s}}_0}$ $\hat b_j^2$ 后,可根据Cholesky分解还原,就能得到纵向数据协方差矩阵 ${{\hat{ \varSigma }}_0}$

2 交替优化算法

为了解决带有约束条件的目标优化函数式(4),采用拉格朗日算法。

定义 ${{\varLambda }} = {({\lambda _1},{{\lambda }}_2^{\rm T},{\lambda _3})^{\rm T}}$

$\begin{split}&\quad {L_0}({{w}},{{\theta }},{{\varLambda }}) =\displaystyle \sum\limits_{k = 1}^N {{w_k}\log (N{w_k})} - \\ &\qquad {\lambda _1}\left( {\displaystyle\sum\limits_{k = 1}^N {{w_k} - 1} } \right) -{{\lambda }}_2^{\rm T}\displaystyle\sum\limits_{k = 1}^N {{w_k}\left( {{y_k} - {{d}}_k^{\rm T}{{\theta }}} \right){{{d}}_k}} -\qquad\qquad\qquad\qquad\qquad \\ &\qquad {\lambda _3}\displaystyle\sum\limits_{k = 1}^N {{w_k}\left\{ {{{\left( {{y_k} - {{d}}_k^{\rm T}{{\theta }}} \right)}^2} - \tilde \sigma _{T}^2} \right\}} \end{split} $

$\partial {L_0}/\partial {w_k} = 0,\; k = 1,2, \cdots, N$ $\partial {L_0}/\partial {\lambda _1} \!= \!0$ ,解得

${w_k} = w_k^*\Biggr/\displaystyle\sum\limits_{k = 1}^N {w_k^*}$

其中,

   $w_k^ * = \exp\left[ {{\lambda }}_2^{\rm T}{{{d}}_{{k}}}{r_k}({{\theta }}) + \right.$ $ \left.{\lambda _3}\left\{ {{r_k}({{\theta }}) - \tilde \sigma _T^2} \right\} \right]$

$\text{∂} {L_0}/\text{∂} \lambda _2^{} = 0,\;\;\;\text{∂} {L_0}/\text{∂}{\lambda _3} = 0,\;\;\;\text{∂} {L_0}/\text{∂} {{\theta }} = 0$ ,可得

$\left\{\begin{split} &\sum\limits_{k = 1}^N {w_k^ * {r_k}({{\theta }}){{{d}}_k} = 0\;} \\ &\sum\limits_{k = 1}^N {w_k^ * \left\{ {r_k^2({{\theta }}) - \tilde \sigma _T^2} \right\} = 0\;} \;\;\\ &\sum\limits_{k = 1}^N {w_k^ * } \left\{ {{{{d}}_k}{{d}}_k^{\rm T}{{{\lambda }}_2} + 2{r_k}({{\theta }}){{{d}}_k}} \right\} = 0\;\;\end{split}\right. $ (9)

由式(9)可以解得 ${{{\lambda }}_2} = {{\textit{0}}}$ ,所以,最终得到关于 ${\hat \lambda _3}$ ${\hat{ \theta }}$ $\hat w_k^ * $ 的方程组

$\left\{\begin{split}& \sum\limits_{k = 1}^N {\hat w_k^ * } \left\{ {r_k^2({\hat{ \theta }}) - \tilde \sigma _T^2} \right\} = 0\\ &\sum\limits_{k = 1}^N {\hat w_k^ * {r_k}({\hat{ \theta }}){{{d}}_k} = 0} \\ &\hat w_k^ * - \exp \left[ {{{\hat \lambda }_3}\left\{ {r_k^2({\hat{ \theta }}) - \tilde \sigma _T^2} \right\}} \right] = 0,\\ &\qquad k = 1,2, \cdots, N\end{split}\right.$ (10)

对于方程组(10),可将其转变成交替极大值和极小值问题。记

$ \begin{split} &\quad J({\lambda _3},{{\theta }}) = (1/N)\displaystyle\sum\limits_k {\exp \left[{\lambda _3}\left\{ {r_k^2({{\theta }}) - \tilde \sigma _T^2} \right\}\right]\;} \; =\qquad\qquad\;\;\\ &\qquad \left\{ \exp ({\lambda _3}\tilde \sigma _T^2)/N \right\}\displaystyle\sum\limits_k {\exp \left\{ {{\lambda _3}r_k^2({{\theta }})} \right\}} \end{split} $

${J_1}({\lambda _3},{{\theta }}) \!=\! \text{∂} J/\text{∂} {\lambda _3}$ ${J_{11}}({\lambda _3},{{\theta }}) \!=\! \text{∂} {J^2}/{\text{∂} ^2}{\lambda _3}$ ${J_2}({\lambda _3},$ ${{\theta }}) \!=\! \text{∂} J/\text{∂} {{\theta }}$

容易验证

$\left\{ \begin{array}{l} {J_1}({\lambda _3},{{\theta }}) = 0 \\ {J_2}({\lambda _3},{{\theta }}) = 0 \\ \end{array} \right.$ (11)

等式(11)与方程(10)等价。

容易求得 ${J_{11}}({\lambda _3},{{\theta }}) > 0$ ${J_1}(0,{{\theta }}) \leqslant 0$ ,所以,当固定 ${{\theta }}$ $J({\lambda _3},{{\theta }})$ 是关于 ${\lambda _3}$ 的凸函数,所以, $J({\lambda _3},{{\theta }})$ 一定能取到极小值。

现固定 ${{\theta }}$ ,求关于 ${\lambda _3}$ 的函数 $J({\lambda _3},{{\theta }})$ 取极小值时, ${\lambda _3}$ 的取值 ${\hat \lambda _3}({{\theta }})$ ,其对应点的一阶导数必为0,即

${J_1}({\hat \lambda _3}({{\theta }}),{{\theta }}) = 0。$ (12)

${\lambda _3} = {\hat \lambda _3}({{\theta }})$ 代入 $J({\lambda _3},{{\theta }})$ ,得到仅与 ${{\theta }}$ 有关的函数 $Q({{\theta }}) = J({\hat \lambda _3}({{\theta }}),{{\theta }})$ ,当 $Q({{\theta }})$ 取极大值时, ${{\theta }}$ 对应的取值记为 ${{\hat{ \theta }}_Q}$ ,则其对应点的一阶导数必为0,即

${J_1}({\hat \lambda _3}({{\hat{ \theta }}_Q}),{{\hat{ \theta }}_Q})*(\partial {\hat \lambda _3}({{\theta }})/\partial {{\theta }}) + {J_2}({\hat \lambda _3}({{\hat{ \theta }}_Q}),{{\hat{ \theta }}_Q}) = 0\!\!\!\!\!\!\!\!$ (13)

根据定义 ${J_1}({\hat \lambda _3}({{\hat{ \theta }}_Q}),{{\hat{ \theta }}_Q}) = 0$ ,式(13)可以表示为

$\left\{ \begin{array}{l} {J_1}({{\hat {\lambda} }_{3Q}},{{{\hat{ {\theta} }}}_Q}) = 0 \\ {J_2}({{{\hat {\lambda}} }_{3Q}},{{{\hat{ {\theta} }}}_Q}) = 0 \\ \end{array} \right.$ (14)

计算步骤总结如下:

a.取定 $\;\tilde \sigma _T^2\; = \min (\hat \sigma _{\rm {OLS}}^2,\;\hat \sigma _{\rm{LTS}}^2)\;$ $\;\hat \sigma _{\rm{LTS}}^2$ 是样本在截尾最小二乘估计下的残差尺度的估计。

b.选取截尾最小二乘估计 ${{\hat{ \theta }}^{(0)}}$ 作为 ${{\theta }}$ 的初始估计值,另外,给定 $\;\varepsilon $ 为一个很小的数。

c. ${{\theta }} = {{\hat{ \theta }}^{(0)}}$ ,求关于 ${\lambda _3}$ 的函数 $J({\lambda _3},{{\hat{ \theta }}^{(0)}})$ 的极小值,即 ${J_1}({\lambda _3},{{\hat{ \theta }}^{(0)}}) = 0$ 时对应的 ${\lambda _3}$ 的取值,记为 $\hat \lambda _3^{(0)} = \arg \min J({\lambda _3},{{\hat{ \theta }}^{(0)}})$

d.再固定 ${\lambda _3} = \hat \lambda _3^{(0)}$ ,求关于 ${{\theta }}$ 的函数 $J(\hat \lambda _3^{(0)},{{\theta }})$ 中的极大值,即 ${J_2}(\hat \lambda _3^{(0)},{{\theta }}) = 0$ 时对应 ${{\theta }}$ 的估计值 ${{\hat{ \theta }}^{(1)}}$ ,此时 ${{\hat{ \theta }}^{(1)}}$ 即为更新后的参数估计值。

e.重复步骤c—d,第 $i$ 次更新后的值记为 ${{\hat{ \theta }}^{(i)}}$ ,直至 $\left| {{{{\hat{ \theta }}}^{(i + 1)}} - {{{\hat{ \theta }}}^{(i)}}} \right| < \varepsilon $ ,或者,循环次数 $i > 1\;000$ 截止,此时 ${{\hat{ \theta }}^{(*)}}$ 记为 ${{{\theta }}^{(i)}}$ 的最终迭代值。

f.计算样本在 ${{\hat{ \theta }}^{(*)}}$ 下的残差 ${r_k}({{\hat{ \theta }}^{(*)}})$ ,则权重 ${w_k} = \exp \left\{ {\hat \lambda _3^{(0)}{r_k}({{{\hat{ \theta }}}^{(*)}})} \right\}\Biggr/\displaystyle\sum\limits_{k=1}^N {\exp \left\{ {\hat \lambda _3^{(0)}{r_k}({{{\hat{ \theta }}}^{(*)}})} \right\}} $ ,最终 ${{\theta }}$ 的估计值由式(7)可得,提取 ${\hat{ \theta }}$ 的前 $p$ 个分量为均值 $\,{{\hat{ \beta }}_0}$ ,最后再利用Cholesky分解,通过对 ${\hat{ \theta }}$ 的反向分解和还原,得到纵向数据协方差矩阵 ${{\hat{ \varSigma }}_0}$

3 模拟研究

根据模型 ${Y_{ij}} = {\beta _1} + {x_{ij1}}{\beta _2} + {x_{ij2}}{\beta _3} + {\varepsilon _{ij}}\;,\;i = 1,2, \cdots, $ $n, \;\;j = 1,2, \cdots, m$ 产生数据,设 $n = 100$ $m = 3$ ${x_{ij1}}\sim$ $ N(0,1),{x_{ij2}}\sim N(0,1)$ ,残差项 ${\varepsilon _{ij}}$ 服从均值为0,协方差结构分别为I单位结构、 $\,\rho = 0.5$ 的可交换结构、 $\,\rho = 0.5$ 的AR(1)和无结构的 $m$ 维正态分布。其中,无结构可定义为

$\left[ \begin{array}{l} \;1\;\;\;\;\;0.7\;\;\;0.4\;\; \\ 0.7\;\;\;\;\;1\;\;\;\;0.2 \\ 0.4\;\;\;\;0.2\;\;\;\;1 \\ \end{array} \right]$

$\,{{{\beta }}_0} = {({\beta _1},{\beta _2},{\beta _3})^{\rm T}} = {(1,1,1)^{\rm T}}$ ,模拟重复次数为100次。 ${{\theta }}$ 的初始估计值选取数据在最小二乘估计下的参数估计值 ${{\hat{ \theta }}_{\rm{OLS}}}$ ,记Xu为本文提出的方法,He是He等[5]提出的稳健广义估计方程的方法,Qin是Qin等[7]提出的稳健方程的估计方法,由于Qin等[7]研究的纵向数据带有缺失,而本文只考虑了污染,故在模拟时将缺失示性函数全部改为1,代表数据不缺失。此外,数据可能存在污染等异常情况,本文 $\,{{{\beta }}_0}$ 是最终估计式(7)中 ${\hat{ \theta }}$ 的前 $p$ 个分量,协方差 ${{\varSigma }}$ 通过返还函数对 ${\hat{ \theta }}$ 反向分解和还原可得,在返还函数中,R表示采用了绝对中位数离差MAD方法,NR表示一般处理方法,以此来验证本文提出的方法是否对异常值具备抵抗性。另外,两种方法中R和NR均代表各自的稳健方法和非稳健方法,表1比较了3种方法得到的 $\,{{\hat{ \beta }}_0}$ 的偏差(BIAS)、标准差(SE)和均方误差(MSE)。


表 1 数据不污染时参数估计对比 Table 1 Parameter estimation and comparison when data are not polluted

为了进一步说明本文提出方法的稳健性,对数据进行3%污染,污染方式为随机挑选样本3%的 ${x_i}$ 将其替换成 ${x_i} - 2$ ,并将 ${x_i}$ 所对应的 $y{}_i$ 替换成 ${y_i} + 2$ ,所得 $\,{{{\beta }}_0}$ 估计结果如表2所示。


表 2 数据3%污染时参数估计对比 Table 2 Parameter estimation and comparison when data are 3% polluted

通过比较可以看出:在数据不污染和3%污染的情况下,本文提出的稳健方法均要优于其他2种方法,参数的偏差和均方误差都比其他2种方法要小,显示估计更准确;当数据不污染时,本文提出的方法与其他方法比较,主要是参数的偏差比较小,从而使得本文方法的有效性比其他2种方法高;当数据被污染时,本文提出的方法的有效性更为显著。相同条件下,本文提出方法的偏差远远小于其他2种方法,以参数 $\,{\beta _3}$ 为例,固定协方差为I单位结构,在数据不污染时,相同稳健的方法,本文 $\,{\beta _3}$ 偏差分别比Qin和He小0.000 2,0.006;在数据3%污染时,相同稳健的方法,本文 $\,{\beta _3}$ 偏差分别比Qin和He小0.001,0.011 1。

此外,本文还对数据进行了5%的污染,结果显示,数据污染程度越大,本文的方法越有效,为了节省篇幅,相关表格在本文没有给出。

为了进一步说明对纵向数据协方差矩阵的估计好坏程度,本文定义指标QLEL来衡量 ${\hat{ \varSigma }}$ ${{\varSigma }}$ 的接近程度,QLEL越小,表明协方差矩阵的估计值越接近真实值。

$QL = {\varDelta _Q}({{\varSigma }},{\hat{ \varSigma }}) = {\left( {{{trace}}({{{\varSigma }}^{ - 1}}{\hat{ \varSigma }} - {{I}})} \right)^2}$
$EL = {\varDelta _E}({\hat{ \varSigma }},{{\varSigma }}) = {{trace}}({{{\varSigma }}^{ - 1}}{\hat{ \varSigma }}) - \log \left| {{{{\varSigma }}^{ - 1}}{\hat{ \varSigma }}} \right|$

式中: $\varDelta $ 为协方差真实值与估计值的差值;trace为矩阵的迹。

在数据不污染和3%污染的情况下,将本文方法和Qin等[7]提出的稳健方程进行比较,对应的QLEL结果如表3所示,其中, ${{EL}}_i,{{QL}}_i(i =1,2,$ $ 3,4)$ 代表误差项服从第 $i$ 种协方差结构下的QLEL指标,对应模型中设定的4种协方差结构。


表 3 ELQL对比表 Table 3 EL and QL comparison table

表3可以看出:在数据不污染时,本文提出的非稳健方法所估计出的协方差的EL指标,一直优于Qin提出的方法;在数据样本受到3%污染时,本文提出的稳健方法所估计的协方差的ELQL指标,都几乎一直优于Qin提出的方法,表明本文方法估计出的协方差结构与真实值最为接近,说明本文提出的方法具有显著稳健性和有效性,不仅能保证参数 $\,{{{\beta }}_0}$ ${{\varSigma }}$ 估计的有效性,而且在样本受到污染时仍能保证其准确性;此处还对数据进行了5%的污染,得到结果与上面类似,并且污染程度越大,本文方法的优越性越显著。

为了进一步比较经验似然方法的优越性,选取误差项服从多元t分布的纵向数据进行模拟研究,研究方法同上,所得结果如表4所示。


表 4 t分布误差下的参数估计对比表 Table 4 Comparison of parameter estimates under t distribution error

表4可知,在误差服从t分布时,经验似然方法的估计结果仍然非常可观。3种方法的稳健结果显示,在4种协方差结构下,本文方法的估计方差要远远小于其他2种方法,说明本文提出的经验似然估计方法比其他2种方法更有效,参数估计更趋于一致。此外,对比误差项服从多元t分布下本文方法和Qin等[7]的4种协方差估计情况如表5所示。


表 5 t分布误差下QLEL Table 5 EL and QL under t distribution error

表5可知,在误差服从多元t分布时,经验似然方法的有效性更为显著。根据QL指标和EL指标显示,本文估计的纵向数据协方差矩阵与真实协方差矩阵相似度更高,与表3对比可看出,当模拟数据的误差服从多元t分布时,其协方差矩阵的估计效果更优于误差服从多元正态分布时的情形,进一步说明了本文方法对误差具有抵抗性。

4 实证分析

为了进一步验证本文方法的实用性,分析了CD4细胞计数纵向数据。关于这个数据集的完整描述见Diggle的主页: http://www.lancs.ac.uk/diggle/. 这里分析所用的数据是完整数据的截取部分, 包含340个样本的1 020次观测值。响应变量 $y$ 是CD4计数的算数平方根, 协变量包括血清转换的时间 ${x_2}$ ,相对于一个起点的年龄 ${x_3}$ ,由吸烟的包数刻画的吸烟状况 ${x_4}$ ,娱乐性药物使用是/否 ${x_5}$ ,性伴侣的个数 ${x_6}$ ,以及流行病中心给出抑郁状态和抑郁程度 ${x_7}$ 。有许多学者研究过这个数据集。Zeger等[10]和Wang等[11]分别对这个数据拟合了半参数模型和部分线性半参数模型, 其中,协变量 ${x_2}$ 作为非线性形式进入模型。

本文分析的主要目的是寻找CD4数量和协变量之间的关系,考虑到可能存在非线性关系,引入了除 ${x_5}$ 以外的平方项,加上截距项 ${x_1}$ ,记

$\quad\;\;{{X}} = {({x_1},{x_2},x_2^2,{x_3},x_3^2,{x_4},x_4^2,{x_5},{x_{6,}}x_6^2,{x_7},x_7^2)^{\rm T}}\!\!\!\!$ (15)

建立线性模型

$\qquad\quad{y_{ij}} = {{{X}}_{ij}}^{\rm T}{{\beta }}\;\;\;,i = 1,2 ,\cdots, 340,\;\;\;j = 1,2,3$ (16)

由于实际数据无法知道真实参数,所以,现利用交互验证的方法来比较各种方法的优劣,交互验证过程中的均方误差

$MS\!E_{CV} = \frac{1}{n}\sum\limits_{i = 1}^n {\left| {{y_i} - {{\hat y}_{( - i)}}} \right|} $ (17)

这里的 $n = 339$ ${\hat y_{( - i)}}$ 是去掉第 $i$ 个个体外,由其余339个个体拟合所得的预测值,数据预处理后,代入算法中,本文方法与Qin和He方法的比较结果如表6所示。

表6说明,本文方法在实际应用中十分有效,CV值要明显小于另外2种方法的CV值,值得推广。


表 6 CV值比较 Table 6 Comparison of CV value
参考文献
[1]
THOMAS D R, GRUNKEMEIER G L. Confidence interval estimation of survival probabilities for censored data[J]. Journal of the American Statistical Association, 1975, 70(352): 865-871. DOI:10.1080/01621459.1975.10480315
[2]
OWEN A. Empirical likelihood ratio confidence regions[J]. The Annals of Statistics, 1990, 18(1): 90-120. DOI:10.1214/aos/1176347494
[3]
OWEN A. Empirical likelihood for linear models[J]. The Annals of Statistics, 1991, 19(4): 1725-1747. DOI:10.1214/aos/1176348368
[4]
BONDELL H D, STEFANSKI L A. Efficient robust regression via two-stage generalized empirical likelihood[J]. Journal of the American Statistical Association, 2013, 108(502): 644-655. DOI:10.1080/01621459.2013.779847
[5]
HE X M, FUNG W K, ZHU Z Y. Robust estimation in generalized partial linear models for clustered data[J]. Journal of the American Statistical Association, 2005, 100(472): 1176-1184. DOI:10.1198/016214505000000277
[6]
QIN G Y, ZHU Z Y, FUNG W K. Robust estimation of generalized partially linear model for longitudinal data with dropouts[J]. Annals of the Institute of Statistical Mathematics, 2016, 68(5): 977-1000. DOI:10.1007/s10463-015-0519-8
[7]
QIN G Y, ZHU Z Y. Robust estimation of mean and covariance for longitudinal data with dropouts[J]. Journal of Applied Statistics, 2015, 42(6): 1240-1254. DOI:10.1080/02664763.2014.999033
[8]
CRESSIE N, READ T R C. Multinomial goodness-of-fit tests[J]. Journal of the Royal Statistical Society. Series B, 1984, 46(3): 440-464.
[9]
CORCORAN S A. Bartlett adjustment of empirical discrepancy statistics[J]. Biometrika, 1998, 85(4): 967-972. DOI:10.1093/biomet/85.4.967
[10]
ZEGER S L, DIGGLE P J. Semiparametric models for longitudinal data with application to CD4 cell numbers in HIV seroconverters[J]. Biometrics, 1994, 50(3): 689-699. DOI:10.2307/2532783
[11]
WANG Y G, LIN X, ZHU M. Robust estimating functions and bias correction for longitudinal data analysis[J]. Biometrics, 2005, 61(3): 684-691. DOI:10.1111/j.1541-0420.2005.00354.x