加速失效模型下大规模生存数据的分布式模型平均

Distributed Model Averaging for Large-Scale Survival Data under Accelerated Failure Time Models

摘要: 大规模分布式数据的出现，给传统统计分析带来了巨大挑战，也使得分布式统计分析近年来成为统计研究中的热点问题。分布式模型平均方法虽能有效处理大规模数据中的模型不确定性问题，但据作者了解，现有文献仅研究了完全数据的场景，尚未涉及生存数据的分布式模型平均。针对大规模分布式右删失生存数据，利用逆概率删失加权技术，本文提出了加速失效模型下的两种分布式Mallows类型的模型平均方法。两种方法分别采用单轮型估计算法与替代最小二乘的迭代型估计算法进行备选模型的回归系数估计，但在权重选择上均统一采用局部最优权重的简单平均策略。本文通过数值模拟实验验证了所提方法的有效性，并将所提方法用到了实际数据分析。结果表明，所提方法在保持与基于全数据方法相近预测精度的同时，显著提升了计算效率。

Abstract: The emergence of massive distributed data imposes significant challenges to traditional statistical approaches, and also makes distributed statistical learning a hot topic in recent years. Although distributed model averaging methods can effectively handle model uncertainty in these large-scale datas, to the best of our knowledge, existing relevant literature has only investigated the scenarios with complete data and has not yet covered survival data. For the massive distributed survival data with rightly censoring responses, we propose two distributed Mallows-type model averaging methods using the technique of inverse probability of censoring weighting under the accelerated failure time model. These two methods differ in the estimation of the coefficient vector of candidate models—one employs a one-round algorithm while the other uses an iterative algorithm based on surrogate least squares—but both adopt a simple averaging strategy based on locally optimal weights for weight aggregation. The effectiveness of the proposed methods is validated through numerical simulation studies, and the methods are also applied to a real data analysis. The results show that the proposed approaches maintain the prediction accuracy similar to that of full-data methods and yield a marked improvement in computational time.

1. 引言

随着科学技术的发展，大规模数据越来越多地出现在电子商务、社交网络、电子健康系统、公共政策、生物医学等诸多领域。与传统数据相比，大规模数据具有样本量大、数据物理分布于不同节点、数据结构复杂、变量维数高、更新速度快等特点。由于计算和存储的限制，一方面，过大的样本量使得人们无法使用传统的统计分析方法和软件分析大规模数据；另一方面，大规模数据经常物理分散在不同节点、数据中心等。因此，大规模数据为统计学的发展带来了新的挑战和机遇。

大规模数据统计建模和推断研究的一个重要内容是如何用传统的、标准的统计模型和方法来拟合大规模数据。这个问题的困难之处有以下三点：第一，由于样本量过大，大规模数据已不能完全存储于单个计算机内存中；第二，分析大规模数据时，传统统计方法需要的计算时间远远超出了可以接受的范围；第三，大规模数据在物理上经常分布式存储于不同节点、数据中心等，且由于隐私保护数据无法共享。分治策略是处理大规模分布式数据的一种非常自然的框架。近年来，基于分治策略的分布式统计学习方法成为统计学研究的热点，且取得了显著进展。相关的研究进展见综述文章Gao等(2022) [1]、Li等(2024) [2]及其中的参考文献。

在实际数据分析中，数据的真实生成模型一般是未知的，实际工作者因此会面临模型不确定性的问题。作为一种集成学习方法，模型平均能有效处理模型的不确定性问题。模型平均方法主要分为贝叶斯模型平均与频率模型平均两大类，其中频率模型平均依据其对模型参数的假定，可进一步划分为固定参数框架、局部渐近框架两种类型。在线性模型下，Hansen (2007) [3]提出了基于Mallows准则的最小二乘模型平均方法，并证明了权重选择的最优性。Hansen (2007) [3]关于频率模型平均方法中最优权重选择的研究是最优模型平均领域的奠基性工作。该工作不仅激发了后续大量相关研究，也系统性地推动了整个模型平均领域的发展。过去的二十年里，最优模型平均方法的适用性已显著拓展，覆盖了广泛的统计学习模型和多种数据类型，包括广义线性模型(Zhang, 2016) [4]、非线性回归模型(Feng等，2022) [5]、分位数回归模型(Lu和Su，2015) [6]、高维数据(Ando和Li，2014) [7]、缺失数据(Wei等，2021) [8]等。生存数据是一类与生存时间密切相关的复杂数据。该类数据的复杂之处在于由于某些原因，例如病人退出试验、试验的终止等，生存时间往往被右删失。生存数据是一类与事件发生时间密切相关的复杂数据类型，其复杂性主要源于观测过程中的右删失现象，例如因受试者退出试验或研究终止等原因，导致确切的生存时间无法被完整观测。近年来，模型平均的研究也被扩展到生存分析领域，例如He等(2020) [9]、Liang等(2022) [10]、He等(2023) [11]等。

近年来，大规模数据的快速增长推动了对大规模数据模型平均方法的研究。Fang等(2018) [12]针对线性模型与广义线性模型开展了分布式模型平均的初步研究，但该研究并没有给出严格的理论性质的证明。在线性回归模型下，Zhang等(2023) [13]引入了两种基于最小二乘 Mallows 准则的模型平均方法，并严格证明了权重选择的渐近最优性等理论性质。Zhang等(2023) [13]工作的局限性在于需要候选模型为嵌套结构。Xia等(2025) [14]进一步拓展了Zhang等(2023) [13]的工作，允许候选模型为非嵌套结构。Song等(2026) [15]研究了非参数多元可加模型下的分布式Mallows模型平均方法。与现有其他文献不同，Zhang等(2025) [16]首次在局部渐近框架下对分布式模型平均进行了研究。尽管已有部分文献对分布式模型平均进行了研究，但据作者所知，目前尚未有文献探讨面向大规模生存数据的分布式模型平均问题。

针对大规模分布式右删失生存数据，本文融合逆概率删失加权技术、分布式估计算法与Mallows模型平均，提出了两种加速失效模型下的分布式模型平均方法。基于逆概率删失加权的最小二乘损失函数，两种方法首先分别采用单轮型算法与基于替代最小二乘的迭代型算法，对备选模型的回归系数进行估计。在此基础上，进一步在局部节点构造逆概率删失加权的Mallows类型权重准则函数，并通过简单平均的方式，对各局部节点得到的最优模型平均权重进行聚合，最终用于估计或预测。数值模拟实验结果表明，本文所提出的分布式模型平均方法在风险水平上与基于全数据的方法近似相等，但计算所需时间显著降低。本文还将新的方法应用到了一个实际数据分析中。

2. 数据与模型

设有$N$ 个独立个体，第$i$ 个个体的生存时间和(可数无穷维)预测变量分别为$T_i$ 和$x_i={\left(x_{i1},x_{i2},\cdots \right)}^{\text{T}}$ ，其服从加速失效模型

$\log \left(T_i\right)={\mu }_i+{\epsilon }_i$ ，$i=1,2,\cdots ,N$ ， (1)

其中${\mu }_i=\mu \left(x_i\right)={\displaystyle {\sum }_{j=1}^{\infty }{\beta }_j{x}_{ij}}$ ，${\beta }_j$ 为回归系数，${\epsilon }_i$ 为随机误差且满足$E\left({\epsilon }_i|x_i\right)=0$ ，$D\left({\epsilon }_i|x_i\right)={\sigma }^2$ 。考虑模型(1)的$M$ 个线性近似模型，第$m$ 个备选近似模型为

$\log \left(T_i\right)={\displaystyle \sum _{j=1}^{p_m}{x}_{ij}}{\beta }_{j,m}+e_{i,m}=x_{i,m}^{\text{T}}{\beta }_m+e_{i,m}$ ，

其中$x_{i,m}={\left(x_{i1,m},x_{i2,m},\cdots ,x_{i{p}_m,m}\right)}^{\text{T}}$ 为第$m$ 个候选模型的协变量，$p_m$ 为协变量的维数，$e_{i,m}$ 为备选模型的随机误差。

在生物医学等领域的具体研究中，由于研究时间有限、受访者失访等原因，个体的生存时间往往被右删失了。记第$i$ 个个体的右删失时间为$V_i$ ，则实际观测时间为生存时间和删失时间中的较小者。为后续叙述的方便，记$Y_i=\log \left(T_i\right)$ ，$C_i=\log \left(V_i\right)$ ，$Z_i=\min \left\{Y_i,C_i\right\}$ 。令${\delta }_i=I\left(Y_i\le C_i\right)$ 为删失指示变量，其中$I(\cdot )$ 为示性函数。记全部观测个体的数据集为$\mathcal{D}=\left\{\left(Y_i,x_i,{\delta }_i\right),i=1,\cdots ,N\right\}$ 。假设数据集$\mathcal{D}$ 中的数据被存储于$K$ 个局部计算机。不失一般性，假设每个局部计算机存储的数据量相同，记为$n$ ，则$N=Kn$ 。记第$k$ 个局部计算机中的数据集为${\mathcal{D}}_k=\left\{\left(Y_{ki},x_{ki},{\delta }_{ki}\right)|i=1,\cdots ,n\right\}$ 。

根据上述的记号，在k个局部计算机上，数据在正确模型和近似模型下的回归关系分别为

$\log \left(T_{ki}\right)={\displaystyle {\sum }_{j=1}^{\infty }{\beta }_j{x}_{kij}}+{\epsilon }_{ki}$

和

$\log \left(T_{ki}\right)=x_{ki,m}^{\text{T}}{\beta }_m+e_{ki,m}$ 。

3. 分布式模型平均

3.1. 基于单轮型估计的分布式模型平均

单轮型分布式估计是处理大规模数据的一种高效策略，各节点仅需与中心机器进行一次通信即可完成参数聚合。本节基于逆概率删失加权最小二乘损失函数，提出基于单轮型(one-shot, OS)估计的分布式模型平均。

对于第$m$ 个候选模型，在第k个局部计算机上，基于逆概率删失加权的最小二乘估计方法，回归系数向量${\beta }_m$ 估计的损失函数为

$Q_k\left({\beta }_m\right)={\displaystyle \sum _{i=1}^n\frac{{\delta }_i}{1-G\left(Y_{ki}\right)}}{\left(Y_{ki}-x_{ki,m}^{\text{T}}{\beta }_m\right)}^2$ ， (2)

其中$G$ 是$C_i$ 的累积分布函数。但是，由于$G$ 是未知的，上述的损失函数并不能直接使用。一个直观的做法是使用第$k$ 个局部计算机上的数据估计$G$ 。注意到，使用全局数据来估计$G$ 只需要一轮向量的传输，通信效率的损失非常少。因此，首先把各局部计算机上的观测时间和删失指示变量传输到中心计算机，其次在中心计算机上构建$G$ 的全局估计${\widehat{G}}_N$ ，最后把${\widehat{G}}_N$ 再传输给各局部计算机。具体地，令${\widehat{G}}_N$ 是$G$ 的Kaplan-Meier估计，其中

${\widehat{G}}_N\left(y\right)=1-{{\displaystyle \prod _{i=1}^N\left(1-\frac{1}{{\displaystyle \sum _{j=1}^{N}I\left\{Y_j\ge Y_i\right\}}}\right)}}^{\left(1-{\delta }_i\right)I\left(Y_i\le y\right)}$ 。 (3)

将(3)式代入(2)即得可行的局部损失函数

$Q_k\left({\beta }_m\right)={\displaystyle \sum _{i=1}^n\frac{{\delta }_i}{1-{\widehat{G}}_N\left(Y_{ki}\right)}}{\left(Y_{ki}-x_{ki,m}^{\text{T}}{\beta }_m\right)}^2$。 (4)

记$n_k$ 为第$k$ 个局部计算机上未删失观测值的数量，不失一般性，假设前$n_k$ 个观测值未删失，则(4)式可改写为

$Q_k\left({\beta }_m\right)={\displaystyle \sum _{i=1}^{n_k}\frac{1}{1-{\widehat{G}}_N\left(Y_{ki}\right)}}{\left(Y_{ki}-x_{ki,m}^{\text{T}}{\beta }_m\right)}^2$。 (5)

令$Y_{k,un}={\left(Y_{k1},Y_{k2},\cdots ,Y_{k{n}_k}\right)}^{\text{T}}$ ，$X_{k,un}={\left(x_{k1},x_{k2},\cdots ,x_{k{n}_k}\right)}^{\text{T}}$ ，${\widehat{D}}_{k,un}=\text{diag}\left({\widehat{d}}_1,{\widehat{d}}_2,\cdots ,{\widehat{d}}_{n_k}\right)$，其中 ${\widehat{d}}_i=1/\left[1-{\widehat{G}}_N\left(Y_{ki}\right)\right]$ ，则${\beta }_m$ 的逆概率删失加权估计为

${\widehat{\beta }}_{k,m}={\left(X_{k,un,m}^{\text{T}}{\widehat{D}}_{k,un}{X}_{k,un,m}\right)}^{-1}{X}_{k,un,m}^{\text{T}}{\widehat{D}}_{k,un}{Y}_{k,un}$，

其中$X_{k,un,m}$ 是由$X_{k,un}$ 的前$m$ 个列构成的子矩阵。对第$m$ 个候选模型，每个局部计算机将回归系数估计传输到中心计算机，通过简单聚合得到回归系数的单轮型估计${\widehat{\beta }}_{os,m}=K^{-1}{\displaystyle {\sum }_{k=1}^K{\widehat{\beta }}_{k,m}}$ 。

给定各候选模型的回归系数估计和新的预测向量$x={\left(x_1,x_2,\cdots \right)}^{\text{T}}$ ，其回归均值的加权平均预测值为

${\widehat{\mu }}_{os}\left(w\right)={\displaystyle \sum _{m=1}^M{w}_m}{\widehat{\mu }}_{os,m}={\displaystyle \sum _{m=1}^M{w}_m}{x}_m^{\text{T}}{\widehat{\beta }}_{os,m}$，

其中${\widehat{\mu }}_{os,m}=x_m^{\text{T}}{\widehat{\beta }}_{os,m}$，$x_m={\left(x_{1,m},x_{2,m},\cdots ,x_{p_m,m}\right)}^{\text{T}}$ ，$w={\left(w_1,w_2,\cdots ,w_M\right)}^{\text{T}}$ 为权重向量，且$w\in \mathcal{X}=\left\{w\in {\left[0,1\right]}^M:{\displaystyle \sum _{m=1}^M{w}_m}=1\right\}$ 。对于${\widehat{\mu }}_{os}\left(w\right)$ ，一个关键的问题是如何获得数据驱动的权重向量，本文借鉴Zhang等(2020) [17]的方法选择权重。

在第$k$ 个局部计算机上，定义如下的Mallows类型的权重选择函数

$L_k\left(w\right)={\Vert Y_{k,un}-{\widehat{\mu }}_k\left(w\right)\Vert }^2+{\varphi }_N{\widehat{\sigma }}_k^2{\displaystyle \sum _{m=1}^M{w}_m}{p}_m$ ，

其中${\widehat{\mu }}_k\left(w\right)={\displaystyle {\sum }_{m=1}^M{w}_m{\widehat{\mu }}_{k,m}}$ ，${\widehat{\mu }}_{k,m}=X_{k,un,m}^{\text{T}}{\widehat{\beta }}_{os,m}$，${\widehat{\sigma }}_k^2={\Vert Y_k-{\widehat{\mu }}_{k,m^{\ast }}\Vert }^2/\left(n-p_{m^{\ast }}\right)$ ，$m^{\ast }=\arg {\max }_{1\le m\le M}{p}_m$ 。基于损失函数$L_k\left(w\right)$ ，计算第$k$ 个局部计算机上的最优权重${\widehat{w}}_k=\arg {\min }_{w\in \mathcal{X}}{L}_k\left(w\right)$ ，并将$K$ 个最优权重传输给中心计算机，进而通过简单加权平均的方式得到权重的全局估计量$\widehat{w}=K^{-1}{\displaystyle {\sum }_{k=1}^K{\widehat{w}}_k}$ 。新观测的预测向量$x$ 的最终预测值为

${\widehat{\mu }}_{os}\left(\widehat{w}\right)={\displaystyle \sum _{m=1}^M{\widehat{w}}_m}{\widehat{\mu }}_{os,m}={\displaystyle \sum _{m=1}^M{\widehat{w}}_m}{x}_m^T{\widehat{\beta }}_{\text{os},m}$。

3.2. 基于梯度增强迭代估计的分布式模型平均

单轮型分布式估计方法虽然具有较高的通讯效率，但是当局部计算机中数据的样本量较小时，局部估计偏差可能影响最终的估计。梯度增强的分布式迭代估计算法是单轮型分布式估计的重要改进，它在局部计算和全局聚合之间迭代。特别地，在该类算法的每一轮迭代中，每个节点计算机都会基于自身样本，以及通过通信获取的来自其它所有局部计算机的梯度信息，最小化一个修正后的损失函数，从而得到该轮迭代的局部估计(Jordan等(2019) [18]，Fan等(2023) [19])。受Jordan等(2019) [18]和Shamir等(2014) [20]的启发，本小节考虑基于梯度增强迭代(gradient-enhanced iterative, GEI)估计的分布式模型平均。

对于第$m$ 个候选模型，在第$k$ 个局部计算机上，从(5)式可知局部损失函数$Q_k\left({\beta }_m\right)$ 的梯度函数为

$\nabla Q_k\left({\beta }_m\right)=-X_{k,un,m}^{\text{T}}{\widehat{D}}_{k,un}\left(Y_{k,un}-X_{k,un,m}{\beta }_m\right)$.

第$m$ 个候选模型回归系数估计的迭代过程如下：给定第$l$ 次迭代的全局估计值${\widehat{\beta }}_m^{\left(l\right)}$ ，中心计算机首先将其传输给各局部计算机，并局部计算梯度函数$\nabla Q_k\left({\widehat{\beta }}_m^{\left(l\right)}\right)$ ；其次，中心计算机收集各局部计算机的梯度函数值，计算$\nabla Q\left({\beta }_m\right)=K^{-1}{\displaystyle {\sum }_{k=1}^K\nabla Q_k\left({\beta }_m\right)}$ ，再将$\nabla Q\left({\beta }_m\right)$ 分别传播至各局部计算机；再次，在第$k$ 个局部计算机上定义梯度增强的替代损失函数

${\tilde{Q}}_k\left({\beta }_m\right)=Q_k\left({\beta }_m\right)-{\beta }_m^{\text{T}}\left[\nabla Q_k\left({\widehat{\beta }}_m^{\left(l\right)}\right)-\nabla Q\left({\widehat{\beta }}_m^{\left(l\right)}\right)\right]$，

接下来求解第$k$ 个局部计算机的下一次迭代估计${\widehat{\beta }}_{k,m}^{\left(l+1\right)}=\arg {\min }_{{\beta }_m}{\tilde{Q}}_k\left({\beta }_m\right)$ ；最后，各局部计算机将局部估计传输给中心计算机，中心计算机计算全局回归系数估计的第$l+1$ 次迭代值，即

$\begin{array}{c}{\widehat{\beta }}_m^{\left(l+1\right)}=\frac{1}{K}{\displaystyle \sum _{k=1}^K{\widehat{\beta }}_{k,m}^{\left(l+1\right)}}\\ ={\widehat{\beta }}_m^{\left(l\right)}-\frac{1}{K}{\displaystyle \sum _{k=1}^K\left[{\left(X_{k,un,m}^{\text{T}}{\widehat{D}}_{k,un}{X}_{k,un,m}\right)}^{-1}\left(\frac{1}{K}{\displaystyle \sum _{k=1}^K{X}_{k,un,m}^{\text{T}}}{\widehat{D}}_{k,un}{\widehat{\epsilon }}_{km}^{\left(l\right)}\right)\right]}\end{array}$，

其中${\widehat{\epsilon }}_{km}^{\left(l\right)}=Y_{k,un}-X_{k,un,m}{\widehat{\beta }}_m^{\left(l\right)}$ 。假设经过$L$ 次迭代后，迭代过程收敛，则第$m$ 个候选模型回归系数的估计为${\widehat{\beta }}_{iter,m}={\widehat{\beta }}_m^{\left(L\right)}$ 。

类似于第2.1节，给定各候选模型的回归系数估计和新的预测向量$x={\left(x_1,x_2,\cdots \right)}^{\text{T}}$ ，其回归均值的加权平均预测值为

${\widehat{\mu }}_{iter}\left(w\right)={\displaystyle \sum _{m=1}^M{w}_m}{\widehat{\mu }}_{iter,m}={\displaystyle \sum _{m=1}^M{w}_m}{x}_m^{\text{T}}{\widehat{\beta }}_{iter,m}$，

其中${\widehat{\mu }}_{iter,m}=x_m^{\text{T}}{\widehat{\beta }}_{iter,m}$，$x_m={\left(x_{1,m},x_{2,m},\cdots ,x_{p_m,m}\right)}^{\text{T}}$ ，$w={\left(w_1,w_2,\cdots ,w_M\right)}^{\text{T}}$ 为权重向量。在第$k$ 个局部计算机上，定义权重选择函数

$L_k^{\ast }\left(w\right)={\Vert Y_{k,un}-{\widehat{\mu }}_k^{\ast }\left(w\right)\Vert }^2+{\varphi }_N{\widehat{\sigma }}_k^2{\displaystyle \sum _{m=1}^M{w}_m}{p}_m$ ，

其中${\widehat{\mu }}_k^{\ast }\left(w\right)={\displaystyle {\sum }_{m=1}^M{w}_m{\widehat{\mu }}_{k,m}^{\ast }}$ ，${\widehat{\mu }}_{k,m}^{\ast }=X_{k,un,m}^{\text{T}}{\widehat{\beta }}_{iter,m}$，${\widehat{\sigma }}_k^2={\Vert Y_k-{\widehat{\mu }}_{k,m^{\ast }}^{\ast }\Vert }^2/\left(n-p_{m^{\ast }}\right)$ ，$m^{\ast }=\arg {\max }_{1\le m\le M}{p}_m$ 。基于损失函数$L_k^{\ast }\left(w\right)$ ，计算第$k$ 个局部计算机上的最优权重${\widehat{w}}_k^{\ast }=\arg {\min }_{w\in \mathcal{X}}{L}_k^{\ast }\left(w\right)$ ，并将$K$ 个最优权重传输给中心计算机，进而通过简单加权平均的方式得到权重的全局估计量${\widehat{w}}^{\ast }=K^{-1}{\displaystyle {\sum }_{k=1}^K{\widehat{w}}_k^{\ast }}$ 。新观测的预测向量$x$ 的最终预测值为

${\widehat{\mu }}_{iter}\left({\widehat{w}}^{\ast }\right)={\displaystyle \sum _{m=1}^M{\widehat{w}}_m^{\ast }}{\widehat{\mu }}_{iter,m}={\displaystyle \sum _{m=1}^M{\widehat{w}}_m^{\ast }}{x}_m^{\text{T}}{\widehat{\beta }}_{iter,m}$。

4. 数值模拟

本节通过数值模拟评估所提出的分布式模型平均方法在有限样本下的表现，并与基于全部数据的模型平均方法进行比较。具体地，考虑如下的四种模型平均方法：1) FULL，基于全部数据的模型平均方法；2) OS，单轮型(OS)分布式模型平均；3) GEI1，基于局部加权最小二乘估计初始化的梯度增强迭代(GEI)分布式模型平均；4) GEI2，基于OS全局估计初始化的梯度增强迭代(GEI)分布式模型平均。

数值模拟数据的生成过程如下：首先，独立生成$N$ 个$p$ 维预测向量$x_i={\left(x_{i1},x_{i2},\cdots ,x_{ip}\right)}^{\text{T}}$ ，其中$p=50$ ，$x_{i1}=1$ ，$x_{ij}$ ，$j=2,\cdots ,p$ ，是独立从标准正态分布生成的随机数；其次，基于下述模型生成生存时间，$\log \left(T_i\right)=x_i^{\text{T}}\beta +{\epsilon }_i$ ，$i=1,2,\cdots ,N$ 。其中$\beta ={\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_p\right)}^{\text{T}}$ ，${\beta }_j=1.5j^{-1}$ ，$j=1,\cdots ,p$ ，误差${\epsilon }_i$ 独立同分布于正态分布$N\left(0,{\sigma }^2\right)$ ；最后，令$C_0$ 是某个正常数，从$\left(0,C_0\right)$ 上的均匀分布生成$N$ 个对数删失时间$C_i$ ，进而得到$Z_i=\min \left\{\log \left(T_i\right),C_i\right\}$ 和删失指示变量${\delta }_i=I\left(Y_i\le C_i\right)$ 。在数据生成过程中，选择不同的${\sigma }^2$ ，使得决定系数$R^2=\mathrm{var}\left(x_i^{\text{T}}\beta \right)/\left[\mathrm{var}\left(x_i^{\text{T}}\beta \right)+{\sigma }^2\right]$ 分别取到0.1，0.2，0.3，0.4，0.5，0.6，0.7，0.8，0.9；选择不同的常数$C_0$ ，使得在不同的情形下删失率达到40%和60%。

在预测向量$x_i={\left(x_{i1},x_{i2},\cdots ,x_{ip}\right)}^{\text{T}}$ 中，我们选取其前$p_1$ 个分量进行建模。具体地，构造$M=p_1$ 个嵌套的子模型，其中第$m$ 个候选模型的协变量由$x_i$ 的前$m$ 个分量构成。注意到，第1个候选模型仅包含截距项；第$M$ 个候选模型包含全部选取的$p_1$ 个分量，也就是候选模型中最大的模型。基于$H=100$ 次的数据重复，我们汇报每个方法平均计算时间(time，单位：秒)和平均均方预测误差(mean squared prediction error,

MSPE) $\text{MSPE}=H^{-1}{\displaystyle {\sum }_{h=1}^H{N}^{-1}{\displaystyle {\sum }_{i=1}^N{\left({\mu }_i-{\widehat{\mu }}_i\right)}^2}}$，其中${\mu }_i$ 和${\widehat{\mu }}_i$ 分别是第$i$ 个个体真实的和估计的条件均值。

(a) 删失率40% (b) 删失率60%

Figure 1. Comparison of MSPE ratios for various methods versus the FULL method under different values of $R^2$

图1. 不同$R^2$ 下各方法与FULL方法的MSPE比值的比较

首先，基于不同决定系数$R^2$ 的设定，比较本文所提方法与FULL方法的有限样本差异，见图1。设定$N=2^{13}$ ，$K=4$ ，$p_1=20$ ，$T=6$ ，图1展示了不同$R^2$ 下本文所提方法的MSPE和FULL方法的MSPE的比值。比值曲线越接近于1表明所提方法的有限样本效果越好。从图1的结果可以看出，在各种$R^2$ 的取值下，本文提出方法都展现了良好的有限样本性质，三种方法的MSPE都非常接近FULL方法的MSPE。与此同时，梯度增强迭代型方法在多数情况下略优于单轮型方法，但并未呈现出稳定且一致的优势。

其次，进一步比较本文所提方法与FULL方法在不同样本量下的有限样本差异，见图2。设定$K=4$ ，$p_1=20$ ，$T=6$ ，$R^2=0.5$ ，图2展示了不同样本量$N$ 下OS、GEI1和GEI2的MSPE和FULL方法的MSPE的比值。从图2的结果可以看出，在各种$N$ 的取值下，本文提出方法都展现了良好的有限样本性质，三种方法的MSPE都非常接近FULL方法的MSPE。

(a) 删失率40% (b) 删失率60%

Figure 2. Comparison of MSPE ratios of various methods relative to the FULL method under different sample sizes

图2. 不同样本量下各方法与FULL方法的MSPE比值的比较

Table 1. Comparison of average computation time for different methods (Unit: seconds)

表1. 各方法平均计算耗时的比较(单位：秒)

$\left(N,K,R^2,p_1\right)$	FULL	OS	GEI1	GEI2
$\left(2^{12},2,0.5,20\right)$	2.366	1.208	1.182	1.855
$\left(2^{13},2,0.5,20\right)$	5.542	2.804	2.773	2.757
$\left(2^{13},4,0.5,20\right)$	5.725	1.418	1.424	1.381
$\left(2^{13},4,0.9,20\right)$	5.714	1.415	1.419	1.379
$\left(2^{13},4,0.9,30\right)$	12.003	3.062	2.996	3.012

最后，本文结合多个不同的参数设置，综合评估各方法的计算耗时，见表1。表1给出了各方法在一些情形下平均计算所需的时间。由表1结果可见，本文所提方法所需的平均计算时间显著低于基于完整数据的模型平均方法。并且，这种计算上的优势会随着样本量和备选模型规模的增大进一步凸显。此外，梯度增强迭代型方法和单轮型方法的计算效率相近，无明显差异。

5. 实例分析

美国国家健康与营养调查(National Health and Nutrition Examination Survey, NHANES)是由美国疾病控制与预防中心下属的国家健康统计中心(National Center for Health Statistics, NCHS)主导的一项旨在评估美国成人和儿童健康与营养状况的权威调查项目。该调查始于20世纪60年代初期，自1999年以来，NCHS持续开展NHANES，并搭建了相应的NHANES数据库。NHANES收集的数据包括参与者的人口学特征、体检指标、重要的生化指标以及死亡随访信息等。本文聚焦于NHANES的2009年至2018年的五个连续调查周期数据，评估所提方法在实际数据分析中的预测能力和计算效率，并与基于完整数据的方法进行比较。

通过对数据进行预处理，我们得到一个包含25313个观测值的数据集，其删失率为92%。除生存时间外，每个观测值包含20个预测变量，按与生存时间的皮尔逊相关系数的绝对值从大到小排序如下：红细胞分布宽度、平均红细胞血红蛋白浓度、平均红细胞血红蛋白含量、年龄、单核细胞计数、收缩压第三次读数、平均收缩压、家庭受访者年龄、收缩压第二次读数、最大充气水平、单核细胞百分比、嗜碱性粒细胞绝对值、平均血小板体积、嗜碱性粒细胞百分比、骨密度状态、种族、总胆固醇、胆固醇、糖化血红蛋白、红细胞计数。基于这些预测变量，类似于数值模拟研究，构建20个嵌套的备选模型用于模型平均。

Table 2. Comparison of four methods under NHANES data

表2. NHANES数据下四种方法的比较

		FULL	OS	GEI1	GEI2
$K=2$	MSPE	0.143	0.148	0.181	0.169
	Times (s)	0.719	0.358	0.348	0.350
$K=3$	MSPE	0.143	0.152	0.267	0.197
	Times (s)	0.719	0.242	0.237	0.238

在方法评估阶段，首先将数据集随机划分为训练集和测试集，其中训练集和测试集的样本量分别为24,313和1000；其次，在训练集上进行模型参数和权重的估计，然后在测试集上进行预测和评估；最后，把这个过程重复100次，计算平均的平方预测误差和计算耗时。表2给出了实际数据分析得到的平均平方预测误差和计算时间。从表2结果可见，本文所提方法的计算效率明显高于基于完整数据的方法；OS方法的MSPE与FULL方法基本一致，GEI方法的MSPE稍高于FULL方法。在梯度增强迭代分布式模型平均方法中，当分组数$K=3$ 时的MSPE要高于$K=2$ 时。究其原因，主要在于数据集的删失率较高，随着分组数的增加，各局部子集的样本量急剧缩减，致使逆概率删失加权估计的不稳定性显著加剧。

6. 结论与展望

针对大规模分布式生存数据，在加速失效模型下，本文提出了两种Mallows类型的分布式模型平均方法：基于单轮型估计的分布式模型平均和基于梯度增强迭代估计的分布式模型平均。本文不仅通过数值模拟验证了所提方法在有限样本下的表现，还进一步利用实际数据检验了其实用性。数值模拟与实际数据验证结果表明，本文提出的两种分布式模型平均方法均具有优良的预测性能。

在估计删失时间的分布函数时，本文采用了全局Kaplan-Meier估计的策略，旨在兼顾估计效率与计算可行性。具体而言，该策略能够基于全体样本提高估计的精度，且其传输向量的计算负担较轻，其所带来的时间成本是可以接受的。然而，在数据分布式存储且无法共享的场景下，全局估计策略将不再适用，需要转而寻求能充分利用分布式特性的估计方法，例如利用局部数据直接估计删失时间的分布，或采用分布式迭代算法等。此外，本文使用了最基本的Kaplan-Meier估计。事实上，通过引入更精细的统计手段，如核光滑技术，可以构造出具有更优良性质的删失时间的分布式估计。

据作者了解，现有针对生存数据的分布式模型平均的研究仍较为匮乏。本文仅考虑了生存时间右删失的情形，然而在实际医学、公共卫生等应用中，生存时间往往伴随更复杂的删失机制，例如信息删失、区间删失等。此外，Cox比例风险率模型是生存分析研究中应用最广泛的模型，但其框架下的分布式模型平均尚未见相关研究。诸如此类的科学问题还很多，因此，生存数据的分布式模型平均还存在诸多亟待解决的科学问题，具备广阔的研究前景。

基金项目

教育部人文社会科学研究规划基金项目(21YJA910002)，国家级大学生创新训练计划项目(202510446018)，国家社会科学基金一般项目(25BTJ038)。

NOTES

^*通讯作者。