因变量缺失下线性回归模型的随机约束估计

Stochastic Restricted Estimation of Linear Regression Models with Missing Responses

摘要: 本文研究了线性回归模型在因变量存在缺失的同时回归系数附加有随机约束条件时的估计问题，基于完整数据方法和单点插补方法，给出了模型系数的两种约束估计，并研究了这两类估计量的渐近正态性。最后通过数值模拟验证了所提估计方法的有效性。

Abstract: This paper discusses estimation of linear regression models in the presence of multicollinearity and there are stochastic linear restrictions on the regression coefficients. Based on the complete-case method and single imputation technique, two stochastic restricted estimators are proposed. Asymptotically, properties of the proposed estimators are shown. Finally, some simulations are conducted to illustrate the proposed methods.

1. 引言

线性回归模型作为最为常用的统计分析之一，已经得到了深入的理论研究和广泛的应用。众所周知，对于线性回归模型，可以采用最小二乘估计方法或极大似然估计方法等。然而，这些估计方法一般只利用了样本总体假定和样本信息。在分析一些实际领域的问题时，常常可以通过该问题的专业背景以及前期研究等渠道获得一些额外信息，这些额外信息经常被转化为模型中回归系数的约束条件，比如某个回归系数按照实际问题来说应该是正值，或者两个或多个回归系数之间存在某种关系，详细介绍可参考Toutenburg (1982) [1]。那么，在模型的估计中如果能够将这些约束条件充分利用，可以提高估计量的有效性。对于一般线性回归模型，如果约束条件为针对系数的线性约束，那么此时采用约束最小二乘估计是比普通最小二乘法更有效。除了线性约束之外，还有随机线性约束以及不等式约束等情况，具体内容可参考Rao和Toutenburg (1999) [2]。

另一方面，实际数据分析中，数据的缺失是非常普遍的，关于缺失数据的详细介绍可参考著作Little和Rubin (1987) [3]。缺失数据下回归模型的研究得到了关注，相关研究可参考Cheng (1994) [4]、Chu & Cheng (1995) [5]、Wang & Rao (2002) [6]，Wang & Rao (2004) [7] 和Qin等(2009) [8] 等有关文献。

目前关于线性回归模型的约束估计，大都是基于无缺失的情形，关于缺失数据下约束估计的研究还很少。杨徐佳等(2011) [9] 构造了因变量缺失下线性回归模型的估计。安佰玲等(2013) [10] 在此基础上讨论了模型的精确约束估计问题。本文将在这两篇论文的基础上考虑因变量缺失下线性回归模型的随机约束估计。

考虑如下的线性回归模型

$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)

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

考虑如下的随机约束条件

$b=A\beta +\eta$ (3)

其中A是 $k\times p$ 维的已知矩阵，且 $rank\left(A\right)=k$，b是 $k\times 1$ 维的已知向量。 $\eta$ 为均值为0，协方差为 ${\sigma }^2\Omega$ 的随机向量，其中 $\Omega$ 为一已知正定矩阵。为了在线性回归模型的估计中考虑随机约束条件(3)，Durbin (1953) [11]，Theil和Goldberger (1961) [12] 以及Theil (1963) [13] 提出了混合估计方法。

本文重点研究线性回归模型(1)在因变量缺失机制(2)和随机约束条件(3)下的估计问题，从而将杨徐佳等(2011) [9] 和安佰玲等(2013) [10] 的结果推广到了随机约束情形。

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

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}}$。

基于最小二乘法，可得 $\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)

下面基于文 [11] [12] [13] 中的混合估计方法来估计随机约束条件(3)下的线性回归模型(5)。构造的如下的辅助函数：

$F\left(\beta \right)={\left(\bar{Y}-\bar{X}\beta \right)}^{\text{T}}\left(\bar{Y}-\bar{X}\beta \right)+{\left(b-A\beta \right)}^{\text{T}}{\Omega }^{-1}\left(b-A\beta \right)$ (7)

用辅助函数对于 $\beta$ 求导数并令其为0，可得：

$\frac{\partial F\left(\beta \right)}{\partial \beta }=-2{\bar{X}}^{\text{T}}\left(\bar{Y}-\bar{X}\beta \right)-2A^{\text{T}}{\Omega }^{-1}{\left(b-A\beta \right)}^{\text{T}}=0.$ (8)

简单整理可得 $\beta$ 基于完整数据方法的随机约束估计

${\hat{\beta }}_C^{SR}={\left({\bar{X}}^{\text{T}}\bar{X}+A^{\text{T}}{\Omega }^{-1}A\right)}^{-1}\left({\bar{X}}^{\text{T}}\bar{Y}+A^{\text{T}}{\Omega }^{-1}b\right).$ (9)

根据文 [9] 中的定理A.18，可得:

${\left({\bar{X}}^{\text{T}}\bar{X}+A^{\text{T}}{\Omega }^{-1}{A}^{\text{T}}\right)}^{-1}={\left({\bar{X}}^{\text{T}}\bar{X}\right)}^{-1}-{\left({\bar{X}}^{\text{T}}\bar{X}\right)}^{-1}{A}^{\text{T}}{\left[\Omega +A{\left({\bar{X}}^{\text{T}}\bar{X}\right)}^{-1}{A}^{\text{T}}\right]}^{-1}A{\left({\bar{X}}^{\text{T}}\bar{X}\right)}^{-1}$ (10)

则由(6)、(9)和(10)， ${\hat{\beta }}_C^{SR}$ 可以等价表示为如下的形式

${\hat{\beta }}_C^{RS}={\hat{\beta }}_C-{\left({\bar{X}}^{\text{T}}\bar{X}\right)}^{-1}{A}^{\text{T}}{\left[\Omega +A{\left({\bar{X}}^{\text{T}}\bar{X}\right)}^{-1}{A}^{\text{T}}\right]}^{-1}\left(A{\hat{\beta }}_C-b\right)$ (11)

下面给出关于 ${\hat{\beta }}_C^{RS}$ 的渐近性质。

定理1. 如果第5节的假设条件成立， ${\hat{\beta }}_C^{RS}$ 是渐近正态的，有

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

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

显然，我们所构造的因变量缺失下的随机约束估计 ${\hat{\beta }}_C^{RS}$ 的渐近性质与杨徐佳等(2011) [9] 中构造的精确约束最小二乘估计的渐近性质相同。

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

为了弥补完整数据方法丢弃数据从而损失信息的不足，下面基于单点插补方法构造模型系数的随机约束估计。基于上一节得到的最小二乘估计 ${\hat{\beta }}_C$，针对因变量存在缺失这一情况，定义如下的新的因变量

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

显然，当 $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.$ (13)

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

对模型(13)使用最小二乘方法，得到 $\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).$ (14)

考虑随机约束条件(3)，构造如下的辅助函数

$F^{*}\left(\beta ,\lambda \right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(Y_i^{*}-X_i^{\text{T}}\beta \right)}^2}+{\left(b-A\beta \right)}^{\text{T}}{\Omega }^{-1}\left(b-A\beta \right).$ (15)

同第2节类似，基于上面的辅助函数，可得 $\beta$ 单点插补的随机约束估计为

${\hat{\beta }}_I^{SR}={\left(X^{\text{T}}X+A^{\text{T}}{\Omega }^{-1}A\right)}^{-1}\left(X^{\text{T}}{Y}^{*}+A^{\text{T}}{\Omega }^{-1}b\right)$ (16)

同样，该估计可以等价表达为

${\hat{\beta }}_I^{RS}={\hat{\beta }}_I-{\left(X^{\text{T}}X\right)}^{-1}{A}^{\text{T}}{\left[\Omega +A{\left(X^{\text{T}}X\right)}^{-1}{A}^{\text{T}}\right]}^{-1}\left(A{\hat{\beta }}_I-b\right).$ (17)

下面给出 ${\hat{\beta }}_I^{RS}$ 的渐近性质。

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

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

需要注意的是 ${\hat{\beta }}_I^{RS}$ 的渐近性质也与杨徐佳等(2013) [9] 中构造的完整数据最小二乘估计的渐近性质相同。

4. 数值模拟

本节我们将通过数值模拟考察前面所提出估计方法的有效性。假设数据产生于如下的模型

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

其中 $x_{1i}~N\left(1,1\right)$， $x_{2i}~U\left(0,1\right)$， $\left({\beta }_1,{\beta }_2\right)$ 的真实值为 $\left(3,2\right)$。为了考察模型误差的分布对估计和检验结果影响，考虑如下两种情况(1) ${\epsilon }_i~N\left(0,{0.5}^2\right)$ ；(2) ${\epsilon }_i~U\left(-\sqrt{3}/2,\sqrt{3}/2\right)$，表1中分别用N和U表示。缺失的机制采用安佰玲等(2013) [10] 中的设置，当 $\left|x_{1i}-1\right|+\left|x_{2i}-0.5\right|\le 1$ 时 $\Delta =p\left(\delta =1|x_1=x_{1i},x_2=x_{2i}\right)=0.8+0.2\left(\left|x_{1i}-1\right|+\left|x_{2i}-0.5\right|\right)$ ；其余情况等于0.9。

随机约束设定为 ${\beta }_1+{\beta }_2=5+\eta$， $E\eta =0$， $Var\left(\eta \right)=0.49$，即 $A=\left(1,1\right)$， $b=5$， $\Omega =0.49$。

针对模型(18)，基于上面的各种缺失机制和误差分布，我们取样本量n分别为50、100和150，在每一种设定下分别求取 $\left({\beta }_1,{\beta }_2\right)$ 的基于完整数据分析方法的估计(表1中用C表示，下同)和随机约束估计(C-SR)，基于单点插补方法的估计(I)和随机约束估计(I-SR)。每种情况重复计算500次，以这些估计量的均方误差(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. 不同估计量的MSE

从模拟结果可以看出：(1) 这四类估计量的均方误差都随着样本量的增加而减小，此外，估计值对于误差分布的改变几乎没有变化。(2) 同样设置下，单点插补方法的均方误差一般小于完整估计方法。(3) 同样设置下，考虑了约束条件的随机约束估计的表现优于没有考虑约束条件时的估计。

5. 定理的证明

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

条件1： $\sum =E\left({\delta }_1{X}_1{X}_1^{\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) [9] 中的定理1。

引理2. 如果前面的假设条件成立，有

$\frac{1}{n}{\bar{X}}^{\text{T}}\bar{X}\stackrel{P}{\to }\Sigma ,\frac{1}{n}{X}^{\text{T}}X\stackrel{P}{\to }E{X}_1{X}_1^{\text{T}}$

定理1的证明：定义 $S_{SR}={\bar{X}}^{\text{T}}\bar{X}+A^{\text{T}}{\Omega }^{-1}A$， $\Delta =\text{diag}\left({\delta }_1,{\delta }_2,\cdots ,{\delta }_n\right)$，则由 ${\hat{\beta }}_C^{RS}$ 的定义可得

$\begin{array}{c}{\hat{\beta }}_C^{RS}=S_{SR}^{-1}\left({\bar{X}}^{\text{T}}\bar{Y}+A^{\text{T}}{\Omega }^{-1}b\right)=S_{SR}^{-1}\left(X^{\text{T}}\Delta Y+A^{\text{T}}{\Omega }^{-1}b\right)\\ =S_{SR}^{-1}\left[X^{\text{T}}\Delta X\beta +X^{\text{T}}\Delta \epsilon +A^{\text{T}}{\Omega }^{-1}\left(A\beta +\eta \right)\right]\\ =S_{SR}^{-1}\left[\left(X^{\text{T}}\Delta X+A^{\text{T}}{\Omega }^{-1}A\right)\beta +X^{\text{T}}\Delta \epsilon +A^{\text{T}}{\Omega }^{-1}\eta \right]\\ =\beta +S_{SR}^{-1}{X}^{\text{T}}\Delta \epsilon +S_{SR}^{-1}{A}^{\text{T}}{\Omega }^{-1}\eta \end{array}$

从而有

$\sqrt{n}\left({\hat{\beta }}_C^{RS}-\beta \right)={\left(\frac{S_{SR}}{n}\right)}^{-1}\frac{1}{\sqrt{n}}{X}^{\text{T}}\Delta \epsilon +{\left(\frac{S_{SR}}{n}\right)}^{-1}\frac{1}{\sqrt{n}}{A}^{\text{T}}{\Omega }^{-1}\eta$

由引理2，结合 $\frac{1}{n}{A}^{\text{T}}{\Omega }^{-1}A=o_p\left(1\right)$，可得：

$\frac{S_{SR}}{n}=\frac{1}{n}\left({\bar{X}}^{\text{T}}\bar{X}+A^{\text{T}}{\Omega }^{-1}A\right)\stackrel{p}{\to }\sum ,\frac{1}{\sqrt{n}}{X}^{\text{T}}\Delta \epsilon \stackrel{D}{\to }N\left(0,{\sigma }^2\sum \right)$

另一方面由 $E{A}^{\text{T}}{\Omega }^{-1}\eta =0$ 和 $Cov\left(A^{\text{T}}{\Omega }^{-1}\eta \right)=A^{\text{T}}{\Omega }^{-1}A$，可得

$\frac{1}{\sqrt{n}}{A}^{\text{T}}{\Omega }^{-1}\eta =o_p(\; 1\; )$

结合上面的结论，由Slutsky定理，可得

$\sqrt{n}\left({\hat{\beta }}_C^{RS}-\hat{\beta }\right)\stackrel{D}{\to }N\left(0,{\sigma }^2{\sum }^{-1}\right)$

定理2的证明：由 $Y_i^{*}={\delta }_i{Y}_i+\left(1-{\delta }_i\right)X_i^{\text{T}}{\hat{\beta }}_C$，可得

$\begin{array}{c}{Y}^{*}=\Delta Y+\left(I_n-\Delta \right)X{\hat{\beta }}_C\\ =\Delta Y+\left(I_n-\Delta \right)X{\left(X^{\text{T}}\Delta X\right)}^{-1}{X}^{\text{T}}\Delta Y\\ =\Delta Y+X{\left(X^{\text{T}}\Delta X\right)}^{-1}{X}^{\text{T}}\Delta Y-\Delta X{\left(X^{\text{T}}\Delta X\right)}^{-1}{X}^{\text{T}}\Delta Y\end{array}$

从而有

$\begin{array}{c}{\beta }_I^{RS}={\left(X^{\text{T}}X+A^{\text{T}}{\Omega }^{-1}A\right)}^{-1}\left(X^{\text{T}}{Y}^{*}+A^{\text{T}}{\Omega }^{-1}b\right)\\ ={\left(X^{\text{T}}X+A^{\text{T}}{\Omega }^{-1}A\right)}^{-1}\left[X^{\text{T}}\left(\Delta Y+X{\left(X^{\text{T}}\Delta X\right)}^{-1}{X}^{\text{T}}\Delta Y-\Delta X{\left(X^{\text{T}}\Delta X\right)}^{-1}{X}^{\text{T}}\Delta Y\right)+A^{\text{T}}{\Omega }^{-1}b\right]\\ ={\left(X^{\text{T}}X+A^{\text{T}}{\Omega }^{-1}A\right)}^{-1}\left[X^{\text{T}}X{\left(X^{\text{T}}\Delta X\right)}^{-1}{X}^{\text{T}}\Delta Y+A^{\text{T}}{\Omega }^{-1}b\right]\\ ={\left(X^{\text{T}}X+A^{\text{T}}{\Omega }^{-1}A\right)}^{-1}\left[X^{\text{T}}X{\left(X^{\text{T}}\Delta X\right)}^{-1}{X}^{\text{T}}\Delta \left(X\beta +\epsilon \right)+A^{\text{T}}{\Omega }^{-1}\left(A\beta +\eta \right)\right]\\ =\beta +{\left(X^{\text{T}}X+A^{\text{T}}{\Omega }^{-1}A\right)}^{-1}\left[X^{\text{T}}X{\left(X^{\text{T}}\Delta X\right)}^{-1}{X}^{\text{T}}\Delta \epsilon +A^{\text{T}}{\Omega }^{-1}\eta \right]\end{array}$

由引理2以及 $\frac{1}{n}{A}^{\text{T}}{\Omega }^{-1}A\stackrel{P}{\to }0$ 可得

$\begin{array}{c}\sqrt{n}\left({\hat{\beta }}_{PM}-\beta \right)={\left(\frac{X^{\text{T}}X+A^{\text{T}}{\Omega }^{-1}A}{n}\right)}^{-1}{X}^{\text{T}}X{\left(X^{\text{T}}\Delta X\right)}^{-1}\frac{1}{\sqrt{n}}{X}^{\text{T}}\Delta \epsilon +{\left(\frac{X^{\text{T}}X+A^{\text{T}}{\Omega }^{-1}A}{n}\right)}^{-1}\frac{1}{\sqrt{n}}{A}^{\text{T}}{\Omega }^{-1}\eta \\ ={\left(X^{\text{T}}\Delta X\right)}^{-1}\frac{1}{\sqrt{n}}{X}^{\text{T}}\Delta \epsilon +o_p(\; 1\; )\end{array}$

从而可得

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

基金项目

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