自适应设计下基于COX回归模型的序贯压缩估计研究

Sequential Shrinkage Estimate Based on COX Regression Model under Adaptive Design

摘要: 在COX回归模型的应用中，我们经常会遇到包含太多变量的数据集，而这些变量中只有少数变量对模型有贡献。因此，在推断过程中估计“无效”变量会浪费大量的样本。在本文中，我们提出一种基于自适应压缩估计的序贯抽样策略来构造“有效”参数的固定长度的置信集，这样在忽略模型中的“无效”变量影响的同时，使用最少样本将模型中的“有效”变量快速地识别出来。最后，在自适应设计下对我们所提出的序贯抽样策略进行数值模拟并且数值模拟达到了预期的效果。

Abstract: In the applications of COX regression models, we always encounter the data sets which contain too many variables that only a few of them contribute to the model. Therefore, it will waste much more samples to estimate the “non-effective” variables in the inference. In this paper, we use a sequential procedure for constructing the fixed size confidence set for the “effective” parameters to the model based on an adaptive shrinkage estimate such that the “effective” coefficients can be efficiently identified with the minimum sample size. Adaptive design is considered for numerical simulation.

1. 引言

COX比例风险模型 [1] 是一种常用的处理具有协变量的生存数据分析方法。它被广泛应用于生物医学研究和工程等许多领域，用来评估协变量效应。在本文中，我们设失效时间为T，条件风险函数为

$h\left(t|x\right)=h_0\left(t\right){\text{e}}^{{\beta }^{\prime }X}$ (1)

其中， $\beta$ ( $p$ 维向量)是回归系数， $h_0\left(t\right)$ (非负函数)称为基准风险函数。Cox (1984) [2] 和Fleming (1991) [3] 对COX比例风险模型做了大量的研究。然而，在生物学、工程学和流行病学等研究中，数据集通常有大量的解释变量，但其中只有少数对模型有贡献。也就是说，在一个 $p$ 维的回归系数中只有 $p_0$ ( $p_0<p$ 且 $p_0$ 未知)个分量是取非零值的。这些取非零值的变量被Wang和Zhang (2013) [4] 称为有效变量。目前有很多方法可以用来识别有效变量，如LASSO [5] 和LARS [6] 等。但另外需要关注的问题是，用多少样本才能既识别出有效变量，同时又能使参数估计达到预定的精度。这对于生物学和流行病学等需要考虑抽样成本的研究具有重要的意义。对于线性回归模型，Wang和Zhang (2013) [4] 提出了一种序贯压缩估计方法来识别有效变量，从而达到参数估计的精度。数值模拟结果表明，与传统的序贯抽样方法相比，序贯压缩估计不仅可以从所有变量中识别出有效变量，而且可以节省大量样本。对于COX回归模型，还没有人提出相应的序贯方法。

本文针对COX回归模型提出了一种基于自适应压缩估计(ASE)来构造有效变量的固定窗宽的置信集的序贯抽样方法，使有效变量能以最小样本量快速识别。本文将在适应性设计(adaptive design)下研究所提出的自适应压缩估计(ASE)的大样本性质，同时在适应性设计下通过数值模拟得到了很好的模拟结果。

2. 序贯自适应压缩估计(ASE)

2.1. 最大偏似然估计(MPLE)

设样本容量为 $n$， $T_i$ 和 $C_i$ 分别是第 $i\left(i=1,2,\cdots ,n\right)$ 个随机样本的失效时间和删失时间。假定 $p$ 维向量 $X_i={\left(X_{i1},X_{i2},\cdots ,X_{in}\right)}^{\text{T}}$ 为第 $i$ 个样本的协变量并且 $T_i$ 和 $C_i$ 关于 $X_i$ 条件独立。令 $Y_i=\min \left\{T_i,C_i\right\}$ 是我们真正能观测到的个体失效时间， ${\delta }_i=I\left\{T_i\le C_i\right\}$ 是第 $i$ 个样本被观测到的示性函数。假设 $t_1<t_2<\cdots <t_n$， $R_j$ 是 $t_j$ 时刻的风险集，即 $R_j=\left\{i:Y_i\ge t_j\right\}$。这样我们就可以定义模型(1)的偏似然函数为

${\displaystyle \underset{j=1}{\overset{N}{\prod }}\frac{\exp \left({\beta }^{\text{T}}{X}_{\left(j\right)}\right)}{{\displaystyle {\sum }_{i\in R_j}\exp \left({\beta }^{\text{T}}{X}_{(\; i\; )}\right)}}}$

由此可得其对数似然函数为

$L\left(\beta \right)={\displaystyle \underset{j=1}{\overset{N}{\sum }}\left\{{\beta }^{\text{T}}{X}_{\left(j\right)}-\log \left[{\displaystyle \underset{i\in R_j}{\sum }\exp \left({\beta }^{\text{T}}{X}_i\right)}\right]\right\}}$

设回归参数 $\beta$ 的最大偏似然估计(MPLE)为 $\tilde{\beta }$，则 $\tilde{\beta }$ 通可以通过解下面的得分方程获得： $U\left(\beta \right)=0$，其中

$U\left(\beta \right)=\frac{\partial L\left(\beta \right)}{\partial \beta }={\displaystyle \underset{j=1}{\overset{N}{\sum }}\left\{X_{\left(j\right)}-\frac{{\displaystyle {\sum }_{i\in R_j}{X}_i\exp \left({\beta }^{\text{T}}{X}_i\right)}}{{\displaystyle {\sum }_{i\in R_j}\exp \left({\beta }^{\text{T}}{X}_i\right)}}\right\}}$

2.2. 自适应压缩估计(ASE)

设 $\kappa =\kappa (n)$，当 $n\to \infty$ 时，存在 $0<\delta <1/2$ 和 $\gamma >0$ 使得 $n^{\frac{1}{2}}\kappa \to 0$， $n^{\frac{1}{2}+\gamma \delta }\kappa \to \infty$。本文中我们需要以下假定：

(A1). 解释变量 $x_i$ 满足条件 ${\sup }_i\Vert x_i\Vert <\infty$，并且残差项 ${\epsilon }_i=\hat{\Lambda }\left(Y_i\right)\exp \left({\hat{\beta }}^T{X}_i\right)$ 具有有限二阶矩，即当 $\zeta >2$ 时， $E{\left|{\epsilon }_i\right|}^{\zeta }<\infty$.其中是 $\hat{\Lambda }$ 累计基准风险函数；

(A2). $\underset{n\to \infty }{\lim }{I}_n\left(\beta \right)/n=\Sigma$，其中 $I_n\left(\beta \right)$ 是回归参数 $\beta$ 的信息矩阵， $\Sigma$ 是一个正定矩阵。

定义2.2.1 设 $\tilde{\beta }$ 为模型(1)的最大偏似然估计，则称 ${\hat{\beta }}_n=In\left(\epsilon \right){\tilde{\beta }}_n$ 为回归系数 ${\beta }_0$ 的自适应压缩估计(ASE)，其中 $I_n\left(\epsilon \right)=diag\left\{I_{n1}\left(\epsilon \right),I_{n2}\left(\epsilon \right),\cdots ,I_{np}\left(\epsilon \right)\right\}$ 是一个 $p\times p$ 维对角阵。同时可以证明 ${\hat{\beta }}_n=In\left(\epsilon \right){\tilde{\beta }}_n$ 满足相合性和渐进正态性。

2.3. 序贯抽样策略

由文献 [7]， [8] 的结论可以证明 $\sqrt{n}\left({\hat{\beta }}_n-{\beta }_0\right),n=1,2,\cdots$ 是依概率一致连续的，由此可得如下定理：

定理2.3.1设随机变量N(t)取正整数值，当 $t\to \infty$ 有 $N\left(t\right)/t$ 依概率收敛于1，且条件(A1)和(A2)成立，则当 $t\to \infty$ 时，

$\sqrt{N\left(t\right)}\left({\hat{\beta }}_{N(t)}-{\beta }_0\right)\to N\left(0,I_0\Sigma I_0{}^{-1}\right)$

由定理2.3.1我们可以构造 ${\beta }_0$ 的置信集和能够决定最小样本量的停止法则的序贯抽样策略。设 $\left\{\left(y_i,x_i\right):i=1,2,\cdots ,k\right\}$ 是最先进入研究的 $k$ 个样本，用 $C_k$ 来表示。在任意给定小正数 $\epsilon$ 下，

${\hat{p}}_0\left(k\right)={\displaystyle {\sum }_{j=1}^p{I}_{kj}(\; \epsilon \; )}$

是回归系数 $p_0$ 基于条件 $C_k$ 的估计量。令 $a_k^2\in R$ 对任意 $\alpha >0$，有 $P\left({\chi }_{{\hat{p}}_0(k)}^2\le a_k^2|C_k\right)=1-\alpha$ 成立。现在定义停时法则 $N_d$ 为

$N=N_d\equiv \inf \left\{k:k\ge n_0\text{and}\frac{d^2}{a_k^2}\ge {\nu }_k\right\}$, (2)

其中 ${\nu }_k$ 是 $k{I}_k\left(\epsilon \right){\left(\Sigma \right)}^{-1}{I}_k\left(\epsilon \right)$ 的最大特征值， $d$ 是置信集的预设精度。在本文的序贯估计策略中，一次只有一个新的观测进入研究直到满足(2)式的停止法则时就停止抽样，此时 ${\beta }_0$ 的置信集为

$R_N=\left\{Z\in R^p:\frac{S_N}{N}\le \frac{d^2}{{\nu }_N}且当I_{N_j}\left(\epsilon \right)=0时,z_j=0,1\le j\le p\right\}$ (3)

其中 $S_N={\left(Z_{N_1}-{\hat{\beta }}_{N_1}\right)}^{\text{T}}{\tilde{\Sigma }}_{11}\left(Z_{N_1}-{\hat{\beta }}_{N_1}\right)$。我们所提出的序贯抽样方法致力于找到有效变量的同时忽略无效变量的影响，这是和传统序贯方法相比我们能够节省大量样本的关键，在下面的定理中我们给出停时 $N_d$ 和置信集 $R_N$ 的相关性质。

定理2.3. 2假定条件(A1)和(A2)都成立，设N是满足(2)式的停时，则：

(i). $\underset{d\to 0}{\lim }\frac{d^{2}N}{a^2\nu }=1$，a.s.；(ii). $\underset{d\to 0}{\lim }P\left({\beta }_0\in R_N\right)=1-\alpha$ ；

(iii). $\underset{d\to 0}{\lim }\frac{d^{2}E\left(N\right)}{a^2\nu }=1$ ； (iv). $\underset{d\to 0}{\lim }{\hat{p}}_0\left(N\right)=p_0$，a.s.且 $\underset{d\to 0}{\lim }E\left({\hat{p}}_0\left(N\right)\right)=p_0$，

其中 $\nu$ 是矩阵 $I_0{\Sigma }^{-1}{I}_0$ 的最大特征值。

3. 数值模拟

我们在固定样本量下用所提方法对随机数据集合进行分析，以此来验证所提出的序贯压缩估计方法的性能。按照停止法则的定义，当抽样停止时，最终的置信集将满足预设精度和覆盖概率，因此我们可以比较分别基于MPLE和ASE的序贯抽样方法的平均停时。由于序贯压缩估计方法忽略无效变量的影响，故理论上平均所需停时应该显著小于不考虑变量选择的序贯方法。如果事先已知有效变量为 $p_0$ 个同时无无效变量，那么只使用这 $p_0$ 个有效变量的序贯方法无疑是效率最高的。所以，为便于比较，我们将所有( $p_0$ 个)变量全部为有效变量的序贯估计方法作为基准线，在此情况下所获得的样本量应该是最小的。

在自适应设计下，随机模拟数据集中的 $x_1$ 仍然由多元标准正态分布生成， $x_j\left(j>1\right)$ 由均值为 ${\displaystyle \underset{i=1}{\overset{j-1}{\sum }}\left[x_j/\left(j-1\right)\right]}$，方差协方差矩阵为单位阵的多元正态分布生成。不失一般性，选择模型(1)中的基准风险函数为 $h_0\left(t\right)=t^2$。回归系数真值取 $\left(-1.2,2.0,0,0,0,0,0,0,0,0\right)$，其中含有八个无效变量，回归系数置信集的预设精度 $d\in \left\{0.3,0.4,0.5,0.6\right\}$，取 $\alpha =0.05$， $\gamma =1$， $\delta =0.45$， $\theta =0.75$。另外当用ASE方法时我们用BIC方法来确定 $\epsilon$，

$BIC=-2\left({\displaystyle \underset{j=1}{\overset{N}{\sum }}\left({\beta }^{\text{T}}{X}_{\left(j\right)}-\log \left({\displaystyle \underset{j\in R_j}{\sum }\exp \left({\beta }^{\text{T}}{X}_i\right)}\right)\right)}\right)+\log \left(n\right)\times df/n$

其中 $df$ 表示 $\beta$ 中非零分量的个数。

Table 1. Results of sequential sampling method based on ASE, MPLE with all variables and MPLE with only p 0 non-zero variables for COX regression model

表1. COX回归模型下分别应用ASE，MPLE和 ${\text{MPLE}}_{p_0}$ 的序贯抽样方法的结果分析

${\kappa }^{*}=d^{2}N/\left(a^2\nu \right)$ ； $C{P}^{+}$ 是95%置信集 $R_N$ 的经验覆盖概率；^**经验标准差在括号内。

表1描述了COX回归模型下的序贯抽样方法的数值模拟结果。在表1中我们列出了最终样本量N (停

时)， ${\kappa }^{*}=d^{2}N/\left(a^2\nu \right)$ 和95%置信集的经验覆盖概率 $R_N$。所有三种情况( ${\text{MPLE}}_{p_0}$, ASE, MPLE)下的 $\kappa$ 值都

非常接近1，并且当d不断减小时经验覆盖概率CP越来越接近95%，正如定理2.3.2描述的一样。然而，应用MPLE方法所得的样本量N比应用ASE方法和 ${\text{MPLE}}_{p_0}$ 都大得多。而应用ASE的抽样策略所需的样本量和应用 ${\text{MPLE}}_{p_0}$ 的抽样策略所需样本量差不多，这说明我们所提方法在变量选择的同时效率和回归参数中只有有效变量无无效变量的情况下的效率非常接近，而比不做变量选择情况下(即MPLE)的抽样效率提高很多。

表2比较了在估计COX回归模型的回归系数时分别应用ASE和MPLE的抽样策略对识别回归系数中的有效变量和无效变量的效率。从结果可以看出应用ASE的抽样策略时不能被正确识别的零变量的平均个数几乎趋向于0，而能被正确识别的非0变量的平均个数和模型中有效变量个数的真值非常接近(2和8)。结果表明基于ASE的序贯抽样策略下 ${\hat{p}}_0$ 是 $p_0$ 的优良估计。而基MPLE的序贯抽样策略不能识别有效变量，因此无法获得 $N_c^{*}$ 和 $N_{ic}^{*}$ 的值。此外，所有参数的估计值和它们的真值都非常接近。

Table 2. Power of variable identification and estimation of nonzero components under sequential sampling method based on ASE and MPLE with COX regression model

表2. COX回归模型下分别应用ASE和MPLE的序贯抽样策略的变量识别和非零参数估计效率

$N_{ic}^{*}$ ： $\beta$ 中零分量(无效变量)被错误识别的平均个数； $N_c^{*}$ ： $\beta$ 中非零分量(有效变量)被正确识别的平均个数。

4. 结论

在COX回归模型下基于自适应压缩估计(ASE)建立的序贯抽样方法不仅能够用最少的样本识别出回归参数中的有效变量，同时可以使回归参数的估计值达到预设的精度。数值模拟结果表明和传统的序贯抽样方法相比，我们提出的方法能够节省大量样本。然而，本文中所提方法涉及到的变量维数是固定的，后期我们将研究当变量维数随样本量变化时的序贯抽样方法的相关性质。

基金项目

1) 新疆师范大学博士科研启动基金项目：“基于广义线性模型的序贯分析研究”XJNUBS1539。

2) 新疆维吾尔自治区高校科研计划项目：“基于COX比例风险回归模型的序贯分析研究”(XJEDU2016I033)。