因变量缺失下线性回归模型的岭估计

Ridge Estimation of Linear Regression Model with Missing Responses

摘要: 本文研究线性回归模型存在多重共线性且因变量存在缺失时的估计问题，分别基于完整数据方法和单点插补方法，给出了回归系数的两种岭估计，给出了估计量的渐近性质。最后通过数值模拟验证了所提方法的有效性。

Abstract: This paper discusses estimation of linear regression models in the presence of multicollinearity and the responses are missing at random. Based on the complete-case method and single imputation technique, two ridge estimators for the unknown regression coefficients are proposed. Asymptotically properties of the proposed estimators are given. Finally, some simulations are conducted to illustrate the proposed methods.

1. 引言

缺失数据是数据分析中遇到的常见问题之一，关于缺失数据的详细介绍可参考著作Little和Ruin (2002) [1]。线性回归模型是最为常用的统计方法之一，缺失数据下回归模型的研究得到了关注，Cheng (1994) [2]、Chu & Cheng (1995) [3]、Wang & Rao (2002) [4] 和Qin等(2009) [5] 基于不同的方法研究了因变量缺失下的回归模型。另一方面，多种共线性是线性回归模型实际使用中经常遇到的问题，这一问题是由于自变量之间存在线性关系导致，严重影响统计推断结果。解决多重共线性的途径有多重，其中构造回归系数的有偏估计是一种常见的方法，有偏估计以牺牲无偏性使估计量的均方误差变小，从而解决了多重共线性下最小二乘估计量方差过大的问题。目前讨论较多的有偏估计是岭估计和主成份估计。

线性回归模型可记为如下形式

$Y_i=X_i^{\text{T}}\beta +{\epsilon }_i,i=1,2,\cdots ,n$ (1)

其中 $Y_i$ 是因变量观测值， $X_i={\left(X_{i1},X_{i2},\cdots ,X_{iq}\right)}^{\text{T}}$ 是相应的自变量观测值， $\beta$ 为 $q\times 1$ 需要估计的未知回归系数，模型误差 ${\epsilon }_i$ 为独立同分布的随机变量，有 $E\left({\epsilon }_i|X_i\right)=0$ 和 $Var\left({\epsilon }_i|X_i\right)={\sigma }^2$。为了克服因变量的缺失，用变量 $\delta$ 作为缺失的标志， ${\delta }_i=1$ 表示 $Y_i$ 的值没有缺失，可以被观测到， ${\delta }_i=0$ 则表示 $Y_i$ 值缺失。假定 $Y_i$ 满足随机缺失的机制，即

$P\left(\delta =1|X,Y\right)=P\left(\delta =1|X\right).$ (2)

该缺失机制是缺失数据分析中常用的假设条件。

针对缺失数据的处理，最简单的处理方法就是完整数据方法，即只利用因变量和自变量都完整的观测值，也就是只考虑因变量不存在缺失的那些观测值，即 ${\delta }_i=1$ 的那些观测数据。显然，这种方法虽然简便，但舍弃了存在缺失的数据，造成了信息的浪费。为了充分利用观测数据的信息，还可以利用插补的方法对模型(1)进行估计。对于模型(1)~(2)，基于完整数据方法和单点插补方法，杨徐佳等(2011) [6] 构造了回归系数的估计，并给出了所提估计量的渐近性质，此外还讨论了回归系数的线性关系检验问题。安佰玲等(2013) [7] 则基于这两种方法讨论了模型的约束估计问题。

针对因变量缺失和自变量存在多重共线性这两个问题，目前的研究大都是单独讨论，将两个问题同时考虑的研究成果很少。为了解决这一问题，本文集中讨论因变量缺失下线性回归模型的岭估计问题。

论文第2节和第3节将分别基于完整数据方法和单点插补方法构造模型系数的岭估计，并给出估计量的渐近性质。第4节将通过数值模拟验证所提方法的有效性。总结将在第5节给出，定理的证明将放在第6节。

2. 基于完整数据方法的岭估计

本节基于完整数据方法构造模型系数的岭估计。设 ${\{Y_i,{\delta }_i,X_i\}}_{i=1}^n$ 为来自模型(1)的观测数据，则有如下的线性回归模型

${\delta }_i{Y}_i={\delta }_i{X}_i^{\text{T}}\beta +{\delta }_i{\epsilon }_i,i=1,2,\cdots ,n.$ (4)

模型(4)的矩阵形式为

$\bar{Y}=\bar{X}\beta +\bar{\epsilon },$ (5)

其中 $\bar{Y}={\left({\delta }_1{Y}_1,{\delta }_2{Y}_2,\cdots ,{\delta }_n{Y}_n\right)}^{\text{T}}$， $\bar{X}={\left({\delta }_1{X}_1,{\delta }_2{X}_2,\cdots ,{\delta }_n{X}_n\right)}^{\text{T}}$， $\bar{\epsilon }={\left({\delta }_1{\epsilon }_1,{\delta }_2{\epsilon }_2,\cdots ,{\delta }_n{\epsilon }_n\right)}^{\text{T}}$。

从而对模型(4)使用最小二乘估计，可得回归系数 $\beta$ 基于完整数据的估计

${\hat{\beta }}_C={\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}{\delta }_i{X}_i{X}_i^{\text{T}}}\right)}^{-1}\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}{\delta }_i{X}_i{Y}_i}\right)={\left({\bar{X}}^{\text{T}}\bar{X}\right)}^{-1}{\bar{X}}^{\text{T}}\bar{Y}.$ (6)

下面在模型(4)或(5)的基础上考虑岭估计的构造。类似于普通线性回归模型的估计，构造的如下的辅助函数：

$F_1\left(\beta \right)={\left(\bar{Y}-\bar{X}\beta \right)}^{\text{T}}\left(\bar{Y}-\bar{X}\beta \right)+k{\beta }^{\text{T}}\beta$ (7)

函数 $F\left(\beta \right)$ 关于 $\beta$ 求偏导数，并另导数等于0，可得

$\frac{\partial F_1\left(\beta \right)}{\partial \beta }=-2{\bar{X}}^{\text{T}}\bar{Y}+2{\bar{X}}^{\text{T}}\bar{X}\beta +2k\beta =0$ (8)

从而可得回归系数基于完整数据方法的岭估计为

${\hat{\beta }}_C\left(k\right)={\left({\bar{X}}^{\text{T}}\bar{X}+k{I}_q\right)}^{-1}{\bar{X}}^{\text{T}}\bar{Y}$ (9)

下面给出 ${\hat{\beta }}_c\left(k\right)$ 的渐近性质。

定理1. 如果第6节的假设条件成立， ${\hat{\beta }}_c\left(k\right)$ 是渐近正态的，有

$\sqrt{n}\left({\hat{\beta }}_c^r-\beta \right)\stackrel{d}{\to }N\left(0,{\sigma }^2{\Omega }^{-1}\right),$

其中 $\Omega =E\left(\delta X{X}^{\text{T}}\right)$。

3. 基于单点插补方法的约束估计及其性质

上一节所用的完整数据方法只是用了因变量和自变量有完整观测的数据，将因变量存在缺失的观测数据舍弃，显然造成了信息的损失。为了弥补这一问题，下面将基于单点插补方法构造模型系数的岭估计。基于上一节得到的最小二乘估计 ${\hat{\beta }}_C$，针对因变量缺失这一问题，构造如下新的因变量

$Y_i^{*}={\delta }_i{Y}_i+\left(1-{\delta }_i\right)X_i^{\text{T}}{\hat{\beta }}_C,$ (10)

显然，当 $Y_i$ 不存在缺失时， $Y_i^{*}=Y_i$，另一方面当 $Y_i$ 存在缺失时，相当于用插补值 $Y_i^{*}=X_i^{\text{T}}{\hat{\beta }}_C$ 代替其缺失的真实值。

基于上面的构造，得到了完整数据集 ${\left(Y_i^{*},X_i\right)}_{i=1}^n$，因此有如下的自变量和因变量都存在的线性模型

$Y_i^{*}=X_i^{\text{T}}\beta +e_i,i=1,2,\cdots ,n.$ (11)

其中 $Y_i^{*}={\delta }_i{Y}_i+\left(1-{\delta }_i\right)X_i^{\text{T}}{\hat{\beta }}_c$。

对于上述模型，利用最小二乘方法可以得到 $\beta$ 的单点插补估计

${\hat{\beta }}_I={\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}{X}_i{X}_i^{\text{T}}}\right)}^{-1}\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}{X}_i{Y}_i^{*}}\right).$ (12)

基于模型(11)，考虑其岭估计，构造如下的辅助函数

$F_2\left(\beta ,\lambda \right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(Y_i^{*}-X_i^{\text{T}}\beta \right)}^2}+k{\beta }^{\text{T}}\beta$

同样，函数 $F_2\left(\beta \right)$ 关于 $\beta$ 求偏导数，并另导数等于0，可得

$\frac{\partial F_2\left(\beta \right)}{\partial \beta }=-2X^{\text{T}}{Y}^{*}+2X^{\text{T}}X\beta +2k\beta =0$ (13)

$Y^{*}={\left(Y_1^{*},Y_2^{*},\cdots ,Y_n^{*}\right)}^{\text{T}}$，从而可得回归系数基于单点插补数据方法的岭估计为

${\hat{\beta }}_I\left(k\right)={\left(X^{\text{T}}X+k{I}_q\right)}^{-1}{X}^{\text{T}}{Y}^{*}$ (14)

下面给出 ${\hat{\beta }}_I\left(k\right)$ 的渐近性质。

定理2. 如果第6节的假设条件成立， ${\hat{\beta }}_I^r$ 是渐近正态的，满足

$\sqrt{n}\left({\hat{\beta }}_I\left(k\right)-\beta \right)\stackrel{d}{\to }N\left(0,{\sigma }^2{\Omega }^{-1}\right).$

从定理1和定理2不难看出，首先基于两种方法的估计量的渐近性质相同，且与不考虑岭估计的最小二乘估计的渐近性质一样。这些结论与杨徐佳等(2011)以及其他文献的结论一致。

4. 数值模拟

本节将通过数值模拟考察前面所提出估计方法的有效性。假设数据服从于如下线性回归模型

$y_i=x_{1i}{\beta }_1+x_{2i}{\beta }_2+{\epsilon }_i,i=1,2,\cdots ,n$ (15)

为了构造多重共线性，运用McDonald和Galerneau (1975) [8] 中生成具有多重共线性自变量的方法生成解释变量 $x_1,x_2$，具体为：

$x_{ij}={\left(1-{\rho }^2\right)}^{1/2}{\omega }_{ij}+\rho {\omega }_{i3},i=1,2,\cdots ,n;j=1,2$

其中 ${\omega }_{ij}$ 是独立的标准正态随机数， $\rho$ 是一个具体的数值以确保任何四个解释变量在理论上是相关的，分别取 $\rho =0.90,0.99,0.999$ 刻画不同程度的复共线性问题。

因变量的缺失的机制

$\Delta =p\left(\delta =1|x_1=x_{1i},x_2=x_{2i}\right)=0.8+0.2\left(\left|x_{1i}\right|+\left|x_{2i}\right|\right)$，当 $\left|x_{1i}\right|+\left|x_{2i}\right|\le 1$ 时；

否则等于0.9。模型误差 ${\epsilon }_i$ 服从如下的正态分布(N)和均匀分布(U)

(1) ${\epsilon }_i\sim N\left(0,{0.5}^2\right)$，(2) ${\epsilon }_i\sim U\left(-\frac{\sqrt{3}}{2},\frac{\sqrt{3}}{2}\right)$。

针对模型(15)，取 $\left({\beta }_1,{\beta }_2\right)$ 的真实值为 $\left(3,2\right)$，基于上面的缺失机制和误差分布，样本量n分别设置为50、100和150，重复500次，在每一种情况下分别计算 $\left({\beta }_1,{\beta }_2\right)$ 基于完整数据分析方法的估计(C)和岭估计(C-R)，基于单点插补方法的估计(I)和岭估计(I-R)。其中岭估计是基于R软件MASS package里的lm.ridge函数，其中k的选取使用GCV方法。以这些估计量的均方误差(EMSE)来衡量其表现，

$\text{EMSE}\left({\beta }^{*}\right)=\frac{1}{500}{\displaystyle \underset{k=1}{\overset{500}{\sum }}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{\left({\beta }_{kj}^{*}-{\beta }_j\right)}^2}}$

其中 ${\beta }_{kj}^{*}$ 是参数 ${\beta }_j$ 的第k次重复时的估计值，模拟结果见表1。

Table 1. EMSEs of the estimators

表1. 不同估计量的EMSEs

从模拟结果可以看出：1) 随着样本量的增大，这四类估计的均方误差值都在变小，与理论性质相一致。2) 误差分布对这四类估计的影响很小。3) 随着共线性程度的增加，岭估计优于其对应的最小二乘估计。4) 完整数据估计和单点插补估计相比，单点插补估计由于充分利用了数据的信息从而表现优于损失了信息的完整数据估计。

5. 总结

因变量缺失问题是使用线性回归模型进行实际问题分析时经常遇到的，目前的研究大都是基于最小二乘法或极大似然估计方法进行模型估计，很少有论文同时讨论自变量的多重共线性问题。本文就是针对这一问题，构造了因变量缺失下的线性回归模型的岭估计，并研究了所提估计量的渐近性质，通过数值模拟验证了所提估计的有效性。本文主要讨论了线性回归模型，论文的方法可以推广到半参数模型等其他类型的回归模型上。

6. 定理的证明

在给出定理的证明之前，我们先给出下面条件。

条件1： $\Omega =E\left(\delta X{X}^{\text{T}}\right)$ 为正定矩阵。

条件2： $E\left(\epsilon |X\right)=0$， $E\left({\left|\epsilon \right|}^3|X\right)<\infty$。

引理1. 如果前面的假设条件成立， ${\hat{\beta }}_c$ 和 ${\hat{\beta }}_I$ 都是渐进正态的，二者都满足

$\sqrt{n}\left(\hat{\beta }-\beta \right)\stackrel{d}{\to }N\left(0,{\sigma }^2{\Omega }^{-1}\right),$

其中 $\hat{\beta }$ 为 ${\hat{\beta }}_c$ 或 ${\hat{\beta }}_I$。

证明：该引理即为杨徐佳等(2011) [6] 中的定理1。

定理1证明：由 ${\hat{\beta }}_c\left(k\right)$ 的定义可得

$\begin{array}{c}{\hat{\beta }}_c\left(k\right)={\left({\bar{X}}^{\text{T}}\bar{X}+k{I}_q\right)}^{-1}{\bar{X}}^{\text{T}}\bar{Y}\\ ={\left(X^{\text{T}}\Delta X+k{I}_q\right)}^{-1}{X}^{\text{T}}\Delta Y\\ =\beta -k{\left(X^{\text{T}}\Delta X+k{I}_q\right)}^{-1}\beta +{\left(X^{\text{T}}\Delta X+k{I}_q\right)}^{-1}{X}^{\text{T}}\Delta \epsilon \end{array}$

基于条件1和2，可得到

$\frac{1}{n}\left(X^{\text{T}}\Delta X+k{I}_q\right)\stackrel{P}{\to }\Omega ,\frac{1}{\sqrt{n}}k\beta \stackrel{P}{\to }0,\frac{1}{\sqrt{n}}{X}^{\text{T}}\Delta \epsilon \stackrel{L}{\to }\Omega ,$

由上述结论，根据Slutsky定理可得

$\sqrt{n}\left({\hat{\beta }}_C\left(k\right)-\beta \right)\stackrel{d}{\to }N\left(0,{\sigma }^2{\Omega }^{-1}\right).$

定理2的证明和定理1类似，在此省略。

基金项目

中国劳动关系学院教育教学改革立项项目(JG1406)；2020年度安徽高等学校自然科学项目(KJ2020A1200)。