在过去10年里分子生物学试验技术的进展给我们带来了丰富的生物医学数据，举例来说，DNA显微序列可以用来测量一个细胞中成千上万的基因。这种类型的数据中样本维度p比样本量n要大得多，对于传统的统计推断方法来说是一个巨大的挑战，有很多经典的推断方法在这种情况下变得不适用。这种情形下有效的变量选择方法就变得尤为重要。比较著名的高维数据变量选择方法有Lasso^[1]，SCAD^[2]和MCP^[3]等。

当研究关于患者生存状态的医疗数据时，将高维的生物医疗数据和患者的生存状态数据结合起来分析是一个很有效的方法。因此近些年来也有很多关于高维生存分析模型的变量选择方法，比如Bradic等^[4]关于高维Cox模型的正则化方法，Gorst-Rasmussen和Scheike^[5]关于高维单指数模型的筛选方法，Lin和Lyu^[6]关于高维可加模型的正则化方法等等。高维生存分析模型还广泛地应用到信用风险分析，比如Fan等^[7]。

由于在实际生活中，我们经常会遇到带有测量误差的数据，所以对于带有测量误差数据的分析方法也是一个重要的研究方向，对于高维线性模型有Loh和Wainwright^[8]以及Datta和Zou^[9]的相关工作；对于变系数模型，有刘智凡等^[10]的工作。对于带有测量误差的生存分析数据的变量选择方法，代表文章有Song和Wang^[11]关于工具变量的工作，Chen和Yi^[12]关于Cox模型左截断右删失数据的工作。高维生存分析模型由于其计算复杂度较高以及理论性质较为复杂，所以对于带有测量误差的高维生存分析数据的工作随着近些年大数据的迅速发展才逐步出现在视野之中。具有代表性的文章有Chen和Yi^[13]关于高维生存分析图模型的工作以及Chen等^[14]关于高维Cox模型利用纠正似然函数的工作。本文选择同样具有重要应用的可加风险模型作为基础，结合处理高维线性模型的正则化方法对带有测量误差的生存分析数据进行分析。

1 研究背景

本文所采用的模型为高维可加风险模型，结合高维线性模型测量误差处理办法对带有测量误差的生存分析数据进行分析。下面对高维可加风险模型和高维线性模型测量误差处理方法分别进行介绍。

1.1 高维可加风险模型

对于生存分析数据的变量选择技术的发展已经不拘泥于Cox模型，可加风险模型便是除Cox模型以外的一种重要替代方式。可加风险模型假设失效时间为T的风险函数和p维的协变量X(·)有如下形式的关系

$ \lambda(t \mid \boldsymbol{X})=\lambda_0(t)+\boldsymbol{\beta}_0^{\mathrm{T}} \boldsymbol{X}(t), $

(1)

其中：λ₀(·)是一个不确定的基线风险函数，β₀是一个p维的回归系数。令C为删失时间，则定义删失失效时间为CFT=C∧T，令CFT=t₁, …, t_n，失效指数定义为δ=I(T≤C)，其中I(·)为指示函数，令X(t)=(X₁(t), …, X_p(t))并且假设给定X观察到的数据为(CFT, δ, X(·))，风险函数由式(1)给出。

采用常用的计数手段，定义观察到的失效计数序列为N_i(t)=I(t_i≤t, δ_i=1)，风险中指数为Y_i(t)=I(t_i≥t)，计数过程鞅为

$ M_i(t)=N_i(t)-\int_0^t Y_i(s) \lambda_0(s)+\boldsymbol{\beta}_0^{\mathrm{T}} \boldsymbol{X}_i(s) \mathrm{d} s, $

(2)

后文也将用N(t), Y(t)和M(t)来代表这些计数过程的广义形式。

Lin和Ying^[15]采用一种有如下形式的伪得分方程来对可加风险模型进行分析：

$ \begin{gathered} \boldsymbol{U}_0(\boldsymbol{\beta})=\frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau \boldsymbol{X}_i(t)-\overline{\boldsymbol{X}}(t) \\ \left\{\mathrm{d} N_i(t)-Y_i(t) \boldsymbol{\beta}^{\mathrm{T}} \boldsymbol{X}_i(t) \mathrm{d} t\right\}, \end{gathered} $

(3)

其中$ \boldsymbol{\beta} \in \mathbb{R}^p$，并且

$ \overline{\boldsymbol{X}}(t)=\sum\nolimits_{j=1}^n Y_j(t) \boldsymbol{X}_j(t) / \sum\nolimits_{j=1}^n Y_j(t), $

(4)

τ是最大的跟踪时间(生存时间和删失时间的最大值)。这个估计函数关于回归系数是线性的，令

$ \boldsymbol{b}_0=\frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau \boldsymbol{X}_i(t)-\overline{\boldsymbol{X}}(t) \mathrm{d} N_i(t), $

(5)

和

$ \boldsymbol{V}_0=\frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau Y_i(t)\left\{\boldsymbol{X}_i(t)-\overline{\boldsymbol{X}}(t)\right\}^{\otimes 2} \mathrm{~d} t, $

(6)

其中v^⊗2=vv^T，通过一些代数变换，可以写出如下等式

$ \boldsymbol{U}_0(\boldsymbol{\beta})=\boldsymbol{b}_0-\boldsymbol{V}_0 \boldsymbol{\beta}. $

(7)

在没有测量误差的情况下，V₀是半正定的，式(7)两边关于β积分就可以得到损失函数

$ L(\boldsymbol{\beta})=\frac{1}{2} \boldsymbol{\beta}^{\mathrm{T}} \boldsymbol{V}_0 \boldsymbol{\beta}-\boldsymbol{b}_0{ }^{\mathrm{T}} \boldsymbol{\beta}, $

(8)

Leng和Ma^[16]以及Martinussen和Scheike^[17]都建议用上述损失函数配合正则化方法对可加风险模型(1)进行变量选择。本文的相关工作也是在此基础上进行。

1.2 高维线性模型测量误差数据的处理方法

为了进一步构建更深层次的讨论，假设观察到的是被污染的协变量矩阵

$ \boldsymbol{Z}(\cdot)=\left(z_{i j}(\cdot)\right)_{1 \leqslant i \leqslant n, 1 \leqslant j \leqslant p} \text {, } $

(9)

而不是真实的协变量矩阵X(·)。有很多种造成测量误差的途径，在加法测量误差设定中，z_{i, j}(·)=x_{i, j}(·)+a_{i, j}，其中A(·)=(a_{i, j})是加法测量误差。在乘法测量误差设定中，z_{i, j}(·)=x_{i, j}(·)m_{i, j}，其中m_{i, j}就是乘法测量误差。缺失数据可以看作乘法测量误差的一个特殊形式，m_{i, j}=I(x_{i, j}(·)没缺失)。

不失一般性，用Lasso算法来举例说明测量误差的影响，对于线性模型y=Xβ+$ \boldsymbol\epsilon$来说，Lasso算法是最小化

$ 1 /(2 n)\|\boldsymbol{y}-\boldsymbol{X} \boldsymbol{\beta}\|_2^2+\lambda\|\boldsymbol{\beta}\|_1, $

(10)

这等价于最小化

$ \frac{1}{2} \boldsymbol{\beta}^{\mathrm{T}} \boldsymbol{\varSigma} \boldsymbol{\beta}-\boldsymbol{\rho}^{\mathrm{T}} \boldsymbol{\beta}+\lambda\|\boldsymbol{\beta}\|_1, $

(11)

其中: $\boldsymbol{\varSigma}=\frac{1}{n} \boldsymbol{X}^{\mathrm{T}} \boldsymbol{X}, \boldsymbol{\rho}=\frac{1}{n} \boldsymbol{X}^{\mathrm{T}} \boldsymbol{y} $，λ为Lasso算法的惩罚参数。假设测量误差是加法误差且服从均值是0、方差是τ²的正态分布，τ是已知的常数。Loh和Wainwright^[8]建议使用无偏估计

$ \hat{\boldsymbol{\varSigma}}=\frac{1}{n} \boldsymbol{Z}^{\mathrm{T}} \boldsymbol{Z}-\tau^2 I, \tilde{\boldsymbol{\rho}}=\frac{1}{n} \boldsymbol{Z}^{\mathrm{T}} \boldsymbol{y}, $

(12)

然后解决下面的优化问题来得到β的估计：

$ \frac{1}{2} \boldsymbol{\beta}^{\mathrm{T}} \hat{\boldsymbol{\varSigma}} \boldsymbol{\beta}-\tilde{\boldsymbol{\rho}}^{\mathrm{T}} \boldsymbol{\beta}+\lambda\|\boldsymbol{\beta}\|_1, $

(13)

然而$ \hat{\boldsymbol{\varSigma}}$在有测量误差的情况下经常会出现负的特征值(高维情形下更加常见)，给优化上述问题带来了困难。为解决这个问题，Loh和Wainwright^[8]采用如下方法(简记为NCL):

$ \hat{\boldsymbol{\beta}} \in \underset{\|\boldsymbol{\beta}\|_1 \leqslant R}{\operatorname{argmin}} \frac{1}{2} \boldsymbol{\beta}^{\mathrm{T}} \hat{\boldsymbol{\varSigma}} \boldsymbol{\beta}-\tilde{\boldsymbol{\rho}}^{\mathrm{T}} \boldsymbol{\beta}+\lambda\|\boldsymbol{\beta}\|_1 . $

(14)

其中R是一个跟稀疏度有关的常数。Datta和Zou^[9]提出一种最近邻正定投影矩阵的算法来解决上述问题，对于任意方阵K:

$ (\boldsymbol{K})_{+}=\underset{\boldsymbol{K}_1 \geqslant 0}{\operatorname{argmin}}\left\|\boldsymbol{K}-\boldsymbol{K}_1\right\|_{\max } . $

(15)

其中$ \|\boldsymbol{K}\|_{\max }=\max \limits_{i, j}\left|K_{i, j}\right|$，本文在附录A中给出了最近邻正定投影矩阵的详细计算方法。定义$\tilde{\boldsymbol{\varSigma}}=(\hat{\boldsymbol{\varSigma}})_{+} $，则它们的CoCoLasso估计量为

$ \hat{\boldsymbol{\beta}}=\underset{\boldsymbol{\beta}}{\operatorname{argmin}} \frac{1}{2} \boldsymbol{\beta}^{\mathrm{T}} \tilde{\boldsymbol{\varSigma}} \boldsymbol{\beta}-\tilde{\boldsymbol{\rho}}^{\mathrm{T}} \boldsymbol{\beta}+\lambda\|\boldsymbol{\beta}\|_1 . $

(16)

由定义可知，$ \tilde{\boldsymbol{\varSigma}} $一直是半正定矩阵，为解决上述优化问题，可以重写优化函数为

$ \hat{\boldsymbol{\beta}}=\underset{\boldsymbol{\beta}}{\operatorname{argmin}} \frac{1}{n}\|\tilde{\boldsymbol{y}}-\tilde{\boldsymbol{Z}} \boldsymbol{\beta}\|_2^2+\lambda\|\boldsymbol{\beta}\|_1 . $

(17)

其中：$ \tilde{\boldsymbol{Z}}$是$ \tilde{\boldsymbol{\varSigma}} $的Cholesky分解，$ \tilde{\boldsymbol{y}}$是$ \tilde{\boldsymbol{Z}}^{\mathrm{T}} \tilde{\boldsymbol{y}}=\boldsymbol{Z}^{\mathrm{T}} \boldsymbol{y}$的解。式(17)就是Lasso算法的优化函数，我们已经有很多现成的优化算法来求解。

2 带有测量误差的高维可加风险模型的变量选择方法 2.1 简化伪得分方程

在第1节中已经介绍了Lin和Ying^[15]的伪得分方程的具体形式，下面将在协变量X期望值为0的前提下简化该伪得分方程，提出一种全新的更加容易计算且符合实际情况的损失函数。首先定义

$ \boldsymbol{b}=\frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau \boldsymbol{X}_i(t) \mathrm{d} N_i(t), $

(18)

以及

$ \boldsymbol{b}=\frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau \boldsymbol{X}_i(t) \mathrm{d} N_i(t), $

(19)

则有

$ \boldsymbol{V}=\frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau Y_i(t)\left(\boldsymbol{X}_i(t)\right)^{\otimes 2} \mathrm{~d} t, $

(20)

接着定义

$ \begin{aligned} \boldsymbol{b}= & \frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau \boldsymbol{X}_i(t) \mathrm{d} N_i(t) \\ = & \frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau Y_i(t)\left(\boldsymbol{X}_i(t)\right)^{\otimes 2} \mathrm{~d} t+ \\ & \frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau \boldsymbol{X}_i(t) \mathrm{d} M_i(t) . \end{aligned} $

(21)

由于X的期望为0，所以容易得到E(U(β))=0，在如上定义的基础上，类似于式(7)，有

$ \boldsymbol{U}(\boldsymbol{\beta})=\boldsymbol{b}-\boldsymbol{V} \boldsymbol{\beta}, $

(22)

式(22)对β积分即可得到期望为0时的损失函数

$ L(\boldsymbol{\beta})=\frac{1}{2} \boldsymbol{\beta}^{\mathrm{T}} \boldsymbol{V} \boldsymbol{\beta}-\boldsymbol{b}^{\mathrm{T}} \boldsymbol{\beta}, $

(23)

综上所述即为简化版本的损失函数，我们将基于这个损失函数进行变量选择。

2.2 两种测量误差数据的变量选择方法 2.2.1 加法测量误差

假设观测到的设计矩阵Z(·)被加法测量误差污染，即z_{i, j}(·)=x_{i, j}(·)+a_{i, j}，其中A(·)=(a_{i, j})。同时假设A的行是独立同分布的，均值是0，协方差矩阵是Σ_A，次高斯参数是τ²。假设Σ_A是已知的，则V和b的无偏估计分别为

$ \hat{\boldsymbol{V}}_{\mathrm{add}}=\frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau Y_i(t)\left\{\left(\boldsymbol{Z}_i(t)\right)^{\otimes 2}-\boldsymbol{\varSigma}_A\right\} \mathrm{d} t, $

(24)

和

$ \hat{\boldsymbol{b}}_{\text {add }}=\frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau \boldsymbol{Z}_i(t) \mathrm{d} N_i(t). $

(25)

可以看出$\hat{\boldsymbol{V}}_{\text {add }} $有负的特征值会导致最小化损失函数成为非凸优化，基于Datta和Zou^[9]的方法，可以得到相应的凸损失函数

$ \tilde{f}_{\text {add }}(\boldsymbol{\beta})=(1 / 2) \boldsymbol{\beta}^{\mathrm{T}} \tilde{\boldsymbol{V}}_{\mathrm{add}} \boldsymbol{\beta}-\hat{\boldsymbol{b}}_{\mathrm{add}}^{\mathrm{T}} \boldsymbol{\beta}+\lambda\|\boldsymbol{\beta}\|_1, $

(26)

其中$ \tilde{\boldsymbol{V}}_{\text {add }}=\left(\hat{\boldsymbol{V}}_{\text {add }}\right)_{+}$，在此基础上，运用第1节提到的方法应用Lasso求解最小化上述损失函数即可。即$ \hat{\boldsymbol{\beta}}_{\text {add }}=\underset{\boldsymbol{\beta}}{\operatorname{argmin}} \tilde{f}_{\text {add }}(\boldsymbol{\beta})$。

2.2.2 乘法测量误差

假设测量误差是乘法测量误差，观测到的数据为z_{i, j}(·)=x_{i, j}(·)m_{i, j}，其中m_{i, j}就是乘法测量误差。在矩阵写法中，有Z(·)=X(·)⊙M，其中⊙代表矩阵或向量的元素相乘，假设M的行是独立同分布的，均值是μ_M，协方差矩阵为Σ_M，次高斯参数为τ²。假设Σ_M+μ_Mμ_M^T是严格正定的，并且采用无偏估计

$ \begin{gathered} \hat{\boldsymbol{V}}_{\text {mult }}=\frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau Y_i(t) \\ \left\{\left(\boldsymbol{Z}_i(t)\right)^{\otimes 2} / /\left(\boldsymbol{\Sigma}_M+\mu_M \mu_M^{\mathrm{T}}\right)\right\} \mathrm{d} t . \end{gathered} $

(27)

以及

$ \hat{\boldsymbol{b}}_{\text {mult }}=\frac{1}{n} \sum\nolimits_{i=1}^n \int_0^\tau \boldsymbol{Z}_i(t) / / \mu_M \mathrm{~d} N_i(t). $

(28)

其中//代表向量或者矩阵对应元素相除。和加法测量误差模型类似，乘法测量误差下无偏估计矩阵也有可能不是正定的，所以基于Datta和Zou^[9]的方法，可以得到相应的凸损失函数：

$ \tilde{f}_{\text {mult }}(\boldsymbol{\beta})=(1 / 2) \boldsymbol{\beta}^{\mathrm{T}} \tilde{\boldsymbol{V}}_{\text {mult }} \boldsymbol{\beta}-\hat{\boldsymbol{b}}_{\text {mult }}^{\mathrm{T}} \boldsymbol{\beta}+\lambda\|\boldsymbol{\beta}\|_1 . $

(29)

其中$ \tilde{\boldsymbol{V}}_{\text {mult }}=\left(\hat{\boldsymbol{V}}_{\text {mult }}\right)_{+}$。然后运用第1节提到的方法应用Lasso求解最小化上述损失函数即可。即$ \hat{\boldsymbol{\beta}}_{\text {mult }}=\underset{\boldsymbol{\beta}}{\operatorname{argmin}} \tilde{f}_{\text {mult }}(\boldsymbol{\beta})$。对于缺失数据，可以将其视为乘法测量误差，具体办法在第1节已经详细介绍。

3 理论性质

在这一节中给出并推导估计量的l₁和l₂误差界。记我们的估计量为CoCo估计量。首先定义近邻条件：

条件3.1  近邻条件：假设$ \hat{\boldsymbol{V}} \text { 和 } \hat{\boldsymbol{b}}$由一系列参数θ确定, 则存在依赖于β_S，θ和σ²的全局常数C和c，正函数ς和ε₀使得对每个ε≤ε₀，$ \hat{\boldsymbol{V}} \text { 和 } \hat{\boldsymbol{b}}$满足如下概率条件：

$ P\left(\left|\hat{V}_{i j}-V_{i j}\right| \geqslant \varepsilon\right) \leqslant \exp \left(-c n \varepsilon^2 \varsigma^{-1}\right), $

(30)

$ P\left(\left|\hat{b}_j-b_j\right| \geqslant \varepsilon\right) \leqslant \exp \left(-c n^3 \varepsilon^2 \varsigma^{-1}\right) . $

(31)

对所有1≤i, j≤p成立。其中集合S={1, 2, …，s}是回归系数β的支撑集。

同样也需要和线性模型下一样的特征值限制条件：

条件3.2  协方差阵特征值限制条件

$ 0<\Omega=\min\limits _{\boldsymbol{x} \neq 0, \left\|\boldsymbol{x}_{S^c}\right\| \leqslant 3\left\|\boldsymbol{x}_S\right\|_1} \frac{\boldsymbol{x}^{\prime} \boldsymbol{V} \boldsymbol{x}}{\|\boldsymbol{x}\|_2^2} \text {. } $

(32)

条件3.2是一个在高维线性模型变量选择中比较常见的假设。下面给出CoCo估计量的统计误差界：

定理3.1  在式(30)、式(31)和式(32)成立的前提下，对于λ≤min(ε₀, 12ε₀‖β_S‖_∞)和ε≤min(ε₀, Ω/64s)，下式至少以概率

$ 1-p^2 C \exp \left(-c n^3 \lambda^2 \varsigma^{-1}\right)-p^2 C \exp \left(-c n \varepsilon^2 \varsigma^{-1}\right) $

成立：

$ \|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}\|_2 \leqslant C \lambda \sqrt{s} / \varOmega, \|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}\|_1 \leqslant C \lambda s / \varOmega . $

(33)

其中

$ \hat{\boldsymbol{\beta}}=\underset{\boldsymbol{\beta}}{\operatorname{argmin}} \frac{1}{2} \boldsymbol{\beta}^{\mathrm{T}} \tilde{\boldsymbol{V}} \boldsymbol{\beta}-\hat{\boldsymbol{b}}^{\mathrm{T}} \boldsymbol{\beta}+\lambda\|\boldsymbol{\beta}\|_1 . $

(34)

对于两种测量误差$ \tilde{\boldsymbol{V}} \text { 和 } \hat{\boldsymbol{b}}$的计算方式在2.2节中已详细介绍。定理3.1为Datta和Zou^[9]文章中定理1的平行定理，同样也给出了估计量l₁和l₂的误差界。本文在固定设计的前提下进行理论部分的展开，即$ \frac{1}{n}\left\|x_j\right\|_2^2=1$对1≤j≤p成立。且假设观察时间的最大值t_n < ∞，加法测量误差|a_ij|≤L。下面将对加法测量误差和乘法测量误差对条件3.1的满足给出理论上的保证。

引理3.1  $ \hat{\boldsymbol{V}}_{\text {add }} \text { 和 } \hat{\boldsymbol{b}}_{\text {add }}$满足条件3.1中的近邻条件式(30)和式(31)，$ \varsigma=\max \left(\tau^2 t_n^2, \tau^4 t_n^2, L^4\right)$且ε₀=cτ²t_n。其中τ为次高斯参数。

引理3.1说明加法测量误差的计算方法满足近邻条件。下面将对乘法测量误差进行说明。为了保证乘法测量误差的计算方法也满足近邻条件，需要添加额外的正则化条件如下：

$ \max \limits_{i, j}\left|X_{i j}\right|=X_{\max }<\infty, \min \mu_M=\mu_{\min }>0, $

(35)

$ \min \limits_{i, j} E\left(m_1 m_1^{\prime}\right)=M_{\min }>0, \max \mu_M=\mu_{\max }<\infty. $

则接下来有

引理3.2  存在依赖于β_S, θ, σ²和正则化条件式(35)的$ \varsigma$和ε₀使得$ \hat{\boldsymbol{V}}_{\text {mult }} \text { 和 } \hat{\boldsymbol{b}}_{\text {mult }}$满足条件3.1中的近邻条件式(30)和式(31)。

引理3.2说明了乘法测量误差的计算方法满足近邻条件。将引理3.1，引理3.2和定理3.1结合有

推论3.1在条件3.2成立的前提下，定理3.1的结果对加法测量误差估计量$\hat{\boldsymbol{\beta}}_{\text {add }} $和乘法测量误差估计量$\hat{\boldsymbol{\beta}}_{\text {mult }} $成立。

推论3.1给出了加法测量误差估计方法和乘法测量误差估计方法的理论保证，确定了估计量l₁和l₂的误差界，下面将通过随机模拟实验和实际数据分析来验证我们的理论结果。

4 实验及结果分析

本文的方法简记为CoCo，Loh和Wainwright^[8]的方法记为NCL，在随机模拟实验和实际数据分析中将对两种方法进行比较。

4.1 随机模拟 4.1.1 加法测量误差模型

从可加风险模型中产生数据，设定λ₀=5，回归系数为

$ \boldsymbol{\beta}=(3, 1.5, 0, 0, 2, \cdots, 0). $

(36)

样本量n=100，样本维度p=200，X的行独立同分布，均值为0，协方差矩阵为Σ_X，考虑两种情形下的Σ_X：自回归(Σ_{X, ij}=0.5^|i-j|)和复合对称(Σ_{X, ij}=0.5+I(i=j)*0.5)，删失时间服从U(0, 2)的均匀分布使得删失率维持在20%左右。首先生成3n×p的X，然后从中选出n个满足λ₀+β^TX > 0的样本作为实验数据。加法测量误差为矩阵A，观测数据由Z=X+A生成，A的行是服从N(0, τ²I)的独立同分布变量，其中τ=0.25、0.5和0.75。

实验中，参照Datta和Zou^[9]的方法，运用5折的交叉验证方法得到惩罚参数λ及CoCo估计量的值。对于NCL的算法，参考Duchi等^[18]的工作，其中给出了详细的计算过程。首先根据式(14)得到一个NCL算法的初始估计量，即为基于$ \tilde{\boldsymbol{y}}$和$ \tilde{\boldsymbol{Z}}$的Lasso估计量。其次，依旧需要‖β_S‖₁的信息来调整参数R。由于真实系数的稀疏度信息无法提前得知，所以用一个简单的5折交叉验证从[R_max/500, 2R_max]的100均分点中选取合适的R，其中R_max代表初始估计量的l₁范数。记录C和IC分别代表选对的系数数量和错误的数量，还记录均方误差(MSE)以及其标准差(se)。总共进行了100次实验取平均数作为最后的结果。在表 1中进行展示。

表 1展示了CoCo和NCL两种方法分别在自回归和复合对称条件下的100次重复实验的结果，可以看出在两种情形下本文方法的选对数量和估计的均方误差方面都比NCL方法要好。

4.1.2 乘法测量误差模型

与加法测量误差模拟类似，依旧从可加风险模型中产生数据，λ₀=5，回归系数，样本量和样本维度都保持不变，X的行独立同分布，均值为0，协方差矩阵为Σ_X，依旧考虑Σ_X在自回归和复合对称两种条件下的情形，并且与加法测量误差中的设定保持一致。删失时间服从U(0, 2)的均匀分布使得删失率维持在20%左右，首先生成3n×p的X，然后从中选出n个满足λ₀+β^TX的作为实验数据。乘法测量误差矩阵为M=((m_{i, j}))，观测数据由Z(·)=X(·)⊙M生成，log(m_{i, j})是服从N(0, τ²I)的独立同分布变量，其中τ=0.25、0.5和0.75。与上一个随机模拟实验一样，依旧采用5折的交叉验证方法来估计CoCo估计量和NCL的参数R。同样记录C和IC分别代表选对的系数数量和错误的数量，还记录均方误差(MSE)以及其标准差(se)。总共进行100次实验取平均数作为最后的结果，在表 2中展示。

表 2展示了乘法测量误差中，CoCo和NCL两种方法分别在自回归和复合对称条件下的100次重复实验结果，可以看出在两种情形下本文方法的选对数量和估计的均方误差都比NCL方法要好。但是随着测量误差变大，CoCo和NCL方法的估计精确度都会有明显下降。

4.2 实际数据分析

在这一节中，会分析van Houwelingen Hans C.等^[19]乳腺癌和基因的关系数据。这个数据包含了从荷兰乳腺研究中心的295个女性乳腺癌患者样本，诊断时间为1984—1995年。295个患者的随访时长中位数为6.7 a。所有的肿瘤都由一个包含24 885个基因的cDNA序列来描述。van't Veer L J等^[20]在这个数据只有78个肿瘤的初始工作中，运用Rosetta误差模型根据p值大小选取了4 919个基因。本文的数据沿用这一方法，并在此基础上根据Gorst-Rasmussen和Scheike的ISIS方法^[5]选择了500个和生存时间相关程度最高的基因，因此样本维度是p=500。类似于前面的随机模拟实验，人工地加入加法测量误差来产生被污染的协变量，控制加法测量误差的参数τ=0.75/$ \sqrt{n}$，这是因为变量经过了标准化处理使得每一列的L₂范数是1。

为了检验我们方法的有效性，将295个样本随机分成包含235个样本的训练集和60个样本的验证集并重复100次，在每一次实验中，都采用随机模拟实验中的两种方法，即CoCo和NCL，用训练集训练模型参数并用验证集来筛选表现最好的估计量。计算

$ \mathrm{SEE}=\left|\frac{1}{2} \boldsymbol{\beta}^{\mathrm{T}} \boldsymbol{V} \boldsymbol{\beta}-\boldsymbol{b}^{\mathrm{T}} \boldsymbol{\beta}\right| $

(37)

作为检验两种方法效果的指标。具体的结果展示在表 3中。从表 3中可以看出我们的方法依旧有比较高的预测精确度，这也和随机模拟实验的结果相符。我们方法的指标相比NCL方法要好一些，并且变量选择的数量上也比较相近。

5 结论

本文提出一种针对高维可加风险模型中带有测量误差情况下的变量选择方法。在已知的生存分析数据相关文献中，尚未有针对测量误差数据的变量选择方法。本文基于高维线性模型测量误差数据的估计方法，重构了高维可加风险模型，并给出了加法和乘法两种测量误差模型的变量选择算法。简化伪得分方程的形式更加简洁且实用性强。随机模拟实验和实际数据分析的相关结果证实了本文方法的有效性和精确性。

在未来的工作中，我们将致力于将简化伪得分方程应用于高维可加风险模型的变量选择中。同时也会对Cox模型，加速失效模型等其他生存分析模型中的测量误差数据利用最近邻半正定投影的方法进行变量选择方面的探索。

附录A 最近邻半正定投影矩阵的计算方法

采用ADMM(alternating direction method of multipliers)算法来解决如下优化问题

$ \hat{\boldsymbol{A}}=\underset{A \geqslant \varepsilon I}{\operatorname{argmin}}\|\boldsymbol{A}-\hat{\boldsymbol{\varSigma}}\|_{\max } \text {. } $

(A.1)

引进一个额外的变量B将上述优化问题重新写作

$ (\hat{\boldsymbol{A}}, \hat{\boldsymbol{B}})=\underset{\boldsymbol A \geqslant \varepsilon \boldsymbol I, \boldsymbol B=\boldsymbol A-\hat{\boldsymbol{\varSigma}}}{\operatorname{argmin}}\|\boldsymbol{B}\|_{\max } . $

(A.2)

上述优化问题可以采用增广拉格朗日方程求解

$ f(\boldsymbol{A}, \boldsymbol{B}, \boldsymbol{\varLambda})=1 / 2|\boldsymbol{B}|_{\max }-\langle\boldsymbol{\varLambda}, \boldsymbol{A}-\boldsymbol{B}-\hat{\boldsymbol{\varSigma}}\rangle+\frac{1}{2 \mu}\|\boldsymbol{A}-\boldsymbol{B}-\hat{\boldsymbol{\varSigma}}\|_F^2 . $

(A.3)

其中: μ是惩罚参数，Λ是拉格朗日矩阵，〈·〉代表内积，‖·‖_F代表Frobenius范数。用下面的算法来求解该优化问题。

算法1：最近邻半正定投影矩阵的ADMM算法
1. 输入μ和初始值B0以及Λ0。
2. 在第i步更新：
2.1 (步骤A)$ \boldsymbol{A}_{i+1}=\left(\boldsymbol{B}_i+\hat{\boldsymbol{\varSigma}}+\mu \boldsymbol{\varLambda}_i\right)_{\varepsilon}$,
2.2 (步骤B)$\left.\boldsymbol{B}_{i+1}=\operatorname{mat}_l\left(\operatorname{vec}_l\left(\boldsymbol{A}_{i+1}-\hat{\boldsymbol{\varSigma}}-\mu \boldsymbol{A}_i\right)\right)-l_1\left(\operatorname{vec}_l\left(\boldsymbol{A}_{i+1}-\hat{\boldsymbol{\varSigma}}-\mu \boldsymbol{\varLambda}_i, \mu\right)\right)\right) $,
2.3 (步骤Λ)$ \boldsymbol{\varLambda}_{i+1}=\boldsymbol{\varLambda}_i-\frac{\boldsymbol{A}_{i+1}-\boldsymbol{B}_{i+1}-\hat{\boldsymbol{\varSigma}}}{\mu}$
3. 重复第2步直到收敛。

其中vec_l(M)代表将对称阵M的下半部分向量化。mat_l(M)代表vec_l(M)的逆过程。如果M=∑_jλ_jp_jp_j^'代表了矩阵的谱分解，则M_ε=∑_jmax(λ_j, ε)p_jp_j^'。

附录B 相关定理的证明

本文在固定设计的前提下来进行理论部分的证明，即$ \frac{1}{n}\left\|x_j\right\|_2^2=1$对1≤j≤p成立。且假设观察时间的最大值t_n < ∞，加法测量误差|a_ij|≤L。

定理3.1的证明   定理3.1为Datta和Zou^[9]中定理1的平行定理，证明过程与Datta和Zou^[9]中定理1的证明过程类似，故此处省略。我们将详细证明两种测量误差的估计量满足引理3.1和引理3.2。

引理3.1的证明   令Σ_A=((σ_{a, ij}))，则

$ \hat{V}_{\mathrm{add}, j k}-V_{j k}=\sum\nolimits_{i=1}^n \frac{t_i}{n} a_{i j} x_{i j}+\sum\nolimits_{i=1}^n \frac{t_i}{n} a_{i k} x_{i k}+\sum\nolimits_{i=1}^n \frac{t_i}{n}\left(a_{i k} a_{i j}-\sigma_{a, j k}\right), $

(B.1)

由于$ \frac{1}{n}\left\|\sum_{i=1}^n t_i x_{i j}\right\|_2^2 \leqslant t_n$，则由引理B.2可知$ |\sum_{i=1}^n \frac{t_i}{n} a_{i j} x_{i j}|, | \sum_{i=1}^n \frac{t_i}{n} a_{i k} x_{i k} \mid$至多以概率Cexp(-cnε²/t_n²τ²)大于ε/3。同理，由引理B.1，$ \sum_{i=1}^n \frac{t_i}{n}\left(a_{i k} a_{i j}-\sigma_{a, j k}\right)$也可以被证明是一个小量，综上所述，$\hat{\boldsymbol{V}}_{\text {add }} $满足近邻条件式(30)，且

$ \varsigma=\max \left(t_n^2 \tau^4, t_n^2 \tau^2\right), \varepsilon_0=c \tau^2 t_n. $

(B.2)

对于$ \hat{b}_{\mathrm{add}, j}-b_j=\sum_{i=1}^n \frac{\delta_i t_t}{n} a_{i j}$, 其中a_ij为相互独立且均值为0的随机变量序列，则由Hoeffding不等式有$\left|\hat{b}_{\text {add }, j}-b_j\right| $至多以概率exp(-cε²n³/L⁴)大于ε。综上所述，$ \hat{\boldsymbol{b}}_{\text {add }}$满足近邻条件式(31)，且

$ \varsigma=\max \left(t_n^2 \tau^4, t_n^2 \tau^2, L^4\right), \varepsilon_0=c \tau^2 t_n. $

综上，引理3.1的结论已论证完毕。

引理3.2的证明   令Σ_M=((σ_{m, ij}))，则有

$ \begin{aligned} \hat{V}_{\text {mult }, j k}-V_{j k}= & \frac{1}{n} \sum\nolimits_{i=1}^n \frac{t_i x_{i j} x_{i k}}{\mu_j \mu_k+\sigma_{m, j k}}\left(m_{i j} m_{i k}-\mu_j \mu_k-\sigma_{m, j k}\right) \\ = & \frac{1}{n} \sum\nolimits_{i=1}^n \frac{t_i x_{i j} x_{i k}}{\mu_j \mu_k+\sigma_{m, j k}}\left(\left(m_{i j}-\mu_j\right)\left(m_{i k}-\mu_k\right)-\sigma_{m, j k}\right)+ \\ & \frac{1}{n} \sum\nolimits_{i=1}^n \frac{t_i x_{i j} x_{i k}}{\mu_j \mu_k+\sigma_{m, j k}}\left(\mu_j\left(m_{i k}-\mu_k\right)+\mu_k\left(m_{i j}-\mu_j\right)\right) . \end{aligned} $

(B.3)

利用式(9)的正则化条件，有

$ \begin{gathered} \left|\hat{V}_{\operatorname{mult}, j k}-V_{j k}\right| \leqslant\left|\frac{1}{n M_{\min }} \sum\nolimits_{i=1}^n t_i x_{i j} x_{i k}\left(\left(m_{i j}-\mu_j\right)\left(m_{i k}-\mu_k\right)-\sigma_{m, j k}\right)\right|+ \\ \frac{\mu_{\max }}{n M_{\min }} \sum\nolimits_{i=1}^n t_i x_{i j} x_{i k}\left(m_{i k}-\mu_{i k}\right)+\frac{\mu_{\max }}{n M_{\min }} \sum\nolimits_{i=1}^n t_i x_{i j} x_{i k}\left(m_{i k}-\mu_{i k}\right) . \end{gathered} $

(B.4)

将上式右边3项分别记为T₁, T₂和T₃。对于T₁，首先有‖v‖_∞≤X_max²，其中v_i=x_ijx_ik，则和引理3.1证明过程类似，继续由引理B.1可以得到下面结论，对于$\varsigma=\max \left(t_n{ }^2 \tau^2 X_{\max }^2 \frac{\mu_{\max }^2}{M_{\min }^2}, t_n{ }^2 \tau^4 \frac{X_{\max }^4}{M_{\min }^2}\right) $和$ \varepsilon_0=c \tau^2 t_n X_{\max }^2$，有

$ P\left(T_1 \geqslant \varepsilon\right) \leqslant C \exp \left(-c n \varepsilon^2 \varsigma^{-1}\right) . $

(B.5)

对于T₂和T₃仿照引理3.1的证明过程即可得到类似结论。由于$ \hat{\boldsymbol{b}}_{\text {mult }}$相比于$ \hat{\boldsymbol{b}}$只是对应元素除去了μ_M，并且在条件(9)中μ_M的各个元素上下界已经做了限制，所以对于|$ \hat{\boldsymbol{b}}_{\text {mult }}$-$ \hat{\boldsymbol{b}}$|的界只需简单应用引理B.2即可直接得到结论。综上所述，引理3.2证明完毕。

引理B.1   令z_i=(x_i, y_i)′代表独立同分布的，均值为0，协方差矩阵为Σ=((σ_ij))，次高斯参数为τ²的向量。则存在绝对常数C和c使得，对于任意的ε≤cτ²‖a‖_∞，有

$ P\left(\frac{1}{n}\left|\sum\nolimits_{i=1}^n a_i\left(x_i y_i-\sigma_{12}\right)\right| \geqslant \varepsilon\right) \leqslant \exp \left(\frac{-n c \varepsilon^2}{\tau^4\|a\|_{\infty}^2}\right) . $

(B.6)

定义B.1(次高斯随机变量)   Z是一个随机变量，若存在有限的κ > 0使得$ \kappa=\sup\limits _{p \geqslant 1} p^{-1 / 2}\left(E|Z|^p\right)^{1 / p}$，则Z被称为次高斯随机变量，κ被称为Z的次高斯范数，记为‖Z‖_φ。

引理B.2   一个次高斯随机变量Z满足如下的尾概率界

$ P(|Z|>t) \leqslant 2 \exp \left(-t^2 / 2 \tau^2\right), \forall t>0. $

(B.7)

其中满足上式的最小的τ被称为次高斯参数，如果ω=(ω₁, …ω_n)′并且所有的ω_i是独立的0均值的随机次高斯变量，则有一个重要的性质

$ \left\|\boldsymbol{v}^{\prime} \boldsymbol{\omega}\right\|_{\varphi}^2 \leqslant K\|\boldsymbol{v}\|_2^2 \max \limits_i\left(\left\|\omega_i\right\|_{\varphi}^2\right) . $

(B.8)

引理B.1和引理B.2均引自Datta和Zou^[9]工作中的证明部分，此处省略相关证明。