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

A scalable model⁃averaging method with right censored responses

摘要

给出了响应变量随机右删失情形下线性模型的可伸缩模型平均方法。该方法通过逆概率删失加权合成新的响应变量，并利用奇异值分解对原模型进行转换，使得所构建的方法仅需考虑p个候选模型。然后利用Mallows准则和Jackknife准则选择最优权重，且权重之和没有等于1的约束条件，并在理论上证明了对应的模型平均估计量的渐近最优性。数值研究进一步表明所提方法在预测精度和计算时间上的优良性。

Abstract

This paper concerns a scalable model⁃averaging method with right⁃censored responses. By using inverse probability censoring weighting to synthesize new responses and a singular value decomposition to transform the original models， this method enables us to estimate the weights by considering maximum of p covariate models. The Mallows criterion and the Jackknife criterion are applied to the selection of weights without the standard constraint that the weights sum to one. The theoretical results show that the calculated weights exhibit asymptotic optimality in the sense of achieving the lowest possible least squares error. Numerical studies demonstrate the excellent averaging performance of the proposed model in terms of both predictive accuracy and computational time.

Graphical abstract

关键词

可伸缩的 / 模型平均 / 响应变量右删失 / 奇异值分解 / 渐近最优性

Key words

scalable / model averaging / right censored responses / singular value decomposition / asymptotic optimality

引言

模型平均是一种处理模型不确定性和提高预测精度的有效方法。该方法将一系列候选模型用适当的权重组合在一起，将每个候选模型中包含的信息纳入模型预测，从而得到更加稳健和精确的估计。模型平均主要分为贝叶斯模型平均和频率模型平均两个方向，其中基于渐近最优性的各种频率模型平均方法受到越来越多的关注和研究，如基于Mallows准则的Mallows model averaging（MMA）方法^［1-2］、基于Jackknife准则的Jackknife model averaging（JMA）方法^［3-4］、异方差数据下稳健的HR{C}_p方法^［5］、高维数据下模型平均方法^［6］等。这些方法大多考虑完整数据的情形，但是在实际问题中经常会遇到不完整数据。如在生存分析中，生存时间经常会因各种原因出现右删失现象。针对删失数据，在局部误设定的框架下，文献［7-9］分别考虑了线性模型、分位数回归模型和部分线性分位数回归模型的模型平均问题。在去掉局部误设定的约束后，当响应变量右删失时，Yan等^［10］研究了高维数据下的模型平均，Liang等^［11］、Dong等^［12］和Hu等^［13］分别研究了线性模型及部分线性模型的模型平均估计，他们所建立方法都具有渐近最优性，且有很好的预测性能。

但是，当模型包含的协变量个数p较大时，上述方法需要考虑的候选模型可能高达{\mathrm{2}}^p个，这给权重选择和模型估计都带来了极大的挑战。为了解决候选模型数量较大时带来的计算困难问题，Zhu等^［14］提出了可伸缩的模型平均估计，该方法利用奇异值分解（SVD） 将模型进行转换，转换后仅需考虑p个候选模型。受文献［14］的启发，本文考虑了响应变量右删失下线性模型的可伸缩模型平均估计。我们通过逆概率删失加权法合成新的响应变量，利用SVD方法对模型进行转换，使得转换后的模型与原模型有相同的预测能力，并且候选模型的个数降为p个，从而在不降低预测效率的前提下，显著地提高了计算速度。进一步通过Mallows准则^［1-2］和Jackknife准则^［3-4］选择权重，证明了所得方法具有渐近最优性，且在渐近意义下其最小二次损失不低于传统的{\mathrm{2}}^p个权重的模型平均估计。最后通过模拟实验和实例分析进一步验证了所提方法的优良性。

1 方法和理论性质

1.1 模型和记号

令\left\{\left(y_i,{\mathit{x}}_i\right),i=\mathrm{1,2},\cdots ,n\right\}是n个相互独立且同分布的样本，其中y_i是响应变量，可能是我们感兴趣的生存时间T_i的单调变换，如y_i=\mathrm{l}\mathrm{n}{T}_i，{\mathit{x}}_i={\left(x_{i\mathrm{1}},x_{i\mathrm{2}},\cdots ,x_{ip}\right)}^{\mathrm{T}}是p维协变量。本文考虑的情形，且假设样本由线性模型生成，即

\mathit{y}=\mathit{\mu }+\mathit{e}=\mathit{X}\mathit{\theta }+\mathit{e}（1）

式中，\mathit{\theta }为未知的p维回归系数，\mathit{e}={\left(e_{\mathrm{1}},e_{\mathrm{2}},\cdots ,e_n\right)}^{\mathrm{T}}为随机误差，且满足E\left(e_i\left|{\mathit{x}}_i\right.\right)=\mathrm{0}，E\left(e_i^{\mathrm{2}}\left|{\mathit{x}}_i\right.\right)={\sigma }_i^{\mathrm{2}},\mathit{y}={\left(y_{\mathrm{1}},y_{\mathrm{2}},\cdots ,y_n\right)}^{\mathrm{T}},\mathit{\mu }={\left({\mu }_{\mathrm{1}},{\mu }_{\mathrm{2}},\cdots ,{\mu }_n\right)}^{\mathrm{T}},\mathit{X}={\left({\mathit{x}}_{\mathrm{1}},{\mathit{x}}_{\mathrm{2}},\cdots ,{\mathit{x}}_n\right)}^{\mathrm{T}}。

在很多生物医学研究中，y_i经常由于各种原因不能被完全观测，一种常见的情形是y_i被删失时间C_i随机右删失。此时可观测到的数据为\left\{\left(z_i,{\delta }_i,{\mathit{x}}_i\right),i=\mathrm{1,2}，\cdots ,n\right\}，其中z_i=\mathrm{m}\mathrm{i}\mathrm{n}\left(y_i,C_i\right)，{\delta }_i=I\left(y_i\le C_i\right)，I\left(\cdot \right)为示性函数。假设F\left(\cdot \right)和G\left(\cdot \right)分别为y_i和C_i的分布函数。基于观测数据，本文旨在给出\mathit{\mu }的一种模型平均估计。

为了处理删失数据，定义合成的响应变量z_{G,i}={\delta }_i{z}_i/\left\{\mathrm{1}-G(z_i-)\right\}，其中G(z_i-)为函数G\left(\cdot \right)在z_i处的左极限。显然

E\left(z_{G,i}\left|{\mathit{x}}_i\right.\right)=E\left(y_i\left|{\mathit{x}}_i\right.\right)={\mu }_i，i=\mathrm{1,2}，\cdots ,n

因此模型（1）可写为

{\mathit{z}}_G=\mathit{\mu }+{\mathit{e}}_G=\mathit{X}\mathit{\theta }+{\mathit{e}}_G（2）

式中，{\mathit{z}}_G={\left(z_{G,\mathrm{1}},z_{G,\mathrm{2}},\cdots ,z_{G,n}\right)}^{\mathrm{T}}，{\mathit{e}}_G={\left(e_{G,\mathrm{1}},e_{G,\mathrm{2}},\cdots ,e_{G,n}\right)}^{\mathrm{T}},随机误差e_{G,i}=z_{G,i}-{\mu }_i，且满足E\left(e_{G,i}\left|{\mathit{x}}_i\right.\right)=\mathrm{0}，\mathrm{v}\mathrm{a}\mathrm{r}\left(e_{G,i}^{\mathrm{2}}\left|{\mathit{x}}_{\mathit{i}}\right.\right)={\sigma }_{G,i}^{\mathrm{2}}。由于G\left(\cdot \right)通常都是未知的，使用Kaplan-Meier估计

{\widehat{G}}_n\left(z\right)=\mathrm{1}-\prod _{z_i\le z}{\left(\frac{n-i}{n-i+\mathrm{1}}\right)}^{\mathrm{1}-{\delta }_i}

去替换G\left(z\right)，上述记号中的下标G对应地替换成\widehat{G}。

1.2 可伸缩的模型平均估计

考虑M个候选模型，其中第m个候选模型定义为

{\mathit{z}}_{\widehat{G}}={\mathit{X}}_m{\mathit{\theta }}_m+{\mathit{e}}_{\widehat{G}},m=\mathrm{1,2}，\cdots ,M（3）

式中，{\mathit{X}}_m是n\times p_m协变量矩阵，{\mathit{\theta }}_m是p_m维回归系数， 记S_m为该模型包含的协变量的指标集。{\mathit{\theta }}_m的最小二乘估计为{\tilde{\mathit{\theta }}}_{\widehat{G},m}={\left({\mathit{X}}_m^{\mathrm{T}}{\mathit{X}}_m\right)}^{-\mathrm{1}}{\mathit{X}}_m^{\mathrm{T}}{\mathit{z}}_{\widehat{G}}，对应得到\mathit{\mu }的估计{\tilde{\mathit{\mu }}}_{\widehat{G},m}={\mathit{X}}_m{\tilde{\mathit{\theta }}}_{\widehat{G},m}。从而得到\mathit{\mu }的模型平均估计量{\tilde{\mathit{\mu }}}_{\widehat{G}}\left(\tilde{\mathit{w}}\right)={\displaystyle \sum _{m=\mathrm{1}}^M}{\tilde{w}}_m{\tilde{\mathit{\mu }}}_{\widehat{G},m}，其中权重向量\tilde{\mathit{w}}=\left({\tilde{w}}_{\mathrm{1}},{\tilde{w}}_{\mathrm{2}},\cdots ,{\tilde{w}}_M\right)\in {\tilde{W}}^M，且

{\tilde{W}}^M=\left\{\tilde{\mathit{w}}\in {\left[\mathrm{0,1}\right]}^M:{\displaystyle \sum _{m=\mathrm{1}}^M}{\tilde{w}}_m=\mathrm{1}\right\}

文献［11］采用Mallows准则^［5］选择权重\tilde{\mathit{w}}， 由此得到的right censored linear model averaging（RCLMA）估计是渐近最优的，且有很好的预测性能。但是，当协变量的个数p较大时，这种基于原模型的模型平均方法可能需要估计{\mathrm{2}}^p个权重，极大地降低了计算速度和预测精度。为了解决候选模型数量较大带来的计算困难，我们尝试采用新的方法去改进。参考文献［14］，引入奇异值分解将协变量矩阵\mathit{X}转换成列正交的协变量矩阵\mathit{U}，转换后候选模型的数量降为p个，且候选模型均为简单的一元回归模型，从而有效地降低了计算压力和时间成本。

假设\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\mathit{X}\right)=p，将\mathit{X}进行SVD分解：\mathit{X}=\mathit{U}\mathit{D}{\mathit{V}}^{\mathrm{T}}，其中{\mathit{U}}^{\mathrm{T}}\mathit{U}={\mathit{V}}^{\mathrm{T}}\mathit{V}={\mathit{I}}_{\mathit{p}}，\mathit{U}为n\times p阶的列正交矩阵，\mathit{D}为对角元素是非负实数的p\times p阶对角矩阵，\mathit{V}是p\times p阶正交矩阵。对应可得到{\mathit{X}}_m={\mathit{U}}_m{\mathit{D}}_m{\mathit{V}}_m^{\mathrm{T}}，则候选模型（3）可转换为

{\mathit{z}}_{\widehat{G}}={\mathit{U}}_m{\mathit{\beta }}_m+{\mathit{e}}_{\widehat{G}},m=\mathrm{1,2}，\cdots ,M（4）

其中{\mathit{\beta }}_m={\mathit{D}}_m{\mathit{V}}_m^{\mathrm{T}}{\mathit{\theta }}_m。记{\mathit{u}}_{(j)}是\mathit{U}的第j列向量，\mathrm{1}\le j\le p，则{\mathit{U}}_m=({\mathit{u}}_{(j_{\mathrm{1}})},{\mathit{u}}_{(j_{\mathrm{2}})},\cdots ,{\mathit{u}}_{(j_{p_m})})，其中j_{\mathrm{1}},j_{\mathrm{2}},\cdots ,j_{p_m}是集合S_m的元素。可得{\mathit{\beta }}_m的最小二乘估计为

{\widehat{\mathit{\beta }}}_{\widehat{G},m}={\left({\mathit{U}}_m^{\mathrm{T}}{\mathit{U}}_m\right)}^{-\mathrm{1}}{\mathit{U}}_m^{\mathrm{T}}{\mathit{z}}_{\widehat{G}}={\mathit{U}}_m^{\mathrm{T}}{\mathit{z}}_{\widehat{G}}={\left({\widehat{\mathit{\beta }}}_{\widehat{G},j_{\mathrm{1}}},{\widehat{\mathit{\beta }}}_{\widehat{G},j_{\mathrm{2}}},\cdots ,{\widehat{\mathit{\beta }}}_{\widehat{G},j_{p_m}}\right)}^{\mathrm{T}}

其中{\widehat{\mathit{\beta }}}_{\widehat{G},j}={\mathit{u}}_{(j)}^{\mathrm{T}}{\mathit{z}}_{\widehat{G}}可看成一元线性回归模型{\mathit{z}}_{\widehat{G}}={\mathit{u}}_{(j)}{\beta }_j+{\mathit{e}}_{\widehat{G}}中参数{\beta }_j的最小二乘估计，j\in S_m。对应可得到\mathit{\mu }的估计

{\widehat{\mathit{\mu }}}_{\widehat{G},m}={\mathit{U}}_m{\widehat{\mathit{\beta }}}_{\widehat{G},m}=\sum _{j\in S_m}{\mathit{u}}_{\left(j\right)}{\widehat{\mathit{\beta }}}_{\widehat{G},j}

这里需注意，由于影子矩阵之间存在关系

{\mathit{X}}_m{\left({\mathit{X}}_m^{\mathrm{T}}{\mathit{X}}_m\right)}^{-\mathrm{1}}{\mathit{X}}_m^{\mathrm{T}}={\mathit{U}}_m{\mathit{U}}^{\mathrm{T}}={\mathit{U}}_m{\left({\mathit{U}}_m^{\mathrm{T}}{\mathit{U}}_m\right)}^{-\mathrm{1}}{\mathit{U}}_m^{\mathrm{T}}

因此转换后的模型（4）与原模型（3）在估计\mathit{\mu }时有相同的预测能力。定义\mathit{\mu }的模型平均估计量

{\widehat{\mathit{\mu }}}_{\widehat{G}}(\tilde{\mathit{w}})={\displaystyle \sum _{m=\mathrm{1}}^M}{\tilde{w}}_m {\widehat{\mathit{\mu }}}_{\widehat{G},m}={\displaystyle \sum _{m=\mathrm{1}}^M}{\tilde{w}}_m\left(\sum _{j\in S_m}{\mathit{u}}_{(j)}{\widehat{\mathit{\beta }}}_{\widehat{G},j}\right)={\displaystyle \sum _{j=\mathrm{1}}^p}\left({\displaystyle \sum _{m=\mathrm{1}}^M}{\tilde{w}}_{m}I(j\in S_m)\right){\mathit{u}}_{(j)}{\widehat{\mathit{\beta }}}_{\widehat{G},j}（5）

这里将w_j={\displaystyle \sum _{m=\mathrm{1}}^M}{\tilde{w}}_{m}I(j\in S_m)\in \left[\mathrm{0,1}\right]设置为权重，则式（5）可写为

{\widehat{\mathit{\mu }}}_{\widehat{G}}(\mathit{w})={\displaystyle \sum _{j=\mathrm{1}}^p}{w}_j {\widehat{\mathit{\mu }}}_{(j)}=\mathit{H}(\mathit{w}){\mathit{z}}_{\widehat{G}}（6）

其中{\widehat{\mathit{\mu }}}_{(j)}={\mathit{u}}_{(j)}{\widehat{\mathit{\beta }}}_{\widehat{G},j},\mathit{H}(\mathit{w})={\displaystyle \sum _{j=\mathrm{1}}^p}{w}_j{\mathit{H}}_j,{\mathit{H}}_j={\mathit{u}}_{(j)}{\mathit{u}}_{(j)}^{\mathrm{T}},注意{\displaystyle \sum _{j=\mathrm{1}}^p}{w}_j ={\displaystyle \sum _{j=\mathrm{1}}^p}{\displaystyle \sum _{m=\mathrm{1}}^M}{\tilde{w}}_{m}I(j\in S_m)未必等于1，权重向量\mathit{w}=(w_{\mathrm{1}},w_{\mathrm{2}},\cdots ,w_p{)}^{\mathrm{T}}\in W=\left\{\mathit{w}\in {\left[\mathrm{0,1}\right]}^p\right\}。

因此，式（6）可看作\mathit{\mu }的可伸缩的模型平均估计量，它把p个一元回归候选模型通过一定的权重组合起来。该估计量一个直接的好处是，只需要选择p个权重，而非传统模型平均方法的{\mathrm{2}}^p个权重，另外，由于模型（3）和（4）有相同的预测能力，因此该方法在不牺牲预测效率的前提下，显著地提高了计算速度。其次，在实际应用中，减小计算成本可能会减少优化过程中额外产生的数值误差，从而有助于提高预测精度。这在模拟实验和实例分析的结果中也将得以印证。

分别采用Mallows准则^［1-2］和Jackknife准则^［3-4］来选择权重\mathit{w}。由于候选模型都只使用了1个协变量，因此可伸缩的模型平均的Mallows准则为

M_{\widehat{G}}\left(\mathit{w}\right)={‖{\widehat{\mathit{\mu }}}_{\widehat{G}}\left(\mathit{w}\right)-{\mathit{z}}_{\widehat{G}}‖}^{\mathrm{2}}+\mathrm{2}{\widehat{\sigma }}_{\widehat{G}}^{\mathrm{2}}{\mathit{w}}^{\mathrm{T}}\cdot \mathrm{1}

式中，{\widehat{\sigma }}_{\widehat{G}}^{\mathrm{2}}={\left(n-p\right)}^{-\mathrm{1}}{‖{\mathit{z}}_{\widehat{G}}-\mathit{P}{\mathit{z}}_{\widehat{G}}‖}^{\mathrm{2}}， ‖\cdot ‖为欧式范数，\mathit{P}=\mathit{U}{\mathit{U}}^{\mathrm{T}}，\mathrm{1}为元素均为1的n维向量。通过最小化上述准则得

{\widehat{\mathit{w}}}_{\widehat{G}}^{(\mathrm{1})}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}{\mathrm{n}}_{\mathit{w}\in W}{M}_{\widehat{G}}\left(\mathit{w}\right)

且定义\mathit{\mu }的响应变量右删失下可伸缩的MMA （RC⁃SMMA）估计量为{\widehat{\mathit{\mu }}}_{\widehat{G}}\left({\widehat{\mathit{w}}}_{\widehat{G}}^{(\mathrm{1})}\right)。

根据删一交叉验证法^［3-4］，可伸缩的模型平均的Jackknife准则为

{\mathit{J}}_{\widehat{G}}\left(\mathit{w}\right)={‖\widehat{\mathit{H}}\left(\mathit{w}\right){\mathit{z}}_{\widehat{G}}-{\mathit{z}}_{\widehat{G}}‖}^{\mathrm{2}}

式中，\widehat{\mathit{H}}(\mathit{w})={\displaystyle \sum _{j=\mathrm{1}}^p}{w}_j{\widehat{\mathit{H}}}_j，{\widehat{\mathit{H}}}_j={\mathit{\Lambda }}_j({\mathit{H}}_j-{\mathit{I}}_n)+{\mathit{I}}_n，h_{j,ii}为{\mathit{{\rm H}}}_j的第i个对角线元素，{\mathit{\Lambda }}_j为对角元素为(\mathrm{1}-h_{j,ii}{)}^{-\mathrm{1}}的n阶对角阵。通过最小化上述准则得

{\widehat{\mathit{w}}}_{\widehat{G}}^{(\mathrm{2})}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}{\mathrm{n}}_{\mathit{w}\in W}{J}_{\widehat{G}}\left(\mathit{w}\right)

定义\mathit{\mu }的响应变量右删失下可伸缩JMA（RC⁃SJMA）估计量为{\widehat{\mathit{\mu }}}_{\widehat{G}}\left({\widehat{\mathit{w}}}_{\widehat{G}}^{(\mathrm{2})}\right)。

1.3 渐近性质

接下来讨论所得估计量的渐近性质。首先定义平方误差损失函数

{\tilde{L}}_G(\tilde{\mathit{w}})={‖{\tilde{\mathit{\mu }}}_G(\tilde{\mathit{w}})-\mathit{\mu } ‖}^{\mathrm{2}},\tilde{\mathit{w}}\in {\tilde{W}}^M

L_G(\mathit{w})={‖{\widehat{\mathit{\mu }}}_G(\mathit{w})-\mathit{\mu } ‖}^{\mathrm{2}},\mathit{w}\in W

和风险函数

R_G(\mathit{w})=E[L_G(\mathit{w})|\mathit{X}]

其中{\widehat{\mathit{\mu }}}_G(\mathit{w})=\mathit{H}(\mathit{w}){\mathit{z}}_G，记{\xi }_G=\mathrm{i}\mathrm{n}{\mathrm{f}}_{\mathit{w}\in W}{R}_G(\mathit{w})。为了证明渐近最优性，假设下列条件成立。

（C1） 令H\left(t\right)=\mathrm{1}-\left(\mathrm{1}-F\left(t\right)\right)\left(\mathrm{1}-G\left(t\right)\right)以及{\tau }_H=\mathrm{i}\mathrm{n}\mathrm{f}\left\{t:H\left(t\right)=\mathrm{1}\right\}，假设\mathrm{1}-G\left({\tau }_H-\right)>\mathrm{0}。

（C2） 存在常数c_{\mathrm{1}},c_{\mathrm{2}}>\mathrm{0}，使得\mathrm{m}\mathrm{a}{\mathrm{x}}_{\mathrm{1}\le i\le n}E({\mathit{e}}_i^{\mathrm{4}})\le c_{\mathrm{1}}和\mathrm{m}\mathrm{a}{\mathrm{x}}_{\mathrm{1}\le i\le n}\left|{\mu }_i\right|\le c_{\mathrm{2}}几乎处处成立。

（C3） {\xi }_G^{ -\mathrm{1}}p\stackrel{\mathrm{a}.\mathrm{s}.}{\to }\mathrm{0}。

（C4） 存在常数c_{\mathrm{3}}>\mathrm{0}，使得\mathrm{m}\mathrm{a}{\mathrm{x}}_{i,j}{u}_{i j}^{\mathrm{2}}\le c_{\mathrm{3}}{n}^{-\mathrm{1}}，其中u_{i j}为矩阵\mathit{U}中的元素。

为了确保Kaplan⁃Meier估计{\widehat{G}}_n\left(\cdot \right)的一致收敛性，条件C1是删失数据分析中常用的条件，条件C2~C4是模型平均估计中常用的条件^［4，14］。

定理1　在条件C1~C3下，当n\to \mathrm{\infty }时，有

\frac{\mathrm{i}\mathrm{n}{\mathrm{f}}_{\mathit{w}\in W}{L}_{\widehat{G}}(\mathit{w})}{\mathrm{i}\mathrm{n}{\mathrm{f}}_{\tilde{\mathit{w}}\in {\tilde{W}}^{{\mathrm{2}}^p}}{\tilde{L}}_{\widehat{G}}(\tilde{\mathit{w}})}\stackrel{P}{\to }\mathrm{1}

其中{\tilde{W}}^{{\mathrm{2}}^p}是当M={\mathrm{2}}^p时的权重向量集合{\tilde{W}}^M。

证明： 参考文献［14］的证明，有

\mathrm{i}\mathrm{n}{\mathrm{f}}_{\tilde{\mathit{w}}\in {\tilde{W}}^{{\mathrm{2}}^p}}{\tilde{L}}_{\widehat{G}}(\tilde{\mathit{w}})\ge {‖\mathit{P}\mathit{\mu }-\mathit{\mu } ‖}^{\mathrm{2}}\begin{array}{l}\mathrm{i}\mathrm{n}{\mathrm{f}}_{\mathit{w}\in W}{L}_{\widehat{G}}(\mathit{w})\le {‖\mathit{P}{\mathit{z}}_{\widehat{G}}-\mathit{\mu } ‖}^{\mathrm{2}}\le {‖\mathit{P} \mathit{\mu }-\mathit{\mu } ‖}^{\mathrm{2}}+\\ {\mathit{e}}_G^{\mathrm{T}}\mathit{P}{\mathit{e}}_G+{\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)}^{\mathrm{T}}\mathit{P} \left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)+\mathrm{2}‖{\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G‖‖\mathit{P}{\mathit{e}}_G‖\end{array}

在条件C1、C2下，由文献［11］引理6.2，有

{‖{\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G ‖}^{\mathrm{2}}=O_P(\mathrm{1})（7）

且存在常数c，使得

\mathrm{m}\mathrm{a}{\mathrm{x}}_{\mathrm{1}\le i\le n}E({\mathit{e}}_{G,i}^{\mathrm{4}})\le c（8）

故有

{\mathit{e}}_G^{\mathrm{T}}\mathit{P}{\mathit{e}}_G=E\left({\mathit{e}}_G^{\mathrm{T}}\mathit{P}{\mathit{e}}_G\right)+O_P\left\{\sqrt[]{\mathrm{V}\mathrm{a}\mathrm{r}({\mathit{e}}_G^{\mathrm{T}}\mathit{P}{\mathit{e}}_G)}\right\}=O_P(p)（9）

{\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)}^{\mathrm{T}}\mathit{P} \left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)\le {\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}(\mathit{P}){‖{\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G ‖}^{\mathrm{2}}\le

{‖{\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G ‖}^{\mathrm{2}}=O_p(\mathrm{1})（10）

则\mathrm{i}\mathrm{n}{\mathrm{f}}_{\mathit{w}\in W}{L}_{\widehat{G}}(\mathit{w})\le {‖\mathit{P} \mathit{\mu }-\mathit{\mu } ‖}^{\mathrm{2}}+O_P(p)，且{\xi }_{\widehat{G}}=\mathrm{i}\mathrm{n}{\mathrm{f}}_{\mathit{w}\in W}E[L_{\widehat{G}}(\mathit{w})]\le {‖\mathit{P} \mathit{\mu }-\mathit{\mu } ‖}^{\mathrm{2}}+O_P(p)\mathrm{。}

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

\begin{array}{l}\mathrm{1}\le \frac{\mathrm{i}\mathrm{n}{\mathrm{f}}_{\mathit{w}\in W}{L}_{\widehat{G}}(\mathit{w})}{\mathrm{i}\mathrm{n}{\mathrm{f}}_{\tilde{\mathit{w}}\in {\tilde{W}}^{{\mathrm{2}}^p}}{\tilde{L}}_{\widehat{G}}(\tilde{\mathit{w}})}\le \frac{{‖\mathit{P} \mathit{\mu }-\mathit{\mu } ‖}^{\mathrm{2}}+O_P(p)}{{‖\mathit{P} \mathit{\mu }-\mathit{\mu } ‖}^{\mathrm{2}}}\le \\ \mathrm{1}+\frac{O_P(p)}{{‖\mathit{P} \mathit{\mu }-\mathit{\mu } ‖}^{\mathrm{2}}}\le \mathrm{1}+\frac{O_P(p)}{{\xi }_G-O_P(p)}=\mathrm{1}+o_P(\mathrm{1})\end{array}

定理1得证。

定理1表明与传统的{\mathrm{2}}^p个权重的模型平均方法相比，可伸缩的模型平均方法的最小平方误差损失是渐近相等的。

定理2 在条件C1~C4下，当n\to \mathrm{\infty }时，有\frac{L_{\widehat{G}}({\widehat{\mathit{w}}}_{\widehat{G}}^{(\mathrm{1})})}{\mathrm{i}\mathrm{n}{\mathrm{f}}_{\tilde{\mathit{w}}\in {\tilde{W}}^{{\mathrm{2}}^p}}{\tilde{L}}_{\widehat{G}}(\tilde{\mathit{w}})}\stackrel{P}{\to }\mathrm{1}，\frac{L_{\widehat{G}}({\widehat{\mathit{w}}}_{\widehat{G}}^{(\mathrm{2})})}{\mathrm{i}\mathrm{n}{\mathrm{f}}_{\tilde{\mathit{w}}\in {\tilde{W}}^{{\mathrm{2}}^p}}{\tilde{L}}_{\widehat{G}}(\tilde{\mathit{w}})}\stackrel{P}{\to }\mathrm{1}\mathrm{。}

证明： 记\mathit{A}(\mathit{w})={\mathit{I}}_n-\mathit{H}(\mathit{w})。且有

\mathit{P}\mathit{H}(\mathit{w})=\mathit{H}(\mathit{w})，R_G(\mathit{w})\ge {‖\mathit{A}(\mathit{w})\mathit{\mu } ‖}^{\mathrm{2}}

由定理1，要使定理2成立，只需证

\frac{L_{\widehat{G}}({\widehat{\mathit{w}}}_{\widehat{G}}^{(\mathrm{1})})}{\mathrm{i}\mathrm{n}{\mathrm{f}}_{\mathit{w}\in W}{L}_{\widehat{G}}(\mathit{w})}\stackrel{P}{\to }\mathrm{1}（11）

\frac{L_{\widehat{G}}({\widehat{\mathit{w}}}_{\widehat{G}}^{(\mathrm{2})})}{\mathrm{i}\mathrm{n}{\mathrm{f}}_{\mathit{w}\in W}{L}_{\widehat{G}}(\mathit{w})}\stackrel{P}{\to }\mathrm{1}（12）

首先，经过计算可得

\begin{array}{l}{M}_{\widehat{G}}\left(\mathit{w}\right)=L_{\widehat{G}}(\mathit{w})+\mathrm{2}{t}_{\widehat{G}}^{(\mathrm{1})}(\mathit{w})+‖{\mathit{z}}_{\widehat{G}}-\mathit{\mu }‖+\mathrm{2}{\mathit{e}}_G^{\mathrm{T}}\mathit{\mu }-\\ \mathrm{2}{\mathit{e}}_G^{\mathrm{T}}\mathit{P}\mathit{\mu }\end{array}

其中

t_{\widehat{G}}^{(\mathrm{1})}(\mathit{w})={\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)}^{\mathrm{T}}\left(\mathit{\mu }-\mathit{H}(\mathit{w}){\mathit{z}}_{\widehat{G}}\right)+{\mathit{e}}_G^{\mathrm{T}}\left\{\mathit{P}\mathit{A}(\mathit{w})\mathit{\mu }-\mathit{H}(\mathit{w})\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)-\mathit{H}(\mathit{w}){\mathit{e}}_G\right\}+{\widehat{\sigma }}_{\widehat{G}}^{\mathrm{2}}{\mathit{w}}^{\mathrm{T}}\mathrm{1}

故{\widehat{\mathit{w}}}_{\widehat{G}}^{(\mathrm{1})}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}{\mathrm{n}}_{\mathit{w}\in W}{L}_{\widehat{G}}(\mathit{w})+t_{\widehat{G}}(\mathit{w})。要使式（11）成立，只需证

\mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}\left|\frac{L_{\widehat{G}}(\mathit{w})}{R_G(\mathit{w})}-\mathrm{1}\right|=o_P(\mathrm{1})（13）

\mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}\frac{\left|t_{\widehat{G}}^{(\mathrm{1})}(\mathit{w})\right|}{R_G(\mathit{w})}=o_P(\mathrm{1})（14）

对于式（13），由于

L_{\widehat{G}}(\mathit{w})-L_G(\mathit{w})={‖\mathit{H}(\mathit{w})\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right) ‖}^{\mathrm{2}}+\mathrm{2}{\left(\mathit{H}(\mathit{w}){\mathit{z}}_G-\mathit{\mu }\right)}^{\mathrm{T}}\mathit{H}(\mathit{w})\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)

故

\left|\frac{L_{\widehat{G}}(\mathit{w})}{R_G(\mathit{w})}-\mathrm{1}\right|\le \left|\frac{L_G(\mathit{w})}{R_G(\mathit{w})}-\mathrm{1}\right|+\frac{{‖\mathit{H}(\mathit{w})\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right) ‖}^{\mathrm{2}}}{R_G(\mathit{w})}+

\frac{\mathrm{2}\{L_G{(\mathit{w})\}}^{\mathrm{1}/\mathrm{2}}‖\mathit{H}(\mathit{w})\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right) ‖}{R_G(\mathit{w})}（15）

注意到\mathit{H}(\mathit{w})\le \mathit{P}，{\mathit{H}}^{ \mathrm{2}}(\mathit{w})\le \mathit{P}，这里矩阵\mathit{A}\le \mathit{B}表示当且仅当\mathit{B}-\mathit{A}为半正定矩阵。类似于文献［2］中（A.3）的证明，有

\mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}\left|L_G(\mathit{w})/R_G(\mathit{w})-\mathrm{1}\right|=o_P(\mathrm{1})（16）

另外，

\begin{array}{l}{‖\mathit{H}(\mathit{w})\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right) ‖}^{\mathrm{2}}={\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)}^{\mathrm{T}}{\mathit{H}}^{\mathrm{2}}(\mathit{w})\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)\le \\ {\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)}^{\mathrm{T}}\mathit{P}\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)\le {‖{\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G ‖}^{\mathrm{2}}\end{array}

由式（7）和条件C3，可得

\mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{‖{\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G ‖}^{\mathrm{2}}/R_G(\mathit{w})=o_P(\mathrm{1})（17）

从而，结合式（15）~（17））的结果，式（13）得证。

对于式（14），有\left|t_{\widehat{G}}^{(\mathrm{1})}(\mathit{w})\right|\le {\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{5}}}{T}_i,其中

T_{\mathrm{1}}=‖\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)‖{\left\{L_{\widehat{G}}(\mathit{w})\right\}}^{\mathrm{1}/\mathrm{2}},T_{\mathrm{2}}=\left|{\mathit{e}}_G^{\mathrm{T}}\mathit{P}\mathit{A}(\mathit{w})\mathit{\mu }\right|

T_{\mathrm{3}}=\left|{\mathit{e}}_G^{\mathrm{T}}\mathit{H}(\mathit{w})\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)\right|，T_{\mathrm{4}}={\mathit{e}}_G^{\mathrm{T}}\mathit{H}(\mathit{w}){\mathit{e}}_G

T_{\mathrm{5}}={\widehat{\sigma }}_{\widehat{G}}^{\mathrm{2}}{\mathit{w}}^{\mathrm{T}}\mathrm{1}

由式（13）和（17），要证式（14）成立，只需证\mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{T}_i/R_G(\mathit{w})=o_P(\mathrm{1}),i=\mathrm{2,3},\mathrm{4,5}。

在式（8）~（10）和条件C1~C3下，有

\begin{array}{l}\mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{T}_{\mathrm{2}}/R_G(\mathit{w})\le \mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{\xi }_G^{-\mathrm{1}}‖\mathit{P}{\mathit{e}}_G‖‖\mathit{A}(\mathit{w})\mathit{\mu }‖\le \\ \mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{\xi }_G^{-\mathrm{1}/\mathrm{2}}‖\mathit{P}{\mathit{e}}_G‖=O_P({\xi }_G^{-\mathrm{1}/\mathrm{2}}{p}^{\mathrm{1}/\mathrm{2}})=o_P(\mathrm{1})\end{array}\begin{array}{l}\mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{T}_{\mathrm{3}}/R_G(\mathit{w})=\mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{\xi }_G^{-\mathrm{1}}\left|{\mathit{e}}_G^{\mathrm{T}}\mathit{P}\mathit{H}(\mathit{w})\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)\right|\le \mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{\xi }_G^{-\mathrm{1}}‖\mathit{P}{\mathit{e}}_G‖‖\mathit{H}(\mathit{w})\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)‖\le \\ \mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{\xi }_G^{-\mathrm{1}}{\left({\mathit{e}}_G^{\mathrm{T}}\mathit{P}{\mathit{e}}_G\right)}^{\mathrm{1}/\mathrm{2}}{\left\{{\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)}^{\mathrm{T}}\mathit{P}\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)\right\}}^{\mathrm{1}/\mathrm{2}}=\\ O_P({\xi }_G^{-\mathrm{1}}{p}^{\mathrm{1}/\mathrm{2}})=o_P(\mathrm{1})\end{array}

\mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{T}_{\mathrm{4}}/R_G(\mathit{w})\le \mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{\xi }_G^{-\mathrm{1}}{\mathit{e}}_G^{\mathrm{T}}\mathit{P}{\mathit{e}}_G=O_P({\xi }_G^{-\mathrm{1}}p)=o_P(\mathrm{1})

\begin{array}{l}{\widehat{\sigma }}_{\widehat{G}}^{\mathrm{2}}=\frac{{‖{\mathit{z}}_{\widehat{G}}-\mathit{P}{\mathit{z}}_{\widehat{G}}‖}^{\mathrm{2}}}{n-p}\le \frac{{‖{\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G‖}^{\mathrm{2}}}{n-p}+\\ \frac{{‖\mathit{P}\left({\mathit{z}}_{\widehat{G}}-{\mathit{z}}_G\right)‖}^{\mathrm{2}}}{n-p}+\frac{{‖{\mathit{z}}_G-\mathit{P}{\mathit{z}}_G‖}^{\mathrm{2}}}{n-p}=O_P(\mathrm{1})\end{array}

\mathrm{s}\mathrm{u}{\mathrm{p}}_{\mathit{w}\in W}{T}_{\mathrm{5}}/R_G(\mathit{w})\le {\widehat{\sigma }}_{\widehat{G}}^{\mathrm{2}}p{\xi }_G^{-\mathrm{1}}=o_P(\mathrm{1})

故式（14）成立，至此式（11）证毕。类似地，可证式（12）成立。

定理2表明被选出的权重{\widehat{\mathit{w}}}_{\widehat{G}}^{(\mathrm{1})}和{\widehat{\mathit{w}}}_{\widehat{G}}^{(\mathrm{2})}所对应的平方误差损失渐近等价于不可获得的最优权重所对应的平方误差损失，这说明所提的RC⁃SMMA、RC⁃SJMA估计在实现最小平方误差损失意义下是渐近最优的。

2 模拟实验与实例分析

本节将通过模拟实验和实例分析考察所提的RC⁃SMMA、RC⁃SJMA方法在有限样本集上的表现，并与模型平均方法RCLMA^［11］、weighted least squares model averaging（WLSMA）^［12］、smoothed Akaike information criterion（SAIC）、smoothed Bayesian information criterion （SBIC）、等权重（equal weight，EW）方法，以及模型选择方法Akaike information criterion（AIC）、Bayesian information criterion （BIC）进行比较。

AIC、BIC、SAIC和SBIC方法的原理可参考文献［15］。对于AIC、BIC方法，第m个候选模型的AIC和BIC评分分别为

{S_{\mathrm{A}\mathrm{I}\mathrm{C}，}}_m=\mathrm{l}\mathrm{n}\left({\widehat{\sigma }}_{\widehat{G},m}^{\mathrm{2}}\right)+\mathrm{2}{n}^{-\mathrm{1}}tr\left\{{\mathit{P}}_m\right\}

{S_{\mathrm{B}\mathrm{I}\mathrm{C},}}_m=\mathrm{l}\mathrm{n}\left({\widehat{\sigma }}_{\widehat{G},m}^{\mathrm{2}}\right)+n^{-\mathrm{1}}tr\left\{{\mathit{P}}_m\right\}\mathrm{l}\mathrm{n}n

式中，{\widehat{\sigma }}_{\widehat{G},m}^{\mathrm{2}}={\left(n-p_m\right)}^{-\mathrm{1}}{‖{\mathit{z}}_{\widehat{G}}-{\mathit{P}}_m{\mathit{z}}_{\widehat{G}}‖}^{\mathrm{2}}。对于第m个候选模型，基于AIC和BIC方法，SAIC和SBIC权重分别为

w_{\mathrm{S}\mathrm{A}\mathrm{I}\mathrm{C},m}=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{\Delta }{S_{\mathrm{A}\mathrm{I}\mathrm{C}，}}_m/\mathrm{2}\right)/{\displaystyle \sum _{j=\mathrm{1}}^M}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{\Delta }{S_{\mathrm{A}\mathrm{I}\mathrm{C}，}}_j/\mathrm{2}\right)w_{\mathrm{S}\mathrm{B}\mathrm{I}\mathrm{C},m}=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{\Delta }{S_{\mathrm{B}\mathrm{I}\mathrm{C}，}}_m/\mathrm{2}\right)/{\displaystyle \sum _{j=\mathrm{1}}^M}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{\Delta }{S_{\mathrm{B}\mathrm{I}\mathrm{C}，}}_j/\mathrm{2}\right)

式中，\mathrm{\Delta }{S_{\mathrm{A}\mathrm{I}\mathrm{C}，}}_m={S_{\mathrm{A}\mathrm{I}\mathrm{C}，}}_m-\mathrm{m}\mathrm{i}{\mathrm{n}}_{m^{\text{'}}}{S_{\mathrm{A}\mathrm{I}\mathrm{C}，}}_{m^{\text{'}}},\mathrm{\Delta }{S_{\mathrm{B}\mathrm{I}\mathrm{C}，}}_m={S_{\mathrm{B}\mathrm{I}\mathrm{C}，}}_m-\mathrm{m}\mathrm{i}{\mathrm{n}}_{m^{\text{'}}}{S_{\mathrm{B}\mathrm{I}\mathrm{C}，}}_{m^{\text{'}}}。在模拟实验中，利用均方误差（mean square error， MSE）作为评价指标，定义

E_{\mathrm{M}\mathrm{S}\mathrm{E}}=\frac{\mathrm{1}}{n}{\displaystyle \sum _{b=\mathrm{1}}^B}{‖{\widehat{\mathit{\mu }}}_{\widehat{\mathit{G}}}^{(b)}\left(\mathit{w}\right)-\mathit{\mu }‖}^{\mathrm{2}}/B

式中，{\widehat{\mathit{\mu }}}_{\widehat{\mathit{G}}}^{(b)}\left(\mathit{w}\right)为第b次实验的模型平均估计量，B=200为实验次数。

例1~例3的数据均产生于y_i={\displaystyle \sum _{j=\mathrm{1}}^p}{\theta }_j{x}_{ij}+{\epsilon }_i，其中x_{i\mathrm{1}}=\mathrm{1}，其余的x_{ij}是独立同分布的标准正态随机变量，{\epsilon }_i~N\left(\mathrm{0},{\sigma }^{\mathrm{2}}\right)，i=\mathrm{1,2}，\cdots ,n，删失时间C_i服从区间［0，k］上的均匀分布。在例1和例2中，样本量取n=100， 200，且取不同的k值使得删失率（cr）大约为0.35和0.50。

例1　考虑全模型是正确的情形。具体地，令p=5，\mathit{\theta }={\left(\mathrm{0.9,0.5,0.3,0.5,0.3}\right)}^{\mathrm{T}}，{\sigma }^{\mathrm{2}}=\mathrm{0.9}\mathrm{。}RCLMA、WLSMA、EW、SAIC和SBIC方法考虑{\mathrm{2}}^{\mathrm{5}}-\mathrm{1}个候选模型，其中最大的候选模型为使用所有协变量的真实模型。图1（a）~（d）分别呈现了不同删失率、不同样本量下9种方法的MSE箱线图，其中方法1~9分别表示RC⁃SJMA、RC⁃SMMA、RCLMA、WLSMA、SAIC、SBIC、EW、AIC、BIC。另外，为了比较计算速度，针对不同的维度p和样本量n，列出了删失率大约为0.35时9种方法的计算时间，模拟结果见表1。

从图1可以看出，本文所提的RC⁃SMMA和RC⁃SJMA在大多数情形下表现出显著的优势。AIC和BIC的表现较差，RCLMA和WLSMA比EW、SAIC和SBIC有更小的MSE。对于大多数的样本量和删失率，RC⁃SMMA和RC⁃SJMA比RCLMA和WLSMA有更小的MSE，这表明即使用以比较的模型平均方法将真实模型纳入候选模型中，可伸缩的模型平均方法在预测精度上依然优于其他比较方法。所提方法不仅有效降低了计算成本，同时提高了预测精度。此外，RC⁃SJMA略优于RC⁃SMMA。

从表1不难看出，随着协变量个数p的增大，RCLMA和WLSMA的计算时间呈指数级增长，当p大于12时，计算机因为维数过大而无法计算。EW、SAIC和SBIC的计算时间也随着p的增大而快速增长。但是，RC⁃SJMA和RC⁃SMMA的计算时间都不超过1 s。即使当n=300，p=20时，RC⁃SJMA和RC⁃SMMA的计算时间仅需要0.328 s和0.051 s，比AIC和BIC的计算速度都更高。由此可见，可伸缩的模型平均方法可以显著提高计算速度。

例2　考虑所有模型都错误指定的情形，即所有候选模型都是错误的。具体地，取p=\mathrm{6}，{\sigma }^{\mathrm{2}}=\mathrm{1}， {\theta }_j=\sqrt[]{\mathrm{2}}/j,j=\mathrm{1,2},\cdots ,\mathrm{5}， {\theta }_{\mathrm{6}}=\mathrm{0.7}生成数据，但是在候选模型构建中仅使用前5个协变量，即第6个协变量不会出现在任何候选模型中。因此，比较模型平均方法考虑{\mathrm{2}}^{\mathrm{5}}-\mathrm{1}个候选模型，本文方法仅考虑5个一元回归候选模型。模拟结果见图2，其中方法1~9与例1中含义相同。

当模型错误指定时，例2得到了和例1相类似的结果。从图2可以看出，RC⁃SMMA和RC⁃SJMA的表现均优于其他比较方法，这表明所提方法在模型错误指定情形下依然有很强的竞争力。

例3　为了从数值模拟的角度验证定理2的结论，取p=\mathrm{3},\mathit{\theta }={\left(\mathrm{0.9,0.6,0.3}\right)}^{\mathrm{T}}，{\sigma }^{\mathrm{2}}=\mathrm{0.4}生成数据， 且取k值使得删失率约为0.50。计算

LR=\frac{L_{\widehat{G}}(\widehat{\mathit{w}})}{\mathrm{i}\mathrm{n}{\mathrm{f}}_{\mathit{w}\in W}{L}_{\widehat{G}}(\mathit{w})}

在样本量n=100，200，400，800，1 600下，重复实验100次。

图3呈现了RC⁃SMMA和RC⁃SJMA方法下LR的均值曲线。从图3中可以明显看出，本文所提的两种方法的LR均值随着样本量的增加逐渐减小且越来越接近于1。这一结果从数值上直观地验证了RC⁃SMMA和RC⁃SJMA方法的渐近最优性。

例4（实例分析）　将RC⁃SMMA、RC⁃SJMA和其他比较方法应用到实际数据分析中。考虑Primary Biliary Cirrhosis （PBC）数据集，该数据集是梅奥（Mayo）诊所收治的424例原发性胆汁性肝硬化患者的生存数据。该PBC数据集可从R语言中的“survival”包中获得，且已被很多文献使用^{［11-12，16］}。本文仅考虑没有缺失值的n=276例患者的数据，每例患者对应的有17个协变量。感兴趣的响应变量是生存时间的对数，其中生存时间是患者从登记入诊到死亡之间的天数。对于接受肝移植或在研究结束时仍然活着的患者，他们的响应变量被认为是右删失的。在这276例患者中，有111例患者死亡，因此删失率约为0.60。

在17个协变量中，Huang等^［16］指出有6个显著的变量：edema（0，无水肿；0.5，未治疗或成功治疗；1，未成功治疗的水肿），bili（血清胆红素，mg/dL），stage（疾病的组织学阶段，分级为1，2，3，4），albumin（白蛋白，g/dL），age（年龄，岁）和alk.phos（碱性磷酸酶，U/L）。利用这p=6个变量构建候选模型。本文方法考虑6个候选模型，其他比较模型平均方法考虑{\mathrm{2}}^{\mathrm{6}}-\mathrm{1}个候选模型。另外，随机选择n_{\mathrm{0}}个患者的数据作为训练集，剩余的n_{\mathrm{1}}=n-n_{\mathrm{0}}个患者的数据作为测试集。使用MSE作为评价指标，即

E_{\mathrm{M}\mathrm{S}\mathrm{E}}=\frac{\mathrm{1}}{n_{\mathrm{1}}}{\displaystyle \sum _{i=n_{\mathrm{0}}+\mathrm{1}}^n}{\left(z_{\widehat{G},i}-{\widehat{\mu }}_i\right)}^{\mathrm{2}}

重复实验100次，这里取n_{\mathrm{0}}=\mathrm{150,200}。表2列出了不同方法在100次重复实验中MSE的均值、中位数和标准差。

从表2可以看出，在n_{\mathrm{0}}和n_{\mathrm{1}}不同的组合中，均值、中位数和标准差的最小值均在所提的RC⁃SMMA或RC⁃SJMA方法中出现，且这两种方法的表现相当。所以RC⁃SMMA和RC⁃SJMA在PBC数据集上的表现总是最优的且结果也是相对较稳定的，这与模拟实验中的结论一致。

3 结束语

本文提出了响应变量随机右删失情形下可伸缩的模型平均方法。该方法仅需考虑p个候选模型，且在不牺牲预测效率的前提下，显著地提高了计算速度。在理论上证明了所提方法的渐近最优性，并通过模拟实验、实例分析验证了该方法在预测精度和计算时间上的优良性，以及从数值角度验证了所提方法的渐近最优性。

在本文中删失指标完全已知，然而由于某些原因，可能会遗漏一些患者的删失指标。因此，未来可考虑将该方法拓展到删失指标缺失的情形。

参考文献

[1] HANSEN B E. Least squares model averaging［J］. Econometrica， 2007， 75（4）： 1175-1189.

[2] WAN A T K， ZHANG X Y， ZOU G H. Least squares model averaging by Mallows criterion［J］. Journal of Econometrics， 2010， 15： 277-283.

[3] HANSEN B E， RACINE J S. Jackknife model averaging［J］. Journal of Econometrics， 2012， 167： 38-46.

[4] ZHANG X Y， WAN A T K， ZOU G H. Model averaging by Jackknife criterion in models with dependent data［J］. Journal of Econometrics， 2013， 174： 82-94.

[5] LIU Q F， OKUI R. Heteroscedasticity⁃robust C_p model averaging［J］. The Econometrics Journal， 2013， 16： 463-472.

[6] ANDO T， LI K C. A model⁃averaging approach for high⁃dimensional regression［J］. Journal of the American Statistical Association， 2014， 109： 254-265.

[7] 孙志猛， 马景义， 苏治.响应变量删失情况下线性模型的FIC模型选择和模型平均［J］.中国科学：数学， 2013，43（7）： 647-661.

[8] SUN Z M， MA J Y， SU Z. FIC model selection and model averaging for linear model with censored response ［J］. Scientia Sinica： Mathematica， 2013，43（7）： 647-661. （in Chinese）

[9] DU J， ZHANG Z Z， XIE T F. Focused information criterion and model averaging in censored quantile regression［J］. Metrika， 2017， 80： 547-570.

[10] SUN Z M， SUN L Q， LU X L， et al. Frequentist model averaging estimation for the censored partial linear quantile regression model［J］. Journal of Statistical Planning and Inference， 2017， 189： 1-15.

[11] YAN X D， WANG H N， WANG W， et al. Optimal model averaging forecasting in high⁃dimensional survival analysis［J］. International Journal of Forecasting， 2021， 37（3）： 1147-1155.

[12] LIANG Z Q， CHEN X L， ZHOU Y Q. Mallows model averaging estimation for linear regression model with right censored data［J］. Acta Mathematicae Applicatae Sinica， English Series， 2022， 38（1）： 5-23.

[13] DONG Q K， LIU B X， ZHAO H. Weighted least squares model averaging for accelerated failure time models［J］. Computational Statistics and Data Analysis， 2023， 184： 107743.

[14] HU G Z， CHENG W H， ZENG J. Optimal model averaging for semiparametric partially linear models with censored data［J］. Mathematics， 2023， 11（3）： 734.

[15] ZHU R， WANG H Y， ZHANG X Y， et al. A scalable frequentist model averaging method［J］. Journal of Business & Economic Statistics， 2023， 41（4）： 1228-1237.

[16] CLAESKENS G， HJORT N L. Model selection and model averaging［M］. Cambridge： Cambridge University Press， 2008.

[17] HUANG J， MA S G， XIE H L. Regularized estimation in the accelerated failure time model with high⁃dimensional covariates［J］. Biometrics， 2006， 62（3）： 813-820.

基金资助

国家自然科学基金(12371219)