基于Cox回归模型的约束估计问题的研究

The Constraint Estimation Based on Cox Regression Model

摘要: 具有约束条件的右删失Cox回归模型的研究经常是基于偏似然函数，而该似然函数忽略了未知的基准风险率函数。文章致力于解决两个问题：基于完全似然函数得到回归参数和累积风险函数的极大似然估计；把有限制条件的优化问题转化为无限制条件的优化问题。ADMM算法能把限制条件引入到目标函数中，从而把有条件的优化问题转化为无条件的优化问题，MM算法在解决优化问题方面具有可分离参数的优点。因此，文章首先利用ADMM算法把有限制条件的优化问题转化为无限制条件的优化问题，然后将MM算法应用于极大化新的目标函数，实现了参数和非参数的分离，进而解决了回归参数和累积风险函数的非参数估计问题。同时，利用不等式放缩把高维优化问题转化为一维优化问题，避开了矩阵求逆的困难。

Abstract: The research on the right-censored Cox regression model with constraints is often based on the partial likelihood function, whereas the likelihood function ignores the unknown baseline hazard rate function. We devote to solving two problems. One is to obtain the maximum likelihood estimations of regression parameters and cumulative hazard function based on the complete likelihood function. Another is to transform the optimization problem with restrictions into the unrestricted optimization problem. The ADMM algorithm can introduce constraints into the objective function so as to transform the conditional optimization problem into an unconditional optimization problem. The MM algorithm has the advantage of separating parameters in solving optimization problems. Therefore, we first use the ADMM algorithm to transform the optimization problem with restrictions into an unconditional optimization problem. Then we apply the MM algorithm to maximize the new objective function in order to separate the parameters and the non-parameters, which is helpful to solve the problems of estimating the regression parameters and nonparametric cumulative hazard function. Meanwhile, the use of inequality which can transform the high-dimensional function into a sum of one-dimensional functions avoids the difficulty of matrix inversion.

1. 引言

Cox回归模型自1972年被Cox [1] 首次提出后，由于该模型具有优良的理论性质，以及直观和解释性强的优点，已被广泛应用于生物医学、可靠性分析和经济领域等。Cox模型是一种半参数生存分析模型，是分析生存数据最为广泛的统计模型之一。但是，在实际的生产实践中，由于实验期限已到、实验对象被移除或者提前退出实验等种种原因，导致观察收集到的生存数据往往又以右删失为主。因此，涌现了大量的科研工作者基于Cox回归模型对右删失数据进行分析研究。比如，王佳等 [2] 在其文章中基于Cox回归模型分析了三组医学类数据：乳腺癌数据、骨髓移植数据和原发性胆汁肝胆硬化数据。洪东跑等 [3] 证实了Cox回归模型可用于分析产品可靠性关于环境因素的动态变化特征。马超群和何文 [4] 基于Cox回归模型就上市公司财务困境问题进行判别能力和稳定性分析，等等。

在实际问题中，若我们提前知道有关Cox模型参数的一些信息，加以利用定会得到更加可靠的统计推断结果。然而，除了Ding等 [5]，目前有关这方面的研究却比较少。但是，该文章是基于Cox模型的偏似然函数进行分析，而忽略了基准风险率函数的非参数估计问题，这显然具有一定的不合理性。因此，本文将基于Cox模型的完全似然函数，对具有限制条件的模型参数进行极大化估计。

显然，本文的难点之一在于极大化完全似然函数时，要考虑到限制条件对估计值的影响。幸而，ADMM算法给出了一个较好的解决办法，即把有条件的优化问题转化为无条件的优化问题，进而避免了在求解过程中要考虑约束条件的限制。事实上，ADMM算法已被广泛应用于解决各种带有约束条件的优化问题。例如，图像提取时的稀疏lasso和凸聚类等问题 [6]，高维稀疏惩罚分位数回归模型 [7]，基于凸函数的子组分析模型 [8]，等等。这些文章都是运用ADMM算法把有限制条件的优化问题转化为无限制条件的优化问题，该算法提高了模型求解的速度。而且，ADMM算法已在理论上被证明了其可靠性 [9]。

再者，转化后的无条件优化问题没有解析解，故需要应用Newton-Raphson算法(以下简称“NR算法”)进行迭代求解。但是，在面向高维数据集时，NR算法需要在矩阵可逆的条件下才具有可行性。由于minorization-maximization算法(以下简称“MM算法”)具有可分离参数的优点 [10]，通过将目标函数进行参数分离，从而绕开了矩阵求逆问题。因此，本文将用两种形式的MM算法来解决基于Cox模型的完全似然函数的无条件优化问题，其主要区别在于处理基准风险率函数的非参数估计问题上的不同 [11]。一种方法运用profile估计方法(Johansen [12] )来估计基准风险率函数，另一种方法则是绕开了profile估计方法，直接应用MM算法把累积风险率从回归参数中分离出来。

本文剩余章节安排如下：第2节介绍受约束的Cox回归模型及其完全似然函数；第3节将ADMM算法和MM算法相结合，从而得到回归参数和非参数的累积风险函数的极大似然估计；第4节进行数值模拟，用于评估两种算法相结合的估计效果；第5节对一组癌症数据进行分析研究；第6节给出本文结论。

2. 受约束Cox回归模型

假设从某个总体中随机地抽取n个样本，并用 $T_i$， $C_i$ 和 $X_i={\left(x_{i1},\cdots ,x_{iq}\right)}^{\prime }$ 分别表示第i个样本的生存时间，删失时间以及q维解释变量值。进一步，假定在 $X_i$ 给定的条件下，删失时间 $C_i$ 独立于生存时间 $T_i$。从而，数据集 $Y_{obs}=\left\{\left(Y_i=\min \left(T_i,C_i\right),I_i,X_i\right),i=1,\cdots ,n\right\}$，其中 $Y_i$， $I_i=I\left(T_i\le C_i\right)$ 和 $I(\cdot )$ 分别表示观测时间，右删失变量以及示性函数。因此，给定 $X_i$ 的风险率函数为：

$\lambda \left(t|X_i\right)={\lambda }_0\left(t\right)\exp \left({X^{\prime }}_i\beta \right)$

其中 ${\lambda }_0(\cdot )$ 是基准风险率函数， $\beta ={\left({\beta }_1,\cdots ,{\beta }_q\right)}^{\prime }$ 是回归参数。

根据Kim等 [13]，可得其完全似然函数为

$L\left(\beta ,\Lambda |Y_{obs}\right)={\displaystyle \underset{i=1}{\overset{n}{\prod }}{\left({\lambda }_0\left(t_i\right)\exp \left({X^{\prime }}_i\beta \right)\right)}^{I_i}\exp \left(-\Lambda \left(t_i\right)\exp \left({X^{\prime }}_i\beta \right)\right)}$

其中 $\Lambda (\cdot )$ 为累积风险函数。因此，其对数似然函数为

$l\left(\beta ,\Lambda |Y_{obs}\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(I_i\ln \left({\lambda }_0\left(t_i\right)\right)+I_i{X^{\prime }}_i\beta -\Lambda \left(t_i\right)\exp \left({X^{\prime }}_i\beta \right)\right)}$ (1)

我们感兴趣的是对受到某些条件限制的回归参数 $\beta$ 的估计，比如边界限制条件或线性不等式的条件限制。在本文中，我们只考虑 $\beta \ge 0$ 的约束条件，故本文最终考虑如下的基于完全似然函数的受约束Cox模型，即：

$\underset{\beta \ge 0}{\text{max}}l\left(\beta ,\Lambda |Y_{obs}\right)$ (2)

3. ADMM + MM算法

根据Boyd等 [9] 在其文章的第5章中提出的理论依据，求解上述的受约束Cox模型等价于求解如下的模型：

$\begin{array}{l}\max l\left(\beta ,\Lambda |Y_{obs}\right)+g\left(z\right)\\ \beta -z=0\end{array}$ (3)

其中 $g\left(z\right)=I\left(z\ge 0\right)$。进而转化为极大化如下的无约束的目标函数：

$l_{\rho }\left(\alpha |Y_{obs}\right)=l\left(\beta ,\Lambda |Y_{obs}\right)+g\left(z\right)-\frac{\rho }{2}{\Vert \beta -z+u\Vert }_2^2$ (4)

其中 $\alpha =\left(\beta ,\Lambda ,z,u\right)$， ${\Vert \cdot \Vert }_2$ 为向量的二范数。再根据Boyd等 [12] 的第5章的内容，可得极大化(4)式的ADMM形式为：

$\left({\beta }^{\left(k\text{+}1\right)},{\Lambda }^{\left(k\text{+}1\right)}\right)=\underset{\left(\beta ,\Lambda \right)}{\arg \max }\left(l\left(\beta ,\Lambda |Y_{obs}\right)-\frac{\rho }{2}{\Vert \beta -z^{\left(k\right)}+u^{\left(k\right)}\Vert }_2^2\right)$, (5)

$z^{\left(k+1\right)}={\left({\beta }^{\left(k+1\right)}+u^{\left(k\right)}\right)}_{+}$, (6)

$u^{\left(k+1\right)}=u^{\left(k\right)}+{\beta }^{\left(k+1\right)}-z^{\left(k+1\right)}$. (7)

观察上述的迭代过程可知，难点在于求解 $\left(\beta ,\Lambda \right)$ 的极大似然估计。若直接利用NR算法，在高维情况下会涉及矩阵求逆问题，且非参数 $\Lambda$ 的估计也是个棘手问题。幸而，近年发展迅速的MM算法为上述问题提供了解决办法。该算法的实现主要依赖于寻找到合适的不等式，进而通过一系列不等式放缩构造目标优化函数 $l_{\rho }\left(\alpha |Y_{obs}\right)$ 的最小化函数(也叫替代函数) $Q\left(\alpha |{\alpha }^{\left(k\right)}\right)$，使其满足

$\{\begin{array}{l}{l}_{\rho }\left(\alpha |Y_{obs}\right)\ge Q\left(\alpha |{\alpha }^{\left(k\right)}\right),\forall \alpha ,{\alpha }^{\left(k\right)}\in \Theta \\ l_{\rho }\left(\alpha |Y_{obs}\right)=Q\left({\alpha }^{\left(k\right)}|{\alpha }^{\left(k\right)}\right)\end{array}$ (8)

其中 ${\alpha }^{\left(k\right)}$ 为第k次对 $\hat{\alpha }$ 的近似值。然后，我们最大化替代函数，可得到如下迭代公式

${\alpha }^{\left(k+1\right)}=\arg \underset{\alpha \in \Theta }{\max }Q\left(\alpha |{\alpha }^{\left(k\right)}\right)$ (9)

当第k + 1次逼近 $\hat{\alpha }$ 时，有

$l_{\rho }\left({\alpha }^{\left(k+1\right)}|Y_{obs}\right)\ge Q\left({\alpha }^{\left(k+1\right)}|{\alpha }^{\left(k\right)}\right)\ge Q\left({\alpha }^{\left(k\right)}|{\alpha }^{\left(k\right)}\right)=l_{\rho }\left({\alpha }^{\left(k\right)}|Y_{obs}\right)$ (10)

其中迭代过程(9)呈现上升趋势，使得目标函数不断增大，从而实现优化转移，即极大化 $Q\left(\cdot |{\alpha }^{\left(k\right)}\right)$ 等价于极大化 $l_{\rho }\left(\cdot |Y_{obs}\right)$。

使用MM算法时，关键且困难的一步是寻找合适的替代函数。在实际过程中，尤其遇到高维的情况，构建一个参数可分离的替代函数尤为重要，因为这样可以避免矩阵求逆的困难。而本文从Huang等 [11] 处得到启发，利用MM算法把参数 $\beta$ 和非参数 $\Lambda$ 进行参数分离，且进一步把 $\beta$ 的各个分量进行分离，从而避免了矩阵求逆问题。

3.1. ADMM + profile MM算法

如上所述，MM算法的关键是构造目标函数 $l\left(\beta ,\Lambda |Y_{obs}\right)-\frac{\rho }{2}{\Vert \beta -z^{\left(k\right)}+u^{\left(k\right)}\Vert }_2^2$ 的替代函数。

首先，观察 $l^A\left(\beta ,\Lambda |{\alpha }^{\left(k\right)}\right)=l\left(\beta ,\Lambda |Y_{obs}\right)-\frac{\rho }{2}{\Vert \beta -z^{\left(k\right)}+u^{\left(k\right)}\Vert }_2^2$ 可得，Johansen [12] 在其文章中用到的profile方法可直接用来计算非参数 $\Lambda$，即令 $\partial l^A\left(\beta ,\Lambda |{\alpha }^{\left(k\right)}\right)/\partial {\lambda }_0\left(t_i\right)=0$ 可得

${\lambda }_0\left(t_i\right)=\frac{I_i}{{\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right)\exp \left({X^{\prime }}_r\beta \right)}},i=1,\cdots ,n$ (11)

然后把(11)式带入 $l^A\left(\beta ,\Lambda |{\alpha }^{\left(k\right)}\right)$，并把与参数 $\beta$ 无关的部分忽略不计，进而整理得到函数

$Q_1\left(\beta |{\alpha }^{\left(k\right)}\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(I_i{X^{\prime }}_i\beta -I_i\ln {\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right)\exp \left({X^{\prime }}_r\beta \right)}\right)}-\frac{\rho }{2}{\Vert \beta -z^{\left(k\right)}+u^{\left(k\right)}\Vert }_2^2$

接下来，把 $Q_1\left(\beta |{\alpha }^{\left(k\right)}\right)$ 看作新的目标函数，对 $\beta$ 的各个分量进行参数分离。根据支撑超平面不等式

$-\ln x\ge -\ln x_0-\frac{x-x_0}{x_0}$

可得

$\begin{array}{l}-\ln {\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right)\exp \left({X^{\prime }}_r\beta \right)}\\ \ge -\ln {\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right)\exp \left({X^{\prime }}_r{\beta }^{\left(k\right)}\right)}-\frac{{\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right)\exp \left({X^{\prime }}_r\beta \right)}-{\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right)\exp \left({X^{\prime }}_r{\beta }^{\left(k\right)}\right)}}{{\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right)\exp ({X^{\prime }}_r{\beta }^{(\; k\; )}}}\end{array}$

从而可得 $Q_1\left(\beta |{\alpha }^{\left(k\right)}\right)$ 的替代函数

$Q_{12}\left(\beta |{\alpha }^{\left(k\right)}\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{p=1}{\overset{q}{\sum }}\left(I_i{x}_{ip}{\beta }_p-\frac{I_i{\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right)\exp \left({X^{\prime }}_r\beta \right)}}{{\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right)\exp \left({X^{\prime }}_r{\beta }^{\left(k\right)}\right)}}\right)}}-\frac{\rho }{2}{\Vert \beta -z^{\left(k\right)}+u^{\left(k\right)}\Vert }_2^2$

然后，基于 $X_i^{\text{T}}\beta$ 构造权重 ${\omega }_{ip}=\left|x_{ip}\right|/{\displaystyle \underset{p=1}{\overset{q}{\sum }}\left|x_{ip}\right|}$，再把 $X_i^{\text{T}}\beta$ 改写为

$X_i^{\text{T}}\beta ={\displaystyle \underset{p=1}{\overset{q}{\sum }}{\omega }_{ip}\left[{\omega }_{ip}^{-1}{x}_{ip}\left({\beta }_p-{\beta }_p^{\left(k\right)}\right)+X_i^{\text{T}}{\beta }^{\left(k\right)}\right]}$

如果 $x_{ip}=0$，则令 $1/{\omega }_{ip}=0$。然后根据离散型Jensen不等式

$\phi \left({\displaystyle \underset{i=1}{\overset{n}{\sum }}{a}_i{x}_i}\right)\ge {\displaystyle \underset{i=1}{\overset{n}{\sum }}{a}_i\phi \left(x_i\right)}$，其中 $a_i\ge 0$，且 ${\displaystyle \underset{i=1}{\overset{n}{\sum }}{a}_i}=1$

得到 $Q_{12}\left(\beta |{\alpha }^{\left(k\right)}\right)$ 的替代函数

$Q_{13}\left(\beta |{\alpha }^{\left(k\right)}\right)\stackrel{\wedge}{=}{\displaystyle \underset{p=1}{\overset{q}{\sum }}{Q}_{13p}\left({\beta }_p|{\alpha }^{\left(k\right)}\right)}$ (12)

其中

$Q_{13p}\left({\beta }_p|{\alpha }^{\left(k\right)}\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(I_i{x}_{ip}{\beta }_p-\frac{I_i{\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right){\omega }_{rp}\exp \left({\omega }_{rp}^{-1}{x}_{rp}\left({\beta }_p-{\beta }_p^{\left(k\right)}\right)+{X^{\prime }}_r{\beta }^{\left(k\right)}\right)}}{{\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right)\exp \left({X^{\prime }}_r{\beta }^{\left(k\right)}\right)}}\right)}-\frac{\rho }{2}{\left({\beta }_p-z_p^{\left(k\right)}+u_p^{\left(k\right)}\right)}^2$ (13)

从(12)和(13)式可以看出，目标函数 $l^A\left(\beta ,\Lambda |{\alpha }^{\left(k\right)}\right)$ 已经被分解成q个单变量函数之和，实现了把高维优化问题转化为低维优化问题。所得的MM算法在其最大化步骤中仅涉及q个单独的单变量优化问题，因此不需要矩阵求逆。令 $\text{d}{Q}_{13p}\left({\beta }_p|{\alpha }^{\left(k\right)}\right)/\text{d}{\beta }_p=0$，再利用NR算法进行迭代求解。

3.2. ADMM + non-profile MM算法

本节提出一种绕过profile方法的non-profile MM算法。

首先，根据算术几何不等式

$-{\displaystyle \underset{i=1}{\overset{n}{\prod }}{x}_i^{a_i}}\ge -{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{a_i}{{\Vert a\Vert }_1}{x}_i^{{\Vert a\Vert }_1}}$

其中 ${\Vert \cdot \Vert }_1$ 为向量的一范数。可得 $l^A\left(\beta ,\Lambda |{\alpha }^{\left(k\right)}\right)$ 的替代函数

$Q_2\left(\alpha |{\alpha }^{\left(k\right)}\right)\stackrel{\wedge}{=}{Q}_{21}\left(\Lambda |{\alpha }^{\left(k\right)}\right)+Q_{22}\left(\beta |{\alpha }^{(\; k\; )}\right)$

其中

$\begin{array}{l}{Q}_{21}\left(\Lambda |{\alpha }^{\left(k\right)}\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(I_i\ln \left({\lambda }_0\left(t_i\right)\right)-\frac{\exp \left({X^{\prime }}_i{\beta }^{\left(k\right)}\right)}{2{\Lambda }^{\left(k\right)}\left(t_i\right)}{\Lambda }^2\left(t_i\right)\right)}\\ Q_{22}\left(\beta |{\alpha }^{\left(k\right)}\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(I_i{X}_i{}^{\prime }\beta -\frac{{\Lambda }^{\left(k\right)}\left(t_i\right)}{2\exp \left({X^{\prime }}_i{\beta }^{\left(k\right)}\right)}\exp \left(2{X^{\prime }}_i\beta \right)\right)}-\frac{\rho }{2}{\Vert \beta -z^{\left(k\right)}+u^{\left(k\right)}\Vert }_2^2\end{array}$

为了得到非参数 $\Lambda$ 的估计，可极大化 $Q_{21}\left(\Lambda |{\alpha }^{\left(k\right)}\right)$，即令 $\text{d}{Q}_{21}\left(\Lambda |{\alpha }^{\left(k\right)}\right)/\text{d}{\lambda }_0\left(t_i\right)=0$。并进一步整理可得 ${\lambda }_0\left(t_i\right)$ 的估计

${\lambda }_0\left(t_i\right)=\frac{I_i}{{\displaystyle \underset{r=1}{\overset{n}{\sum }}I\left(t_r\ge t_i\right)\exp \left({X^{\prime }}_r{\beta }^{\left(k\right)}\right)}}$ (14)

类似于上一节的方法技巧，对 $\beta$ 的各个分量进行参数分离。首先把 $2X_i^{\text{T}}\beta$ 改写为

$2X_i^{\text{T}}\beta ={\displaystyle \underset{p=1}{\overset{q}{\sum }}{\omega }_{ip}\left[2{\omega }_{ip}^{-1}{x}_{ip}\left({\beta }_p-{\beta }_p^{\left(k\right)}\right)+2X_i^{\text{T}}{\beta }^{\left(k\right)}\right]}$，然后把离散型Jensen不等式运用到函数 $-\exp (\cdot )$，从而得 $Q_{22}\left(\beta |{\alpha }^{\left(k\right)}\right)$ 的替代函数为

$Q_{23}\left(\beta |{\alpha }^{\left(k\right)}\right)\stackrel{\wedge}{=}{\displaystyle \underset{p=1}{\overset{q}{\sum }}{Q}_{23p}\left({\beta }_p|{\alpha }^{\left(k\right)}\right)}$ (15)

其中

$Q_{23p}\left({\beta }_p|{\alpha }^{\left(k\right)}\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(I_i{X}_{ip}{\beta }_p-\frac{{\Lambda }^{\left(k\right)}\left(t_i\right){\omega }_{ip}}{2\exp \left({X^{\prime }}_i{\beta }^{\left(k\right)}\right)}\exp \left(2{\omega }_{ip}^{-1}{x}_{ip}\left({\beta }_p-{\beta }_p^{\left(k\right)}\right)+2{X^{\prime }}_i{\beta }^{\left(k\right)}\right)\right)}-\frac{\rho }{2}{\left({\beta }_p-z_p^{\left(k\right)}+u_p^{\left(k\right)}\right)}^2$ (16)

观察(15)和(16)式可以看出，我们得到了与(12)和(13)式类似的结果，即目标函数已经实现参数分离，在接下来的最大化步骤中绕开了矩阵求逆。同理，令 $\text{d}{Q}_{23p}\left({\beta }_p|{\alpha }^{\left(k\right)}\right)/\text{d}{\beta }_p=0$，再利用NR算法进行迭代求解即可。

综上所述，我们得到如下的迭代算法：

第一步：给定初始值 $\left({\beta }^{\left(0\right)},{\Lambda }^{\left(0\right)},z^{\left(0\right)},u^{\left(0\right)}\right)$ ；

第二步：基于式(11)或式(14)计算 ${\Lambda }^{\left(k+1\right)},k=0,1,2,\cdots$ ；

第三步：对应于第二步，利用NR算法对式(13)或式(16)进行迭代求解 ${\beta }_p^{\left(k+1\right)}$，其中 $p=1,\cdots ,q;k=0,1,2,\cdots$ ；

第四步：基于式(6)计算 $z^{\left(k+1\right)},k=0,1,2,\cdots$ ；

第五步：基于式(7)计算 $u^{\left(k+1\right)},k=0,1,2,\cdots$ ；

第六步：重复第二到五步，直到算法收敛。

4. 数值模拟

模拟研究设定为：参数向量 $\beta ={\left(1,2\right)}^{\prime }$，则 $q=2$，且q个解释变量独立同分布，其数值由正态分布 $N\left(-1,1\right)$ 随机产生。基准风险率函数设置为 ${\lambda }_0\left(t\right)=2$，则对应的累积风险函数 $\Lambda \left(t\right)=2t$。在整个模拟实验中，我们通过控制样本量n和删失率r来评估Cox回归模型分别在无约束和有约束的情况下，对回归参数和累积风险率 $\left(\Lambda \left(0.1\right),\Lambda \left(0.5\right)\right)=\left(0.2,1\right)$ 的估计效果。其中，样本量n分别取50和100；删失率r分别达到20%，50%以及70%。

整个模拟由R软件实现，所有实验的收敛条件均为

$\frac{\left|l\left({\beta }^{\left(k+1\right)}|Y_{obs}\right)-l\left({\beta }^{\left(k\right)}|Y_{obs}\right)\right|}{\left|l\left({\beta }^{\left(k\right)}|Y_{obs}\right)\right|+1}<{10}^{-8}$

每种情况都进行 $R=500$ 次的模拟研究。具体模拟结果见表1和表2。

Table 1. The simulations of the unconstrained and constrained Cox model under the profile MM algorithm

表1. profile MM算法下的无约束条件和有约束条件的Cox模型的模拟结果

表1是基于profile MM算法进行500次模拟研究得到的结果，其中有约束条件的Cox模型的极大化估计则结合了ADMM算法，具体过程见3.1。表格中的数据来源于对500次模拟得到的估计值与真值的偏差BIAS和均方误差MSE求平均，以及估计值的样本差S.D.。从表1可以看出，无论样本量和删失率达到多少，对比BIAS数值可以看出，考虑受约束条件 $\beta \ge 0$ 的模拟结果比无约束条件的模拟结果总体上更接近真值。而在约束条件下得到的S.D.和MSE数均比无约束条件下的小，故提前知道回归参数的部分信息有助于得到更可靠的估计值。其次，MM算法和ADMM + MM算法在估计回归参数 $\beta$ 和非参数 $\Lambda$ 方面都表现出良好的估计效果。再者，通过对比观察可知，固定删失率，增加样本量，两种情况下模拟得到的估计值与真值之间的偏差越来越小，RMSE和S.D.的值也越来越小，即估计效果和模拟的稳定性越来越好。固定样本量，当删失率越小，评价指标BIAS、RMSE和S.D.的值也越来越小，这显然符合实际情况。

表2是基于non-profile MM算法的模拟结果，其中有约束条件的Cox模型的极大化估计则结合了ADMM算法，具体过程见3.2。其数据来源于对500次模拟得到的估计值与真值的偏差BIAS和MSE求平均，以及估计值的样本差S.D.。对比表2中的BIAS数值，有约束条件下的估计值比无约束条件下的估计值更接近真值，且S.D.和MSE数值也相对更小。这表示极大化有约束条件的Cox模型的完全似然函数更能得到可靠的结果。类似地可以看出，MM算法和ADMM + MM算法在估计回归参数 $\beta$ 和非参数 $\Lambda$ 方面都表现出良好的估计效果。以及减小删失率或增加样本量，都使得估计值与真值之间的差异性越来越小，模拟效果越来越好。

Table 2. The simulations of the unconstrained and constrained Cox model under the non-profile MM algorithm

表2. non-profile MM算法下的无约束条件和有约束条件的Cox模型的模拟结果

综上所述，无论样本量n取50还是100，删失率r为20%，50%还是70%，极大化有约束条件的Cox回归模型的完全似然函数时，都表现出更好的估计效果。特别地，观察本实验中最不理想的组合 $\left(n,r\right)=\left(50,70\%\right)$，其模拟结果依然是令人满意的。由此可见，ADMM算法和MM算法相结合的效果良好。进一步，该模拟结果告诉我们，可以通过增加样本量、减少删失率来改善估计效果，使其更加接近实际。

5. 实例分析

安装R语言中的survival程序包，调用内置的cancer数据集。该数据集共记录了228个来自中北部癌症治疗组的癌症晚期患者的相关数据 [14]，其中包括机构代码、生存天数、患者生存状态(死亡或者实验数据删失)、年龄、性别、由医生评定的GCOG表现评分、医生评定的Karnofsky表现评分、由患者自己评定的Karnofsky表现评分、进食时消耗的卡路里以及最近6个月的体重下降数。

在本例中，我们只研究患者最近6个月的体重下降数与死亡风险率之间的关系，即回归变量X只取一维。从而，建立如下的Cox回归模型：

$\lambda \left(t|X_i\right)={\lambda }_0\left(t\right)\exp \left(X_i\beta \right)$

然后根据对以上模型进行参数估计的结果可知(见表3中无约束条件下的估计结果)，本文可对回归参数做如下限制：

$\beta \ge 0$。

从而，针对有限制条件的Cox模型的估计，可用本文提出的ADMM + MM (profile和non-profile)算法进行回归参数估计和累积风险函数的估计。由于该模型具有大样本性质，故可用bootstrap方法 [5] 求得样本估计值的标准差和95%置信区间。最终得到如表3所示的数值分析结果。从表3中可以看出，两种ADMM + MM算法的估计结果在精确到小数点三位时是一致的，该结果表明最近6个月患者的体重下降越多，越会增加患者的死亡率，且增加到 ${\text{e}}^{0.015}\approx 1.02$ 倍，但在有约束条件下的95%置信区间的区间长度更小，这表明该估计结果更可靠。因此，我们只画出有约束条件的模型的估计累积风险函数，如图1所示，两种算法下的估计累积风险函数的数值相差不大，走势一致，这表明两种ADMM + MM算法表现相当。

Table 3. The estimations by two proposed ADMM + MM algorithms and the confidence intervals based on bootstrap method

表3. 两种ADMM + MM算法的估计结果以及bootstrap方法的置信区间

Figure 1. The estimated cumulative hazard function for patients

图1. 患者的累积死亡风险估计函数

6. 结论

正如引言所言，若提前知道部分有关未知参数的信息，加以利用必会得到更可靠的估计结果。而针对有约束条件的Cox回归模型的研究很少，而Ding等的研究却是基于该模型的偏似然函数，忽略非参数部分的估计，则无法预知个体即将面临的死亡等风险率，这显然是不合理的。再者，基于该模型的完全似然函数的有约束条件的研究却不曾有过，因此，本文基于右删失型的Cox比例风险模型，在约束条件为 $\beta \ge 0$ 的情况下极大化该模型的完全似然函数，从而得到回归参数估计和非参数估计。

针对有约束条件的Cox模型，由于ADMM算法在有约束条件优化问题方面具有优良的理论性质而得到快速发展，本文应用ADMM算法把基于右删失Cox模型的有约束条件优化问题转化为无条件优化问题，从而减低了极大化其完全似然函数的难度。再者，极大化该完全似然函数时没有解析解，且涉及到非参数估计，因此我们把具有参数分离特点的MM算法应用到极大化过程中，该算法把非参数部分从回归参数中分离出来，从而逐一击破。此外，在处理高维数据集时，MM算法可以把回归参数彼此分离，实现把高维优化问题转化为低维优化问题，从而在Newton-Raphson迭代过程中避免了矩阵求逆。模拟实验的估计结果表明，有约束条件的模型估计结果更可靠，且进一步表明ADMM和MM算法相结合具有良好的收敛性，且通过增加样本量、减少删失率，可以改善估计结果。而对真实数据集，即乳腺癌数据的分析，则进一步表明考虑约束条件对结果的估计是更可靠的。

NOTES

^*通讯作者。