响应变量缺失下可伸缩的模型平均估计

Scalable Model Averaging Estimation with Missing Responses

摘要: 当响应变量缺失时，本文研究了线性模型的可伸缩的模型平均估计问题。基于逆概率加权方法和奇异值分解对原模型进行转换，所构建的模型平均估计最多只需要选择p个候选模型的权重。利用Jackknife准则选择权重，且权重之和不必限制为1。在一定的条件下证明了所提方法在实现最小二次损失意义下是渐近最优的，并通过数值模拟进一步验证了所提方法在预测效率和计算成本方面的优良性。

Abstract: In this paper, a scalable model averaging method is developed for the linear regression models with response missing. By using inverse probability weighted method and the singular value decomposition to transform the original models, this method enables us to find the optimal weights by considering at most p candidate models. The weights can be selected by Jackknife criterion, and the sum of weights is not necessarily restricted to one. Under some mild conditions, it is shown that the proposed method is asymptotically optimal in the sense of achieving the lowest possible squared error. Some simulation studies are conducted to illustrate the superiority of the proposed method in terms of both predictive efficiency and computational cost.

1. 引言

模型平均是当代统计分析中用于处理模型不确定性和提高预测精度的一种流行且有效的方法。其中，针对频率模型平均(Frequency model Averaging, FMA)的方法和理论研究获得了广泛的关注。Hansen [1] 开创性地提出了Mallows模型平均(MMA)方法，该方法利用极小化Mallows准则来选择嵌套模型的权重，并证明了该方法的渐近最优性。Wan等 [2] 将MMA方法拓展到连续权重和非嵌套模型的情形。当随机误差项为异方差时，Hansen和Racine [3] ，Liu和Okui [4] 分别提出了JMA (Jackknife Model Averaging)和HRCp (Heteroscedasticity Robust Cp)模型平均方法，Zhang等 [5] 进一步将JMA方法拓展到相关数据的情形。基于这些重要工作，很多学者将这些方法扩展到不同的数据类型和不同的模型中。

其中，缺失数据是在很多实际问题研究中经常会遇到的数据类型，因此缺失数据下模型平均估计也得到很多的关注和研究。当响应变量缺失时，Wei和Wang [6] ，Wei等 [7] 提出了线性模型的模型平均估计，并分别利用JMA和HRCp准则选择权重。Zeng等 [8] 扩展到变系数偏线性模型的情形。但是，当协变量的个数p较大时，上述方法可选的候选模型的个数可高达 $2^p$ ，这给每个候选模型的估计和权重的选择都带来了极大的困难和挑战。

为了提高估计效率和计算速度，Zhu等 [9] 提出了可伸缩的模型平均估计，该方法利用奇异值分解对模型进行转换，所构建的模型平均估计只需要估计p个权重，且其最小二次损失在渐近意义下不差于传统的模型平均估计。因此，受 [9] 的启发，本文考虑了线性模型在响应变量缺失下的可伸缩的模型平均估计。在响应变量随机缺失(missing at random, MAR)的前提下，我们基于逆概率加权方法和奇异值分解对模型进行转换，使得转换后的模型只需估计p个候选模型，并通过极小化Jackknife准则选择p个权重，从而有效地减轻了计算压力，提高了模型平均的预测效率。而且，在一定的条件下证明了所提模型平均估计方法的渐近最优性，并通过数值模拟实验进一步验证了所提方法的优良性。

2. 模型及估计

2.1. 模型假设和估计方法

考虑线性回归模型：

$y_i={\mu }_i+e_i=x_i^{\text{T}}\theta +e_i={\displaystyle {\sum }_{j=1}^p{\theta }_j{x}_{i\text{}j}}+e_i,i=1,\cdots ,n,$ (1)

其中 $y_i$ 是响应变量， $x_i={\left(x_{i1},\cdots ,x_{i\text{}p}\right)}^{\text{T}}$ 是预测变量， $\theta ={\left({\theta }_1,\cdots ,{\theta }_p\right)}^{\text{T}}$ 是回归系数， $e_1,\cdots ,e_n$ 是相互独立且为异方差的随机误差，满足 $E\left(e_i|x_i\right)=0$ ， $E\left(e_i^{\text{}\text{}2}|x_i\right)={\sigma }_i^2$ 。我们考虑 $p<n$ 的情形。这里假定预测变量 $x_i$ 的值可以被完全观测，响应变量 $y_i$ 是随机缺失的，令 ${\delta }_i$ 为缺失示性变量，即当 $y_i$ 缺失时 ${\delta }_i=0$ ，否则 ${\delta }_i=1$ ，且缺失机制满足MAR条件，即：

$P\left({\delta }_i=1|y_i,x_i\right)=P\left({\delta }_i=1|x_i\right)=\pi \left(x_i\right)$ 。

假设 $\left\{\left(y_i,x_i,{\delta }_i\right),i=1,\cdots ,n\right\}$ 是独立同分布的。当选择概率函数 $\pi \left(x_i\right)$ 已知时，记 $z_{\pi ,i}={\delta }_i{y}_i/\pi \left(x_i\right)$ ，在MAR假设下，易得 $E\left(z_{\pi ,i}|x_i\right)={\mu }_i=x_i^{\text{T}}\theta$ ，因此有

$z_{\pi ,i}=x_i^{\text{T}}\theta +e_{\pi ,i}={\displaystyle {\sum }_{j=1}^p{\theta }_j{x}_{i\text{}j}}+e_{\pi ,i},i=1,\cdots ,n,$ (2)

其中 $e_{\pi ,i}={\delta }_i\left({\mu }_i+e_i\right)/\pi \left(x_i\right)-{\mu }_i$ ， $E\left(e_{\pi ,i}|x_i\right)=0$ ， $E\left(e_{\pi ,i}^2|x_i\right)={\sigma }_{\pi ,i}^2=\left({\sigma }_i^2+{\mu }_{\text{}i}^2\right)/\pi \left(x_i\right)-{\mu }_{\text{}i}^2$ 。模型(2)可表示为矩阵形式：

$z_{\pi }=\mu +e_{\pi }=X\theta +e_{\pi },$ (3)

其中 $z_{\pi }={\left(z_{\pi ,1},\cdots ,z_{\pi ,n}\right)}^{\text{T}}$ ， $\mu ={\left({\mu }_1,\cdots ,{\mu }_n\right)}^{\text{T}}$ ， $e_{\pi }={\left(e_{\pi }{}_{,1},\cdots ,e_{\pi ,n}\right)}^{\text{T}}$ ， $X={\left(x_1,\cdots ,x_n\right)}^{\text{T}}$ 。

对数据生成过程(1)，考虑M个候选模型，相应地，对(3)式，也有M个候选模型：

$z_{\pi }=X_m{\theta }_m+e_{\pi },m=1,\cdots ,M,$ (4)

其中 $X_m$ 第m个模型对应的 $n\times k_m$ 阶预测矩阵， ${\theta }_m$ 是 $k_m$ 维回归系数，记 $S_m$ 为该模型包含的预测变量的指标集，则 $S_m$ 的维度为 $k_m$ 。

由于 $\pi \left(x_i\right)$ 一般都是未知的，故这里对 $\pi \left(x_i\right)$ 假设一个参数模型 $\bar{\pi }\left(x_i;\alpha \right)$ ，其中 $\bar{\pi }\left(\cdot ;\alpha \right)$ 是已知的函数， $\alpha$ 是未知的参数向量。记 $\hat{\alpha }$ 为 $\alpha$ 的一个估计量，如极大似然估计量，记 $\hat{\pi }\left(x_i\right)=\bar{\pi }\left(x_i;\hat{\alpha }\right)$ 。因此，当需要估计 $\pi \left(x_i\right)$ 时，用 $\hat{\pi }\left(x_i\right)$ 替代上述模型(2)~(4)中的 $\pi \left(x_i\right)$ ，相应地，记号的下标 $\pi$ 替换成 $\hat{\pi }$ ，如： $z_{\hat{\pi },i}={\delta }_i{y}_i/\hat{\pi }\left(x_i\right)$ ， $z_{\hat{\pi }}={\left(z_{\hat{\pi },1},\cdots ,z_{\hat{\pi },n}\right)}^{\text{T}}$ ， $e_{\hat{\pi }}={\left(e_{\hat{\pi }}{}_{,1},\cdots ,e_{\hat{\pi },n}\right)}^{\text{T}}$ 等。

对第m个模型，利用最小二乘(OLS)估计可得到 ${\theta }_m,\mu$ 的估计 ${\tilde{\theta }}_{\hat{\pi },m}={\left(X_m^{\text{T}}{X}_m\right)}^{-1}{X}_m^{\text{T}}{z}_{\hat{\pi }}$ ， ${\tilde{\mu }}_{\hat{\pi },m}=P_m{z}_{\hat{\pi }}$ ，其中 $P_m=X_m{\left(X_m^{\text{T}}{X}_m\right)}^{-1}{X}_m^{\text{T}}$ ， $m=1,\cdots ,M$ 。令 $\tilde{w}={\left({\tilde{w}}_1,\cdots ,{\tilde{w}}_M\right)}^{\text{T}}$ 为集合 ${\tilde{W}}^M=\left\{\tilde{w}\in {\left[0,1\right]}^M:{\displaystyle {\sum }_{m=1}^M{\tilde{w}}_m}=1\right\}$ 中的权重向量，可得 $\mu$ 的一个模型平均估计 ${\tilde{\mu }}_{\hat{\pi }}\left(\tilde{w}\right)={\displaystyle {\sum }_{m=1}^M{\tilde{w}}_m{\tilde{\mu }}_{\hat{\pi },m}}=P\left(\tilde{w}\right)z_{\hat{\pi }}$ ，其中 $P\left(\tilde{w}\right)={\displaystyle {\sum }_{m=1}^M{\tilde{w}}_m{P}_m}$ 。进一步，利用Jackknife准则 [3] 来选择权重，令

${\tilde{J}}_{\hat{\pi }}\left(\tilde{w}\right)={\Vert z_{\hat{\pi }}-\tilde{P}\left(\tilde{w}\right)z_{\hat{\pi }}\Vert }^2$ ，

其中 $\tilde{P}\left(\tilde{w}\right)={\displaystyle {\sum }_{m=1}^M{\tilde{w}}_m{\tilde{P}}_m}$ ， ${\tilde{P}}_m=D_m\left(P_m-I_n\right)+I_n$ ， $I_n$ 为n阶单位矩阵， $D_m$ 为对角元素为 ${\left(1-p_{m,ii}\right)}^{-1}$ 的n阶对角阵， $p_{m,ii}$ 为 $P_m$ 的第i个对角线元素。通过极小化Jackknife准则得到权重向量 ${\tilde{w}}_{\hat{\pi }}=\arg {\min }_{\tilde{w}\in {\tilde{W}}^M}{\tilde{J}}_{\hat{\pi }}\left(\tilde{w}\right)$ ，进一步得到 $\mu$ 的响应变量缺失下Jackknife Model Averaging (M-JMA)估计 ${\tilde{\mu }}_{\hat{\pi }}\left({\tilde{w}}_{\hat{\pi }}\right)$ 。由文献 [3] [6] ，在一定的条件下，该M-JMA估计在实现最小二次损失意义下是渐近最优的，因此，该方法得到了广泛的应用，且有很好的预测性能。

但是，当协变量的个数p较大时，可选的候选模型的个数M可高达 $2^p$ ，这给每个模型的估计和权重的选择都带了极大的困难和挑战，从而降低了估计效率和计算速度。因此，我们需要尝试新的方法去改进上述模型平均估计。参考文献 [9] ，本文提出了响应变量缺失下可伸缩的模型平均估计，基于逆概率加权和奇异值分解(SVD)将模型进行转换，转换后所构建的模型平均估计只需估计p个一元回归模型和选择p个权重，从而有效地减轻了计算压力，提高了预测精度。

假设 $\text{rank}\left(X\right)=p$ ，由奇异值分解(SVD)， $X=UD{V}^{\text{T}}$ ，使得 $U^{\text{T}}U=V^{\text{T}}V=I_p$ ，其中 $U$ 是 $n\times p$ 阶的列正交矩阵， $V$ 是p阶正交矩阵， $D$ 是对角线为非负实数的p阶对角矩阵。相应地， $X_m$ 可分解为 $X_m=U_m{D}_m{V}_m^{\text{T}}$ ，则候选模型可表示为

$z_{\hat{\pi }}=U_m{\beta }_m+e_{\hat{\pi }},m=1,\cdots ,M,$ (5)

其中 ${\beta }_m=D_m{V}_m^{\text{T}}{\theta }_m$ ，且 ${\beta }_m$ 的OLS估计为

${\hat{\beta }}_{\hat{\pi },m}={\left(U_m^{\text{T}}{U}_m\right)}^{-1}{U}_m^{\text{T}}{z}_{\hat{\pi }}=U_m^{\text{T}}{z}_{\hat{\pi }}={\left({\hat{\beta }}_{\hat{\pi },j_1},\cdots ,{\hat{\beta }}_{\hat{\pi },j_{k_m}}\right)}^{\text{T}}$ ，

其中 $j_1,\cdots ,j_{k_m}$ 是 $S_m$ 中的元素， ${\hat{\beta }}_{\hat{\pi },j}={\left(u_{\left(j\right)}^{\text{T}}{u}_{\left(j\right)}\right)}^{-1}{u}_{\left(j\right)}^{\text{T}}{z}_{\hat{\pi }}=u_{\left(j\right)}^{\text{T}}{z}_{\hat{\pi }}$ 是一元回归模型 $z_{\hat{\pi }}=u_{\left(j\right)}{\beta }_j+e_{\hat{\pi }}$ 的OLS估计， $u_{\left(j\right)}$ 是 $U$ 的第j列向量， $1\le j\le p$ ，即 ${\beta }_m$ 可以通过 $z_{\hat{\pi }}$ 关于 $U_m$ 的每一列做一元线性回归进行估计。从而， $\mu$ 可估计为 ${\hat{\mu }}_{\hat{\pi },m}=U_m\text{}{\hat{\beta }}_{\hat{\pi },m}={\displaystyle {\sum }_{j\in S_m}{u}_{\left(j\right)}{\hat{\beta }}_{\hat{\pi },j}}$ ，故 $\mu$ 的模型平均估计为

${\hat{\mu }}_{\hat{\pi }}\left(\tilde{w}\right)={\displaystyle {\sum }_{m=1}^M{\tilde{w}}_m{\hat{\mu }}_{\hat{\pi },m}}={\displaystyle {\sum }_{m=1}^M{\tilde{w}}_m\left({\displaystyle {\sum }_{j\in S_m}{u}_{\left(j\right)}{\hat{\beta }}_{\hat{\pi },j}}\right)}={\displaystyle {\sum }_{j=1}^p\left({\displaystyle {\sum }_{m=1}^M{\tilde{w}}_m{I}_{S_m}\left(j\right)}\right)u_{\left(j\right)}{\hat{\beta }}_{\hat{\pi },j}}$ ， (6)

其中 $I_{S_m}\left(j\right)$ 为示性函数，即当 $j\in S_m$ 时 $I_{S_m}\left(j\right)=1$ ，否则为0。若定义 $w_j={\displaystyle {\sum }_{m=1}^M{\tilde{w}}_m{I}_{S_m}\left(j\right)}\in \left[0,1\right]$ 为权重，则式(6)可看成为p个一元回归模型的预测平均。

令 $w={\left(w_1,\cdots ,w_p\right)}^{\text{T}}$ 为集合 $W=\left\{w_j\in \left[0,1\right]:1\le j\le p\right\}$ 中的权重向量，构造 $\mu$ 的模型平均估计量

${\hat{\mu }}_{\hat{\pi }}\left(w\right)={\displaystyle {\sum }_{j=1}^p{w}_j{\hat{\mu }}_{\left(j\right)}}={\displaystyle {\sum }_{j=1}^p{w}_j{u}_{\left(j\right)}{\hat{\beta }}_{\hat{\pi },j}}=H\left(w\right)z_{\hat{\pi }}$ ，

其中 $w\in W$ ， ${\hat{\mu }}_{\left(j\right)}=u_{\left(j\right)}{\hat{\beta }}_{\hat{\pi },j}$ ， $H\left(w\right)={\displaystyle {\sum }_{j=1}^p{w}_j{H}_j}$ ， $H_j=u_{\left(j\right)}{u}_{\left(j\right)}^{\text{T}}$ ，注意权重之和 ${\displaystyle {\sum }_{j=1}^p{w}_j}$ 无须等于1。进一步，利用Jackknife准则 [3] 来选择权重 $w$ ，令

$J_{\hat{\pi }}\left(w\right)={\Vert z_{\hat{\pi }}-\hat{H}\left(w\right)z_{\hat{\pi }}\Vert }^2$ ，

其中 $\hat{H}\left(w\right)={\displaystyle {\sum }_{j=1}^p{w}_j{\hat{H}}_j}$ ， ${\hat{H}}_j={\Lambda }_j\left(H_j-I_n\right)+I_n$ ， $h_{j,ii}$ 为 ${{\rm H}}_j$ 的第i个对角线元素， ${\Lambda }_j$ 为对角元素为 ${\left(1-h_{j,ii}\right)}^{-1}$ 的n阶对角阵。通过极小化Jackknife准则得到权重向量 ${\hat{w}}_{\hat{\pi }}=\arg {\min }_{w\in W}{J}_{\hat{\pi }}\left(w\right)$ ，进一步得到 $\mu$ 的响应变量缺失下Scalable Jackknife Model Averaging (M-SJMA)估计 ${\hat{\mu }}_{\hat{\pi }}\left({\hat{w}}_{\hat{\pi }}\right)$ 。该估计量一个直接的好处是，所选权重个数完成了从到 $2^p$ 到p的大幅缩小，极大地降低计算成本，减少计算中的数值误差，从而提高模型平均的预测效率，这在第3节的数值模拟结果中得到了很好的验证。

2.2. 渐近性质

本小节将讨论所提出的M-SJMA估计的渐近性质。为了说明的方便，先引入一些记号和定义。记 ${\Theta }_{\alpha }$ 为 $\alpha$ 的参数空间， $l\left(\alpha \right)=E\left[\delta \log \bar{\pi }\left(x;\alpha \right)+\left(1-\delta \right)\log \left(1-\bar{\pi }\left(x;\alpha \right)\right)\right]$ 。定义二次损失函数 ${\tilde{L}}_{\pi }\left(\tilde{w}\right)={\Vert {\tilde{\mu }}_{\pi }\left(\tilde{w}\right)-\mu \Vert }^2$ ， $\tilde{w}\in {\tilde{W}}^M$ ， $L_{\pi }\left(w\right)={\Vert {\hat{\mu }}_{\pi }\left(w\right)-\mu \Vert }^2$ ， $w\in W$ ，相应的风险函数 $R_{\pi }\left(w\right)=E\left[L_{\pi }\left(w\right)|X\right]$ ，其中 ${\hat{\mu }}_{\pi }\left(w\right)=H\left(w\right)z_{\pi }$ 。记 ${\xi }_{\pi }={\inf }_{w\in W}{R}_{\pi }\left(w\right)$ 。

接下来，给出证明渐近性质所需要的条件。需要指出的是，所有的大样本条件和结论均是在 $n\to \infty$ 的意义上的。

(C1) ${\Theta }_{\alpha }$ 是有界闭集， $l\left(\alpha \right)$ 在 ${\Theta }_{\alpha }$ 中关于 $\alpha$ 有唯一的极大值点 ${\alpha }_0$ ，且 ${\alpha }_0$ 为 ${\Theta }_{\alpha }$ 的内点。 $\bar{\pi }\left(x;\alpha \right)$ 关于 $\alpha$ 二阶连续可微，且 $E\left\{{\left[\partial \bar{\pi }\left(x;\alpha \right)/\partial \alpha \cdot \partial \bar{\pi }\left(x;\alpha \right)/\partial {\alpha }^{\text{T}}\right]|}_{\alpha ={\alpha }_0}\right\}$ 为正定矩阵。 ${\inf }_{\alpha \in {\Theta }_{\alpha }}\bar{\pi }\left(x;\alpha \right)>0$ 。

(C2)存在常数 $c_1,c_2$ ，使得 ${\max }_{1\le i\le n}E\left(e_i^4\right)\le c_1$ 和 ${\max }_{1\le i\le n}\left|{\mu }_i\right|\le c_2$ 几乎处处成立。

(C3) ${\xi }_{\pi }^{-1}p\stackrel{a.s.}{\to }0$ 。

(C4)存在常数 $c_3$ ，使得 ${\max }_{i,j}{u}_{ij}^2\le c_3{n}^{-1}$ ，其中 $u_{ij}$ 为矩阵 $U$ 中的元素。

这里，条件(C1)是极大似然估计量 $\hat{\alpha }$ 具有相合性和渐近正态性所需要的条件。条件(C2)~(C4)是模型平均估计中常用的条件，比如文献 [5] [9] 中都有类似的条件。

定理1：在条件(C1)~(C3)下，当 $n\to \infty$ 时，有 $\frac{{\inf }_{w\in W}{L}_{\hat{\pi }}\left(w\right)}{{\inf }_{\tilde{w}\in {\tilde{W}}^{2^p}}{\tilde{L}}_{\hat{\pi }}\left(\tilde{w}\right)}\stackrel{P}{\to }1$ ，其中 ${\tilde{W}}^{2^p}$ 是当 $M=2^p$ 时的权重向量集合 ${\tilde{W}}^M$ 。

定理1表明，可伸缩的模型平均方法与传统的 $2^p$ 个权重的模型平均方法相比，它们的最小二次损失是渐近相等的。

定理2：在条件(C1)~(C4)下，当 $n\to \infty$ 时，有 $\frac{L_{\hat{\pi }}\left({\hat{w}}_{\hat{\pi }}\right)}{{\inf }_{w\in W}{L}_{\hat{\pi }}\left(w\right)}\stackrel{P}{\to }1$ 。

定理3：在条件(C1)~(C4)下，当 $n\to \infty$ 时，有 $\frac{L_{\hat{\pi }}\left({\hat{w}}_{\hat{\pi }}\right)}{{\inf }_{\tilde{w}\in {\tilde{W}}^{2^p}}{\tilde{L}}_{\hat{\pi }}\left(\tilde{w}\right)}\stackrel{P}{\to }1$ 。

根据定理1和定理2的结论，定理3显然成立。该定理表明被选出的p个权重 ${\hat{w}}_{\hat{\pi }}$ 所对应的二次损失渐近等价于不可获得的 $2^p$ 个最优权重所对应的二次损失，这说明所提出M-SJMA估计 ${\hat{\mu }}_{\hat{\pi }}\left({\hat{w}}_{\hat{\pi }}\right)$ 在实现最小二次损失意义下是渐近最优的。

3. 数值模拟

本节将通过数值模拟来检验所提M-SJMA方法的有限样本性质。为了更好地评估，我们选取另外两种方法作为比较对象。第一种方法是基于完整数据(Complete Case, CC)分析的SJMA方法( [9] )，记此方法为CC-SJMA。第二种方法是响应变量缺失下的M-JMA方法( [3] [6] )。

考虑数据生成过程为

$y_i={\mu }_i+e_i={\displaystyle {\sum }_{j=1}^p{\theta }_j{x}_{ij}}+e_i,$

其中 $p=5$ ， $\theta ={\left(0.3,0.6,0.3,0.3,0.3\right)}^{\text{T}}$ ， $x_{i1}=1$ ，其余的 $x_{ij}$ 是独立同分布的标准正态随机变量， $e_i~0.5N\left(0,{\sigma }^2\right)+0.5N\left(0,3{\sigma }^2\right)$ ， $i=1,\cdots ,n$ ， $j=2,3,4,5$ 。通过选取合适的 ${\sigma }^2$ 将 $R^2$ 控制在集合 $\left\{0.1,0.2,\cdots ,0.9\right\}$ 上变动，其中：

$R^2=\frac{n^{-1}{\displaystyle {\sum }_{i=1}^n{\left({\mu }_i-\bar{\mu }\right)}^2}}{{\sigma }^2+n^{-1}{\displaystyle {\sum }_{i=1}^n{\left({\mu }_i-\bar{\mu }\right)}^2}},\bar{\mu }=n^{-1}{\displaystyle {\sum }_{i=1}^n{\mu }_i}$ 。

选择概率函数 $\pi \left(x_i\right)$ 设定为

$\pi \left(x_i\right)=\Phi \left({\alpha }_1+{\alpha }_2{x}_{i2}\right)$ ，

其中 $\Phi (\cdot )$ 为标准正态分布的分布函数。选取两种不同的 $\alpha$ 的值： $\alpha ={\left(0.6,0.5\right)}^{\text{T}}$ 和 $\alpha ={\left(0,0.5\right)}^{\text{T}}$ ，使得相应的平均缺失率分别大约为0.3和0.5。对M-JMA方法，考虑所有包含截距项的候选模型，故共有 $M=2^5$ 个候选模型。

为了评估上述三种方法的表现，采用如下的均方误差(Mean Square Error, MSE)作为度量标准，定义如下：

$\text{MSE}={\displaystyle {\sum }_{d=1}^D{\Vert {\mu }^{\left(d\right)}\left(w\right)-{\mu }^{\left(d\right)}\Vert }^2}/D$ ，

其中D是试验的重复次数， ${\mu }^{\left(d\right)}$ 表示 $\mu$ 在第d次试验中的值， ${\mu }^{\left(d\right)}\left(w\right)$ 是相应于 ${\mu }^{\left(d\right)}$ 的模型平均估计。重复次数D设定为200，考虑的样本量为n = 50,150。为了便于比较，考虑基于完全数据的SJMA方法，并将其视为“黄金标准”用于比较中，记此方法为CD-SJMA。我们分别计算了CC-SJMA，M-JMA和所提M-SJMA方法的MSE除以CD-SJMA方法的MSE的比值，并报告最终的比值，模拟结果如图1所示。

Figure 1. The MSE ratios of different model averaging estimators

图1. 不同模型平均估计所对应的MSE比值

Table 1. The computational times (seconds) of different model averaging estimators with different feature dimension p

表1. 不同特征维度p下模型平均估计的计算时间(秒)

从图1中可以看出，本文所提出的M-SJMA方法的MSE在绝大多数情况下都是最小的。在考虑的所有情况下，M-SJMA的表现明显优于CC-SJMA的表现。并且在大部分情况下，M-SJMA都要比M-JMA表现更出色，特别当样本量较大时，如n = 150。这表明随着样本量的增大，M-JMA会因为计算成本增加产生更多的计算误差，从而导致预测精度降低。与M-JMA相比，CC-SJMA的表现大多数都很差。但是当缺失率较低，样本量较大(如n = 150，缺失率 = 0.3)时，M-JMA可能比CC-SJMA表现更差。这也说明当计算量增大时，可伸缩的模型平均方法在预测精度上是优于传统的模型平均方法的。

另外，针对不同的预测变量个数p和样本量n，我们比较了三种方法的计算时间，结果见表1。从表1中可以明显看出，随着p的增大，M-JMA的计算时间呈指数级增长，特别当p大于13时，由于权重向量维度过大，计算机无法进行优化计算。但是，对可伸缩的模型平均方法，如M-SJMA和CC-SJMA，它们的计算时间都非常短。即使当n = 300，p = 25时，M-SJMA和CC-SJMA的计算时间都小于273.3和318.7秒。

总之，模拟结果表明所提M-SJMA方法的有限样本性质要优于另外两种对比的方法。M-SJMA方法在预测精度和计算成本方面均表现优越。

4. 结论

基于逆概率加权，奇异值分解和Jackknife准则，本文研究了响应变量缺失下线性模型的可伸缩的模型平均估计问题。在理论上证明了所提M-SJMA方法的渐近最优性，并通过数值模拟进一步验证了所提方法的有限样本性质。未来可考虑将该方法推广到半参数模型的情况。

附录

定理1的证明：令 $P=U{U}^{\text{T}}$ ，通过使用与文献 [9] 中定理1类似的证明技巧，我们有

${\inf }_{\tilde{w}\in {\tilde{W}}^{2^p}}{\tilde{L}}_{\hat{\pi }}\left(\tilde{w}\right)\ge {\Vert P\mu -\mu \Vert }^2$ ，

${\inf }_{w\in W}{L}_{\hat{\pi }}\left(w\right)\le {\Vert P{z}_{\hat{\pi }}-\mu \Vert }^2\le {\Vert P\mu -\mu \Vert }^2+e_{\pi }^{\text{T}}P{e}_{\pi }+{\left(z_{\hat{\pi }}-z_{\pi }\right)}^{\text{T}}P\left(z_{\hat{\pi }}-z_{\pi }\right)+2\Vert z_{\hat{\pi }}-z_{\pi }\Vert \Vert P{e}_{\pi }\Vert$ 。

在条件(C1)、(C2)下，存在常数c，使得 ${\max }_{1\le i\le n}E\left(e_{\pi ,i}^4\right)\le c$ ，则

$e_{\pi }^{\text{T}}P{e}_{\pi }=E\left(e_{\pi }^{\text{T}}P{e}_{\pi }\right)+O_P\left\{\sqrt{Var\left(e_{\pi }^{\text{T}}P{e}_{\pi }\right)}\right\}=O_P\left(p\right)$ 。 (A.1)

另外，由文献 [7] 的引理1，有 ${\Vert z_{\hat{\pi }}-z_{\pi }\Vert }^2=O_P\left(1\right)$ ，可得

${\left(z_{\hat{\pi }}-z_{\pi }\right)}^{\text{T}}P\left(z_{\hat{\pi }}-z_{\pi }\right)\le {\lambda }_{\max }\left(P\right){\Vert z_{\hat{\pi }}-z_{\pi }\Vert }^2\le {\Vert z_{\hat{\pi }}-z_{\pi }\Vert }^2=O_p\left(1\right)$ ，

则有 ${\inf }_{w\in W}{L}_{\hat{\pi }}\left(w\right)\le {\Vert P\mu -\mu \Vert }^2+O_P\left(p\right)$ ，且 $$ 。

结合上述结果和条件(C3)，有

$1\le \frac{{\inf }_{w\in W}{L}_{\hat{\pi }}\left(w\right)}{{\inf }_{\tilde{w}\in {\tilde{W}}^{2^p}}{\tilde{L}}_{\hat{\pi }}\left(\tilde{w}\right)}\le \frac{{\Vert P\mu -\mu \Vert }^2+O_P\left(p\right)}{{\Vert P\mu -\mu \Vert }^2}\le 1+\frac{O_P\left(p\right)}{{\Vert P\mu -\mu \Vert }^2}\le 1+\frac{O_P\left(p\right)}{{\xi }_{\pi }-O_P\left(p\right)}=1+o_P\left(1\right)$ ，

故定理1得证。

定理2的证明：记 $A\left(w\right)=I_n-H\left(w\right)$ ， ${\Sigma }_{\pi }=E\left(e_{\pi }{e}_{\pi }^{\text{T}}\right)$ 。经过简单的计算，有

$PH\left(w\right)=H\left(w\right)$ ， $R_{\pi }\left(w\right)={\Vert A\left(w\right)\mu \Vert }^2+tr\left(H\left(w\right){\Sigma }_{\pi }H\left(w\right)\right)\ge {\Vert A\left(w\right)\mu \Vert }^2$ ，

$L_{\pi }\left(w\right)-R_{\pi }\left(w\right)={\Vert H\left(w\right)e_{\pi }\Vert }^2-2e_{\pi }^{\text{T}}H\left(w\right)A\left(w\right)\mu -tr\left(H\left(w\right){\Sigma }_{\pi }H\left(w\right)\right)$ 。 (A.2)

首先， $J_{\hat{\pi }}\left(w\right)={\Vert z_{\hat{\pi }}-\hat{H}\left(w\right)z_{\hat{\pi }}\Vert }^2={\Vert \left(z_{\hat{\pi }}-\mu \right)+\left(\mu -H\left(w\right)z_{\hat{\pi }}\right)-\left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}\Vert }^2$ 可表示为

$J_{\hat{\pi }}\left(w\right)=L_{\hat{\pi }}\left(w\right)+t_{\hat{\pi }}\left(w\right)+{\Vert z_{\hat{\pi }}-\mu \Vert }^2+2e_{\pi }^{\text{T}}\mu -2e_{\pi }^{\text{T}}P\mu$ ，

其中

$\begin{array}{c}{t}_{\hat{\pi }}\left(w\right)={\Vert \left(\hat{H}(w)-H(w)\right)z_{\hat{\pi }}\Vert }^2-2{\left(\mu -H\left(w\right)z_{\hat{\pi }}\right)}^{\text{T}}\left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}\\ -2{\left(z_{\hat{\pi }}-z_{\pi }\right)}^{\text{T}}\left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}-2e_{\pi }^{\text{T}}\left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}\\ +2{\left(z_{\hat{\pi }}-z_{\pi }\right)}^{\text{T}}\left(\mu -H\left(w\right)z_{\hat{\pi }}\right)+2e_{\pi }^{\text{T}}\left\{PA\left(w\right)\mu -H\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)-H\left(w\right)e_{\pi }\right\},\end{array}$

故可知 ${\hat{w}}_{\hat{\pi }}=\arg {\min }_{w\in W}{L}_{\hat{\pi }}\left(w\right)+t_{\hat{\pi }}\left(w\right)$ 。据此，要证定理2成立，只用证：

${\sup }_{w\in W}\left|\frac{L_{\hat{\pi }}\left(w\right)}{R_{\pi }\left(w\right)}-1\right|=o_P\left(1\right)$ ， (A.3)

${\sup }_{w\in W}\frac{\left|t_{\hat{\pi }}\left(w\right)\right|}{R_{\pi }\left(w\right)}=o_P\left(1\right)$ 。 (A.4)

对(A.3)，由于 $L_{\hat{\pi }}\left(w\right)={\Vert H\left(w\right)z_{\hat{\pi }}-\mu \Vert }^2=L_{\pi }\left(w\right)+2{\left(H\left(w\right)z_{\pi }-\mu \right)}^{\text{T}}H\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)+{\Vert H\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)\Vert }^2$ ，则有

$\left|\frac{L_{\hat{\pi }}\left(w\right)}{R_{\pi }\left(w\right)}-1\right|\le \left|\frac{L_{\pi }\left(w\right)}{R_{\pi }\left(w\right)}-1\right|+\frac{2{\left\{L_{\pi }\left(w\right)\right\}}^{1/2}\Vert H\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)\Vert }{R_{\pi }\left(w\right)}+\frac{{\Vert H\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)\Vert }^2}{R_{\pi }\left(w\right)}$ ， (A.5)

由(A.2)，若(i) ${\sup }_{w\in W}{\Vert H\left(w\right)e_{\pi }\Vert }^2/R_{\pi }\left(w\right)=o_P\left(1\right)$ ；(ii) ${\sup }_{w\in W}\left|e_{\pi }^{\text{T}}H\left(w\right)A\left(w\right)\mu \right|/R_{\pi }\left(w\right)=o_P\left(1\right)$ ，(iii) ${\sup }_{w\in W}tr\left(H\left(w\right){\Sigma }_{\pi }H\left(w\right)\right)/R_{\pi }\left(w\right)=o_P(\; 1\; )$

成立，则有

${\sup }_{w\in W}\left|L_{\pi }\left(w\right)/R_{\pi }\left(w\right)-1\right|=o_P\left(1\right)$ 。 (A.6)

事实上，由于 $H\left(w\right)={\displaystyle {\sum }_{j=1}^p{w}_j{u}_{\left(j\right)}{u}_{\left(j\right)}^{\text{T}}}\le P$ ， $H^2\left(w\right)={\left\{{\displaystyle {\sum }_{j=1}^p{w}_j{u}_{\left(j\right)}{u}_{\left(j\right)}^{\text{T}}}\right\}}^2={\displaystyle {\sum }_{j=1}^p{w}_j^2{u}_{\left(j\right)}{u}_{\left(j\right)}^{\text{T}}}\le P$ ，这里的不等号是Loewner序：即矩阵 $A\le B$ 表示当且仅当 $B-A$ 为半正定矩阵。据此，类似于文献 [9] 中引理1中的证明，由条件(C2)，(C3)和结论(A.1)，可知(i)~(iii)成立。

另外，注意到

${\Vert H\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)\Vert }^2={\left(z_{\hat{\pi }}-z_{\pi }\right)}^{\text{T}}{H}^2\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)\le {\left(z_{\hat{\pi }}-z_{\pi }\right)}^{\text{T}}P\left(z_{\hat{\pi }}-z_{\pi }\right)\le {\Vert z_{\hat{\pi }}-z_{\pi }\Vert }^2$ , (A.7)

由条件(C2)，(C3)和 ${\Vert z_{\hat{\pi }}-z_{\pi }\Vert }^2=O_p\left(1\right)$ ，不难证明

${\sup }_{w\in W}{\Vert z_{\hat{\pi }}-z_{\pi }\Vert }^2/R_{\pi }\left(w\right)=o_P\left(1\right)$ 。 (A.8)

从而，结合(A.5)~(A.8)的结果(A.3)得证。

对(A.4)，我们有

$\begin{array}{c}\left|t_{\hat{\pi }}\left(w\right)\right|\le {\Vert \left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}\Vert }^2+2{\left\{L_{\hat{\pi }}\left(w\right)\right\}}^{1/2}\Vert \left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}\Vert \\ +2\Vert z_{\hat{\pi }}-z_{\pi }\Vert \Vert \left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}\Vert +2\left|e_{\pi }^{\text{T}}\left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}\right|\\ +2{\Vert \left(z_{\hat{\pi }}-z_{\pi }\right)\Vert }^{\text{T}}{\left\{L_{\hat{\pi }}\left(w\right)\right\}}^{1/2}+2\left|e_{\pi }^{\text{T}}PA(w)\mu \right|+2\left|e_{\pi }^{\text{T}}H\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)\right|+2e_{\pi }^{\text{T}}H\left(w\right)e_{\pi },\end{array}$

这一结果与(A.3)、(A.8)相结合，不难看出，要证(A.4)成立，我们只用证

${\sup }_{w\in W}{\Vert \left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}\Vert }^2/R_{\pi }\left(w\right)=o_P\left(1\right)$ ， (A.9)

${\sup }_{w\in W}\left|e_{\pi }^{\text{T}}\left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}\right|/R_{\pi }\left(w\right)=o_P\left(1\right)$ ， (A.10)

${\sup }_{w\in W}\left|e_{\pi }^{\text{T}}PA\left(w\right)\mu \right|/R_{\pi }\left(w\right)=o_P\left(1\right)$ ， (A.11)

${\sup }_{w\in W}\left|e_{\pi }^{\text{T}}H\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)\right|/R_{\pi }\left(w\right)=o_P\left(1\right)$ ， (A.12)

${\sup }_{w\in W}{e}_{\pi }^{\text{T}}H\left(w\right)e_{\pi }/R_{\pi }\left(w\right)=o_P\left(1\right)$ 。 (A.13)

记 $Q_j$ 为对角元素为 $h_{j,ii}{\left(1-h_{j,ii}\right)}^{-1}$ 的n阶对角阵， ${{\rm A}}_j={{\rm I}}_n-{{\rm H}}_j$ ，则 ${\hat{H}}_j=H_j-Q_j{A}_j$ , $j=1,\cdots ,p$ 。由条件(C4)可知 $h_{j,ii}=u_{ij}^2\le c_3{n}^{-1}$ ，且在条件(C1)和(C2)下，有 $n^{-1}{\Vert z_{\hat{\pi }}\Vert }^2\le n^{-1}\left({\Vert z_{\hat{\pi }}-z_{\pi }\Vert }^2+{\Vert z_{\pi }\Vert }^2\right)=O_P\left(1\right)$ ， $n^{-1}{\Vert e_{\pi }\Vert }^2=O_P\left(1\right)$ 。因此，我们有

$\begin{array}{l}{\sup }_{w\in W}{\Vert \left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}\Vert }^2/R_{\pi }\left(w\right)\\ ={\sup }_{w\in W}{\displaystyle {\sum }_{j=1}^p{\displaystyle {\sum }_{k=1}^p{w}_j{w}_k}}{z}_{\hat{\pi }}^{\text{T}}{A}_j{Q}_j{Q}_k{A}_k{z}_{\hat{\pi }}/R_{\pi }\left(w\right)\\ ={\sup }_{w\in W}{\displaystyle {\sum }_{j=1}^p{\displaystyle {\sum }_{k=1}^p{w}_j{w}_k}}{z}_{\hat{\pi }}^{\text{T}}\left(\frac{A_j{Q}_j{Q}_k{A}_k+A_k{Q}_k{Q}_j{A}_j}{2}\right)z_{\hat{\pi }}/R_{\pi }\left(w\right)\\ \le {\sup }_{w\in W}{\displaystyle {\sum }_{j=1}^p{\displaystyle {\sum }_{k=1}^p{w}_j{w}_k}}{\lambda }_{\max }\left(A_j{Q}_j{Q}_k{A}_k\right){\Vert z_{\hat{\pi }}\Vert }^2/R_{\pi }(w)\\ \le {\xi }_{\pi }^{-1}{p}^2{\left(\frac{c_3{n}^{-1}}{1-c_3{n}^{-1}}\right)}^2{\Vert z_{\hat{\pi }}\Vert }^2=O_P\left({\xi }_{\pi }^{-1}p\right)=o_P(\; 1\; )\end{array}$

$\begin{array}{l}{\sup }_{w\in W}\left|e_{\pi }^{\text{T}}\left(\hat{H}\left(w\right)-H\left(w\right)\right)z_{\hat{\pi }}\right|/R_{\pi }\left(w\right)\\ \le {\sup }_{w\in W}{\displaystyle {\sum }_{j=1}^p{w}_j\left|e_{\pi }^{\text{T}}{Q}_j{A}_j{z}_{\hat{\pi }}\right|}/R_{\pi }\left(w\right)\\ \le {\sup }_{w\in W}{\displaystyle {\sum }_{j=1}^p{w}_j\Vert Q_j{e}_{\pi }\Vert \Vert A_j{z}_{\hat{\pi }}\Vert }/R_{\pi }\left(w\right)\\ \le {\xi }_{\pi }^{-1}p\frac{c_3{n}^{-1}}{1-c_3{n}^{-1}}\Vert e_{\pi }\Vert \Vert z_{\hat{\pi }}\Vert =O_P\left({\xi }_{\pi }^{-1}p\right)=o_P(\; 1\; )\end{array}$

$\begin{array}{l}{\sup }_{w\in W}\left|e_{\pi }^{\text{T}}PA\left(w\right)\mu \right|/R_{\pi }\left(w\right)\\ \le {\sup }_{w\in W}\Vert P{e}_{\pi }\Vert \Vert A\left(w\right)\mu \Vert /R_{\pi }\left(w\right)\\ \le {\sup }_{w\in W}\Vert P{e}_{\pi }\Vert /R_{\pi }^{1/2}\left(w\right)=O_P\left({\xi }_{\pi }^{-1/2}{p}^{1/2}\right)=o_P(\; 1\; )\end{array}$

$\begin{array}{l}{\sup }_{w\in W}\left|e_{\pi }^{\text{T}}H\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)\right|/R_{\pi }\left(w\right)\\ ={\sup }_{w\in W}\left|e_{\pi }^{\text{T}}PH\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)\right|/R_{\pi }\left(w\right)\\ \le {\sup }_{w\in W}\Vert P{e}_{\pi }\Vert \Vert H\left(w\right)\left(z_{\hat{\pi }}-z_{\pi }\right)\Vert /R_{\pi }\left(w\right)\\ \le {\sup }_{w\in W}{\left(e_{\pi }^{\text{T}}P{e}_{\pi }\right)}^{1/2}{\left\{{\left(z_{\hat{\pi }}-z_{\pi }\right)}^{\text{T}}P\left(z_{\hat{\pi }}-z_{\pi }\right)\right\}}^{1/2}/R_{\pi }\left(w\right)\\ =O_P\left({\xi }_{\pi }^{-1}{p}^{1/2}\right)=o_P(\; 1\; )\end{array}$

${\sup }_{w\in W}{e}_{\pi }^{\text{T}}H\left(w\right)e_{\pi }/R_{\pi }\left(w\right)\le {\sup }_{w\in W}{e}_{\pi }^{\text{T}}P{e}_{\pi }/R_{\pi }\left(w\right)=O_P\left({\xi }_{\pi }^{-1}p\right)=o_P(\; 1\; )$

故(A.9)~(A.13)成立，至此(A.4)证毕。进而定理2得证。