纵向数据部分线性模型的一种稳健估计方法

A Robust Estimation Method for the Partial Linear Models with Longitudinal Data

摘要: 本文利用指数平方损失函数，针对纵向数据部分线性模型的参数部分提出了一种稳健的估计方法。首先通过核估计以及变换将非参部分消去，然后利用指数平方损失函数降低异常值的影响，构造出稳健的广义估计方程，最终推导出稳健的经验似然比函数，由此进行估计和统计推断。在一定正则条件下该估计具有相合性以及渐近正态性，蒙特卡洛模拟显示本文的方法具有稳健性以及较高的效率。

Abstract: Based on the exponential squared loss function, we propose a robust estimation method for the parametric parts in the partial linear models. By using the kernel approximation and transfor-mation, we first absorb the nonparametric part. Then by using the exponential squared loss func-tion to deduce the influence of outliers, we propose robust generalized estimation equations and empirical likelihood ratio function for estimation and inference. Under some regularity conditions, we show that the resulting estimators are consistent and asymptotic normality. Simulation results also illustrate the robustness and efficiency advantages of our method.

1. 引言

部分线性模型是一类非常重要的半参数模型，比单纯的线性回归模型或非参数模型有更大的适应性以及更强的解释能力，在医学和计量经济学等领域有着广泛的应用。而在这些领域中常常用到的一种数据类型为纵向数据，即每个个体经过多次测量得到的数据集。考虑来自n个个体的数据，其中第i个个体有 $m_i$ 次观测， $i=1,\cdots ,n$ ，总的观测数为 $N={\displaystyle {\sum }_{i=1}^n{m}_i}$ 。设 $Y_{ij}$ 和 $\left(X_{ij},T_{ij}\right)$ 分别是是第i个个体第j次观测的响应变量和协变量，其中 $X_{ij}$ 是 $p\times 1$ 维向量， $Y_{ij}$ 和 $T_{ij}$ 都是一维变量。纵向数据的部分线性模型为：

$Y_{ij}=X_{ij}^{\text{T}}\beta +g\left(T_{ij}\right)+{\epsilon }_{ij},i=1,2,\cdots ,n;j=1,2,\cdots ,m_i$ (1)

其中 $\beta$ 是 $p\times 1$ 维未知参数向量， $g(\cdot )$ 是未知的基准函数， ${\epsilon }_{ij}$ 为随机误差项，假定 $E\left[{\epsilon }_{ij}|T_{ij}\right]=\text{0}$ ， $E\left[{\epsilon }_{ij}^{\text{2}}|T_{ij}\right]={\sigma }_{\epsilon }^{\text{2}}$ ，不失一般性，不妨假定 $T_{ij}\in \left[0,1\right]$ 。一般假定来自不同个体的观测是相互独立的(称之为组间独立)，而来自同一个个体的观测具有一定的相关性(称之为组内相关)。

Zeger和Diggle [1] (1994)最早对模型(1)进行了研究，利用后移算法(back-fitting)，结合核估计以及最小二乘估计将参数 $\beta$ 和函数 $g(\cdot )$ 估计出来。Lin和Carroll [2] (2001)利用profile估计方程估计参数部分，并结合估计方程与核估计方法估计非参部分。He等 [3] (2002)利用B-样条逼近非参部分，然后选择适当的损失函数用M估计方法估计参数 $\beta$ 以及样条系数。Fan和Li [4] (2004)利用局部多项式逼近非参部分，在假定工作独立协方差矩阵下，提出两种估计方法：差分估计和profile最小二乘估计。Xue和Zhu [5] (2008)利用核估计以及差分法消去未知的基准函数 $g(\cdot )$ ，并利用经验似然(Owen, 1988) [6] 方法估计 $\beta$ 以及进行统计推断。以上方法都没有考虑纵向数据模型中的组内相关性以及估计的稳健性，而在对实际问题的研究中得到的数据常常存在异常值或与假定的分布不符，于是估计的稳健性变得尤为重要。

本文首先利用广义估计方程(Liang, 1986) [7] 的思想，引入工作相关矩阵将组内相关性考虑进来，然后利用指数平方损失函数(Wang, 2013) [8] 降低异常值的影响，并在估计方程中加W权重减轻杠杆值的影响，得到稳健的广义估计方程，并以此构造出稳健的经验似然函数，得到稳健的经验似然估计，并进行统计推断。

2. 稳健经验似然的构造

首先对(1)式两边取给定 $T_{ij}$ 下的条件期望，得

$E\left[Y_{ij}|T_{ij}\right]=E\left[X_{ij}^{\text{T}}|T_{ij}\right]\beta +g\left(T_{ij}\right)$ ，

两式相减得

$Y_{ij}-E\left[Y_{ij}|T_{ij}\right]={\left(X_{ij}-E\left[X_{ij}|T_{ij}\right]\right)}^{\text{T}}\beta +{\epsilon }_{ij}$ ，

令 $Y_{ij}^{*}=Y_{ij}-E\left[Y_{ij}|T_{ij}\right]$ ， $X_{ij}^{*}=X_{ij}-E\left[X_{ij}|T_{ij}\right]$ ，则模型变为

$Y_{ij}^{*}=X_{ij}^{*\text{T}}\beta +{\epsilon }_{ij},i=1,2,\cdots ,n;j=1,2,\cdots ,m_i$ 。

为了稳健，本文引入指数平方损失函数 ${\phi }_{\gamma }\left(r\right)=-\frac{2r}{\gamma }\exp \left(-\frac{r^2}{\gamma }\right)$ ，其中 $\gamma >0$ 为调节参数，当 $\gamma$ 非常大时， $1-\exp \left(-r^2/\gamma \right)\approx r^2/\gamma$ ，于是得到的估计与最小二乘估计类似；当 $\gamma$ 较小时，指数平方损失函数能够降低较大的残差对参数估计的影响，即选取较小的 $\gamma$ 能够有效的控制异常值的影响，提高估计的稳健性，最优 $\gamma$ 值的选取见第三节。

令 $Y_i^{*}={\left(Y_{i1}^{*},Y_{i2}^{*},\cdots ,Y_{i{m}_i}^{*}\right)}^{\text{T}}$ ， $X_i^{*}={\left(X_{i1}^{*},X_{i2}^{*},\cdots ,X_{i{m}_i}^{*}\right)}^{\text{T}}$ ， ${\epsilon }_i={\left({\epsilon }_{i1},{\epsilon }_{i2},\cdots ,{\epsilon }_{i{m}_i}\right)}^{\text{T}}$ ， $r_i^{*}\left(\beta \right)=Y_i^{*}-X_i^{*}\beta$ ， $V_i=Var\left({\epsilon }_i\right)={\sigma }_{\epsilon }^2{R}_i\left(\rho \right)$ ， $R_i\left(\rho \right)$ 为工作相关矩阵， $\rho$ 为未知的参数，注意在随后的矩条件 $E\left[Z\left({\beta }_0\right)\right]=0$ 中可将方差 ${\sigma }_{\epsilon }^2$ 消去，于是不妨直接假定 $V_i=R_i\left(\rho \right)$ 。构造辅助随机向量

$Z_i\left(\beta \right)=X_i^{*\text{T}}{V}_i^{-1}{h}_i^{\gamma }\left(r_i^{*}\left(\beta \right)\right),i=1,\cdots ,n$ 。

其中 $h_i^{\gamma }\left(r_i^{*}\left(\beta \right)\right)=C_i\left[{\phi }_{\gamma }\left(r_i^{*}\left(\beta \right)\right)-E\left[{\phi }_{\gamma }\left(r_i^{*}\left(\beta \right)\right)\right]\right]$ ，减去 ${\phi }_{\gamma }$ 的期望是为了确保Fisher一致性(He [9] , 2005)。权重矩阵 $C_i=diag\left(c_{i1},\cdots ,c_{i{m}_i}\right)$ 用来约束协变量空间中杠杆点(leverage points)的影响，常用的权重是Mallows-type权函数

$c_{ij}=c\left(x_{ij}\right)=\min \left\{1,{\left\{\frac{b_x}{{\left(x_{ij}-m_x\right)}^{\text{T}}{S}_x^{-1}\left(x_{ij}-m_x\right)}\right\}}^{\frac{\xi }{2}}\right\}$ ，

其中 $\xi \ge 1$ ， $m_x$ 和 $S_x$ 分别是对 $x_{ij}$ 位置参数和尺度参数的稳健估计， $b_x$ 是 ${\chi }^2\left(p\right)$ 的95%分位数。

利用辅助向量 $Z_i\left(\beta \right)$ 可以定义 $\beta$ 的稳健经验似然比函数，然而此函数不能直接用于估计 $\beta$ 以及进行统计推断，因为 $Z_i\left(\beta \right)$ 中含有未知的函数 $E\left[Y_{ij}|T_{ij}\right]$ 和 $E\left[X_{ij}|T_{ij}\right]$ 以及未知的相关性参数 $\rho$ ，需要用各自的估计量替代。对于 $\rho$ 可以用稳健的矩估计方法，例如，对于可交换工作相关结构， $\rho$ 的估计为

$\hat{\rho }=\frac{1}{n{H}^2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{1}{m_i\left(m_i-1\right)}}{\displaystyle \underset{j\ne k}{\sum }{\phi }_{\gamma }\left(e_{ij}\right){\phi }_{\gamma }\left(e_{ik}\right)}$ ，

对于一阶自回归工作相关结构， $\rho$ 的估计为

$\hat{\rho }=\frac{1}{n{H}^2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{1}{m_i-1}}{\displaystyle \underset{j\le m_i-1}{\sum }{\phi }_{\gamma }\left(e_{ij}\right){\phi }_{\gamma }\left(e_{i,j+1}\right)}$ ，

其中 $H^2=E\left[{\phi }_{\gamma }^2\left(e_{ij}\right)\right]$ 。

对于 $E\left[Y_{ij}|T_{ij}\right]$ 和 $E\left[X_{ij}|T_{ij}\right]$ 的估计可用非参数统计方法，本文选用核估计

$\hat{E}\left[Y_{ij}|T_{ij}=t\right]={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}{W}_{ij}\left(t\right)Y_{ij}}}$ 和 $\hat{E}\left[X_{ij}|T_{ij}=t\right]={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}{W}_{ij}\left(t\right)X_{ij}}}$ ，其中 $W_{ij}\left(t\right)=K_h\left(T_{ij}-t\right)/{\displaystyle \underset{k=1}{\overset{n}{\sum }}{\displaystyle \underset{l=1}{\overset{m_k}{\sum }}{K}_h\left(T_{kl}-t\right)}}$ ， $K_h(\cdot )=K\left(\cdot /h\right)$ ， $K(\cdot )$ 是一个核函数。

用估计值分别替代 $Z_i\left(\beta \right)$ 中的 $\rho$ 、 $E\left[Y_{ij}|T_{ij}\right]$ 和 $E\left[X_{ij}|T_{ij}\right]$ ，得到 $Z_i\left(\beta \right)$ 的估计量 ${\hat{Z}}_i\left(\beta \right)$ ，以此构建 $\beta$ 的稳健似然比函数

$L\left(\beta \right)=\max \left\{{\displaystyle \underset{i=1}{\overset{n}{\prod }}n{p}_i}|p_i\ge 0,{\displaystyle \underset{i=1}{\overset{n}{\sum }}{p}_i=1,{\displaystyle \underset{i=1}{\overset{n}{\sum }}{p}_i{\hat{Z}}_i\left(\beta \right)=0}}\right\}$ ，

利用拉格朗日乘子法可以得出， $L\left(\beta \right)$ 的最大值在 $p_i=\frac{1}{n\left(1+{\lambda }^{\text{T}}{\hat{Z}}_i\left(\beta \right)\right)}$ ， $i=1,\cdots ,n$ 处取到，

其中 $\lambda$ 满足

${\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{{\hat{Z}}_i\left(\beta \right)}{1+{\lambda }^{\text{T}}{\hat{Z}}_i\left(\beta \right)}}=0$ ，

于是 $\beta$ 的对数经验似然比函数为

$l\left(\beta \right)=-{\displaystyle \underset{i=1}{\overset{n}{\sum }}\log \left(1+{\lambda }^{\text{T}}{\hat{Z}}_i\left(\beta \right)\right)}$ 。

最大化似然比函数可以求出 $\beta$ 的估计 ${\hat{\beta }}_{\text{MELE}}$ ，根据Xue [5] 和He [9] ，在一定正则条件下，有

$\sqrt{n}\left({\hat{\beta }}_{MELE}-{\beta }_0\right)\stackrel{d}{\to }N\left(0,{\Sigma }_0^{-1}{\Omega }_0{\Sigma }_0^{-1}\right)$ ，

以及

$-2l\left(\beta \right)\stackrel{d}{\to }{\chi }^2\left(p\right)$ ，

其中 ${\chi }^2\left(p\right)$ 表示自由度为p的卡方分布，于是本文估计具有相合性以及渐近正态性，且 $\beta$ 的置信域为

$I_{\alpha }=\left\{\beta |-2l\left(\beta \right)\le {\chi }_{1-\alpha }^2\left(p\right)\right\}$ ，

其中 $\alpha$ 为显著性水平。

3. 算法

下面给出具体的算法：

Step 1：用核估计计算 $\hat{E}\left[Y_{ij}|T_{ij}\right]$ 和 $\hat{E}\left[X_{ij}|T_{ij}\right]$ ，对数据进行变换， ${\hat{X}}_{ij}^{*}=X_{ij}-\hat{E}\left[X_{ij}|T_{ij}\right]$ ；

Step 2：用GEE估计作为参数估计的初始值 ${\beta }^{\left(0\right)}$ ，令t=1；

Step 3：在估计值为 ${\beta }^{\left(t\right)}$ 时，求出相关性参数的估计 $\hat{\rho }$ ，并代入工作相关矩阵得到 ${\hat{V}}_i=R_i\left(\hat{\rho }\right)$ 。通过最小化 $\det \left(Cov\left({\beta }^{\left(t\right)}\right)\right)$ 得到调节参数 $\gamma$ 的估计 ${\gamma }_{opt}^{\left(t\right)}$ ，其中 $\det (\cdot )$ 表示行列式运算；

$Cov\left({\beta }^{\left(t\right)}\right)={\left[\hat{\Sigma }\left({\mu }_i\left({\beta }^{\left(t\right)}\right),\hat{\rho }\right)\right]}^{-1}\hat{\Omega }\left({\mu }_i\left({\beta }^{\left(t\right)}\right),\hat{\rho }\right){\left[\hat{\Sigma }\left({\mu }_i\left({\beta }^{\left(t\right)}\right),\hat{\rho }\right)\right]}^{-1}$ ，

$\hat{\Sigma }\left({\mu }_i\left({\beta }^{\left(t\right)}\right),\hat{\rho }\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{X}_i^{\text{T}}{\hat{V}}_i^{-1}}E\left[\frac{\partial h_i^{\gamma }\left({\mu }_i\left({\beta }^{\left(t\right)}\right)\right)}{\partial {\mu }_i\left({\beta }^{\left(t\right)}\right)}\right]X_i$ ，

$\hat{\Omega }\left({\mu }_i\left({\beta }^{\left(t\right)}\right),\hat{\rho }\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{X}_i^{\text{T}}}{\hat{V}}_i^{-1}\left[h_i^{\gamma }\left({\mu }_i\left({\beta }^{\left(t\right)}\right)\right)\right]{\left[h_i^{\gamma }\left({\mu }_i\left({\beta }^{\left(t\right)}\right)\right)\right]}^{\text{T}}{\hat{V}}_i^{-1}{X}_i$ ；

Step 4：最小化 $-2l\left(\beta ,\hat{\rho }\right)$ 得到新的估计值 ${\beta }^{\left(t+1\right)}$ ；

Step 5：令t = t + 1，迭代Step 2~Step 3直到收敛，得到最终的估计值 ${\hat{\beta }}_{MELE}$ 和 ${\hat{\gamma }}_{opt}$ 。

4. 模拟

本节通过蒙特卡洛模拟，比较本文方法与Xue提出方法的稳健性。模型设定为

$y_{ij}=x_{ij}\beta +g\left(t_{ij}\right)+{\epsilon }_{ij},i=1,2,\cdots ,n;j=1,2,\cdots ,m$ ，

其中 $x_{ij}~Uniform\left(-1,1\right)$ ， $t_{ij}~Uniform\left(0,1\right)$ ， $g\left(t\right)=3\sin \left(\frac{\text{π}}{2}t\right)$ ， ${\left({\epsilon }_{i1},\cdots ,{\epsilon }_{iT}\right)}^{\text{T}}~N\left({\left(0,\cdots ,0\right)}^{\text{T}},{\sigma }^{2}R\left(\rho \right)\right)$ ，令 ${\sigma }^{\text{2}}=1$ ， $R\left(\rho \right)$ 分别取可交换结构(Exch)和一阶自相关结构(AR-1)，参数 $\rho$ 选取0.3和0.7两个数值，取n = 50，100；m = 3，10。估计时所用的工作相关性矩阵分别选取独立结构(Ind)、可交换结构(Exch)和一阶自相关结构(AR-1)。令数据以概率 $\delta$ 受到来自分布 ${\chi }^{\text{2}}\left(5\right)-5$ 的污染，考虑 $\delta =0.05,0.1$ 和0.2三种不同的污染程度。由n、m、 $\rho$ 和 $\delta$ 的不同取值以及真实相关矩阵的不同结构共得到24种不同的情形，每种情形分别进行1000次模拟，计算出对应的偏(Bias)、均方误差(MSE)、覆盖率(CP)以及置信区间宽度(Width)，结果列在下文的表格和图中。

表1所列为四种估计的偏(Bias)与均方误差(MSE)的计算数值，其中 $\hat{\beta }\left(\text{Xue}\right)$ 表示用Xue的方法得到

Table 1. The bias and mean squared error of estimators

表1. 估计的偏(Bias)与均方误差(MSE)

的估计， $\hat{\beta }\left(\text{Exch}\right)$ 、 $\hat{\beta }\left(\text{AR-1}\right)$ 和 $\hat{\beta }\left(\text{Ind}\right)$ 表示本文方法利用不同工作相关性矩阵结构得到的估计，第一行的Exch和AR-1表示真实的相关性矩阵的结构。通过比较可以发现，在多数情形下，本文估计量(不论利用何种工作相关矩阵结构)的偏差要更加接近于0；而在所有情形下，本文估计量的均方误差要远远小于Xue的估计，大多不到其一半，这表明本文构建的估计量有较强的稳健性。当污染率增加时，四种估计的均方误差也随之增加，这是由于污染数据的比例增加导致估计结果的波动增大，而相对于Xue的方法本文估计量的均方误差仍处于较低的水平。当n或m增大时，四种估计的均方误差都减小，这是由于样本量增大导致估计的方差减小。

图1表示真实的相关性矩阵结构为Exch时的不同情形下四种估计的均方误差，其中圆圈符号代表Xue的估计，三角形代表本文提出的估计，显然在所有情形下本文方法的均方误差都远小于Xue的方法。图1中的前3个三角表示在n = 50，m = 3， $\rho =0.3$ ，污染率为0.05时分别利用Exch、AR-1和Ind三种工作矩阵结构得到的估计，可以发现工作矩阵用Exch结构得到的估计的均方误差最小，用AR-1结构得到的结果与之相近，用Ind结构(即不考虑组内相关性)得到的估计的均方误差最大。前3个圆圈表示对应情形下Xue的估计，由于Xue的估计没有用到工作相关矩阵，因此得到的估计相同，于是这三个圆圈处

于相同的高度。从圆圈的分布可以看出3个圆圈为一个阶梯，第1~3个阶梯所处情形的区别在于污染率分别为0.05、0.1和0.2，于是均方误差呈现上升的趋势；前3个阶梯与随后的3个阶梯所处情形的区别在于相关性参数的不同，因此这两段阶梯形状和高度相似；前6个阶梯与随后的6个阶梯所处情形的区别在于m分别取3和10，m取10时样本量增加了大约两倍，估计的均方误差更小，于是第7~12个阶梯要明显低于前6个阶梯。对三角形的走势进行同样的分析也能得出类似的结果，而且在每个三角形的阶梯中，总是第1个三角形最低，第3个三角形最高，这表明当模型中存在组内相关性而在估计时忽略相关性(用Ind工作矩阵)得到的估计效率最低，而选取的工作矩阵结构为真实的结构时估计的效率最高。

表2列出的是n = 50，m = 3时不同情形下四种方法求出的置信水平为95%的置信区间的覆盖率(CP)

Figure 1. Mean squared error of four estimators under different conditions

图1. 不同情形下四种估计的均方误差

Table 2. Coverage probabilities and width of the confidence interval

表2. 置信区间的覆盖率(CP)和宽度(Width)

和宽度(Width)，第一行的Exch和AR-1表示真实的相关性矩阵的结构。通过比较可以发现，在几乎所有情形下，与用Xue的方法求出的置信区间相比，本文方法(不论用何种工作相关矩阵结构)得到的置信区间的覆盖率更加接近真实的置信水平95%，而且宽度更小，这表明用本文方法进行统计推断能够得到更加稳健的结果。通过比较利用三种不同的工作相关结构得到的结果发现，当模型中存在组内相关性而在统计推断时忽略相关性(用Ind工作矩阵)得到的结果最差，而选取的工作矩阵结构为真实的结构时统计推断的结果最好。从表2中还可以看出，污染率由0.05增大到0.2时，置信区间的宽度也相应增加，这是由于数据污染越严重，估计量的波动越大，达到相同的置信率所需要的置信区间也随之扩大。

5. 结论

本文利用指数平方损失函数，针对纵向数据部分线性模型的参数部分提出了一种稳健的估计方法。首先通过核估计以及变换将非参部分消去，然后利用广义估计方程(GEE)的思想，通过假定相关性结构将组内相关性考加入到参数估计过程中；同时为了提高估计的稳健性，引入指数平方损失函数降低异常值的影响，并利用W权重减轻杠杆值的影响，构造出稳健的广义估计方程，最终推导出稳健的经验似然比函数，通过最大化似然比函数可得到参数的稳健估计，在一定正则条件下该估计具有相合性以及渐近正态性，且似然比函数服从渐近卡方分布，由此可进行稳健的统计推断。

本文在生成的数据中加入一定比例的污染，考虑了不同样本容量、不同相关性结构以及不同的污染率，对24种情形分别进行蒙特卡洛模拟，并将本文提出的方法与Xue的方法进行比较。结果表明，多数情形下本文估计量(不论利用何种工作相关结构)的偏差要更加接近于0，且在所有情形下的均方误差远小于Xue估计的均方误差，这表明本文构建的估计量有较强的稳健性。将所有情形得到的结果放在一张图上可以明显看出，若存在组内相关性而估计时忽略，估计的有效性最差；而工作相关矩阵用可交换结构(Exch)和一阶自相关结构(AR-1)得到的估计的差别相对较小；若假定的工作相关性矩阵的结构与真实结构相同，估计的效果最好。与用Xue方法求得的置信区间相比，不论工作相关性矩阵的结构是否准确，本文方法得到的置信区间有着更小的宽度以及更加接近真实置信水平的覆盖率，这表明用本文方法进行统计推断能够得到更加稳健的结果。

基金项目

本项目受国家自然科学基金(11231005; 71673171)和山东省自然科学基金(ZR2017BA002)资助。