生存分析主要探究影响事件发生时间的暴露因素, 是统计研究的热点领域之一.比例风险模型是研究生存数据中应用最为广泛的半参数模型之一.生存分析涉及领域非常广泛, 例如, 生物学、医学、金融学和管理学等领域.对比例风险模型中回归参数的估计常常基于其偏似然函数^[1], 通过求解相应的得分方程得到.由于此得分方程没有解析解, 在实际应用中, 需要通过数值方法来计算待估参数的数值解.Newton-Raphson算法作为最常用的数值计算方法之一, 有着诸多优点, 例如:在最优值附近有很快的二次收敛速度; 在求解过程中可以同时给出渐近方差矩阵等.但是, Newton-Raphson算法在求解过程中也常常面临一些问题, 例如:算法涉及二阶导函数矩阵求逆; 可能收敛到局部极值点或者鞍点; 可能不收敛等等^[2].因此, 寻找更加稳健的数值计算方法具有重要的理论意义和应用价值.

近年来, 随着深度学习与机器学习领域研究的不断深入, 自适应优化算法得到了迅速发展.自适应优化算法属于梯度下降优化算法, 其主要特点是:在迭代过程中引入学习率作为迭代步长, 可沿着负梯度方向不断迭代以接近最优值点, 因而能够在较大程度上节省计算成本和计算空间.梯度下降优化算法的缺点在于收敛速度较慢, 早期算法如BGD算法、SGD算法等采用固定学习率作为迭代步长, 主要有两方面的问题:一是学习率的设定及调整问题, 早期算法收敛速度过度依赖学习率的设定.学习率过小导致收敛速度过慢, 学习率过大导致在到达最优点后过度优化, 进而收敛到次优点或者局部最优点.二是迭代过程对所有参数采用相同的学习率, 当数据充分稀疏且特征分布不均时, 相同学习率的算法收敛效率和收敛性较差, 在数据分析中被证明是不可取的(见文献[3]).

针对早期算法的问题, 人们对自适应优化算法的探索主要基于两个方向:一个方向是基于动量方向的优化算法.文献[4]提出Momentum算法, Momentum算法引入动量的概念, 将当期梯度与前期梯度进行加权求和, 如果两期梯度方向一致, 则加权求和的过程不断累加, 累加过程使得学习率不再依赖于初始学习率的设置.后期发展的Nesterov算法(Nesterov Momentum) (见文献[5]), RMSprop算法(Root Mean Square Prop) (见文献[6])都是在Momentum算法的基础上进行的改进.

另一个方向是基于$ L_2 $范数下的优化算法.文献[3]提出了Delta-bar-Delta方法, 该方法基于一个较为明显的想法:在给定模型的条件下, 当损失函数对于某参数的偏导持续保持相同的符号, 则增加其更新步长, 反之减少. Delta-bar-Delta方法是最早自适应优化算法的尝试.由此文献[7]提出了Adagrad算法(Adaptive Gradient Algorithm), Adagrad算法引进累积梯度平方对学习率进行缩放, 从而为每个参数提供了与自身相关的学习率.然而, 由于累积梯度平方随着迭代过程递增, Adagrad算法存在后期学习率衰减的问题.因此, 文献[6]提出了RMsprop算法(Root Mean Square Prop), RMsprop算法对累积梯度平方进行了指数滑动平均以丢弃较远的历史数据, 解决了Adagrad算法的学习率衰减问题.进而, 文献[8]在RMSprop算法的基础上提出了Adam算法(Adaptive Moment Estimation), Adam算法将指数滑动平均应用到梯度的一阶矩估计当中, 使参数估计更具鲁棒性, 从而实现更快的收敛速度.

据我们所知, 将自适应优化算法应用到比例风险模型的研究尚且较少.本文将探究比例风险模型下计算参数估计的自适应优化算法及其改进算法, 从分析自适应优化算法在梯度下降良好的收敛性出发, 使用Adagrad算法、RMsprop算法、Adam算法及其改进算法替代Newton-Raphson算法应用于比例风险模型参数估计的求解问题, 探究其合理性和优良性.在下节中, 先介绍三种自适应优化算法:Adagrad算法、RMSprop算法与Adam算法.

本节介绍比例风险模型下, 三种自适应优化算法—Adagrad算法、RMSprop算法、Adam算法求解回归参数的极大似然估计数值解的步骤, 然后, 对自适应优化算法中Adam算法进行改进, 使其具有更快的收敛速度.

在生存过程中, 用$ \widetilde{T_i} $表示第$ i $个个体的死亡时间, $ C_i $表示第$ i $个个体的删失时间, $ T_i = $ Min$ (\widetilde{T_i}, C_i) $表示第$ i $个个体的观测时间, $ \triangle_i = I (T_i<C_i) $表示右删失示性变量, $ Z(t) $表示$ t $时刻的协变量.文献[1]提出如下比例风险率模型

$ \begin{equation} \lambda(t|Z) = \lambda_0(t)\exp\left(\beta 'Z(t)\right), \end{equation} $

(2.1)

其中$ \lambda_0(t) $是基准风险率函数, $ \beta $是$ p $维待估参数, 其参数空间为$ {B} $.对参数$ \beta $的估计方法, 常见的有矩法估计、最小二乘估计、极大似然似然估计、贝叶斯估计等.在极大似然估计方法下, 文献[1]给出对$ \beta $的统计推断的偏似然函数

$ L(\beta ) = \prod\limits_{i = 1}^n {{{[\frac{{{e^{\beta '{Z_i}({T_i})}}}}{{\sum\limits_{l \in R({T_i})} {{e^{\beta '{Z_l}({T_i})}}} }}]}^{\triangle_i}}} $

(2.2)

其中$ R(t) = {\left\{j:T_j\geq{t}\right\}} $表示$ t $时刻的风险集, 通常对似然函数取对数

$ \begin{equation} l(\beta) = \sum\limits_{i = 1}^n\triangle_i[\beta'Z_i(T_i)-\log\{\sum\limits_{l\in{R(T_i)}}e^{\beta'Z_l(T_i)}\}], \end{equation} $

(2.3)

因此$ \beta $的极大似然估计定义为

$ \begin{equation} \hat{\beta} = \text{arg}\max\limits_{\beta\in{B}}l(\beta). \end{equation} $

(2.4)

极大似然估计具有较好的估计性质, 文献[9]给出了极大似然估计的相合性证明以及渐近正态性证明.在实际问题求解过程中, 一般是将极大似然估计的极值问题(2.4)转换为得分方程求解问题

$ \begin{equation} \triangledown l(\beta) = \sum\limits_{i = 1}^n\triangle_i [Z_i(T_i)-\frac{\sum\limits_{l\in R(T_i)}{Z_l(T_i)}e^{\beta'Z_l(T_i)}} {\sum\limits_{l\in R(T_i)}e^{\beta'Z_l\left(T_i\right)}}] = 0. \end{equation} $

(2.5)

由于此得分方程没有解析解, 在实际问题中, 需要采用数值计算方法求解. Newton-Raphon算法是常用的数值计算方法之一, 在求解过程中需要计算二阶导函数矩阵$ \triangledown^2 l(\beta_t) $, 其形式如下

$ \begin{equation} \triangledown^2 l(\beta_t) = \sum\limits_{i = 1}^{n}\triangle_i [\frac{\sum\limits_{l\in R(T_i)}Z(T_i)^{\otimes2}e^{\beta'_tZ_l(T_i)}} {\sum\limits_{l\in R(T_i)}e^{\beta'_tZ_l(T_i)}} -(\frac{\sum\limits_{l\in R(T_i)}Z_le^{\beta'_tZ_l(T_i)}} {\sum\limits_{l\in R(T_i)}e^{\beta'_tZ_l(T_i)}} ) ^{\otimes2} ], \end{equation} $

(2.6)

其中$ \otimes $表示向量间的运算, 对于向量$ x $, 定义$ x^{\otimes2} = xx' $.

因此, 得出比例风险模型下, Newton-Raphon算法的具体步骤如下

步骤$ 1 $ 给定初值$ \beta_0 $以及计算精度$ \epsilon $;

步骤$ 2 $ 基于如下迭代公式, 将第$ t-1 $步估计值$ \beta_{t-1} $更新为第$ t $步估计值$ \beta_{t} $

$ \begin{equation} \begin{aligned} \beta_{t} = \beta_{t-1}-(\triangledown^2 l(\beta_t))^{-1} \cdot \triangledown l(\beta_t); \end{aligned} \end{equation} $

(2.7)

步骤$ 3 $ 停止准则$ \parallel{\beta_{t}-\beta_{t-1}}\parallel<\epsilon $, 即:不满足时, 令$ t = t+1 $, 返回步骤$ 2 $; 满足时, 则输出估计值$ \beta_{t} $, 迭代结束.

在实际问题求解中, 二阶导函数矩阵可能是奇异矩阵, 无法求逆进而导致Newton-Raphon算法迭代终止.因而, 需要寻找更稳定的数值计算方法.

首先引入Adagrad算法进行求解, Adagrad算法核心思想是以累积梯度平方$ R_{t} $对迭代步长进行缩小或放大, 即将迭代公式(2.7)调整为

$ \begin{equation} \beta_{t} = \beta_{t-1}- \frac{\alpha}{\sqrt{R_t}+\tau} \cdot \triangledown l(\beta_t), \end{equation} $

(2.8)

其中$ \tau $为平滑常数, 其目的是为了防止出现分母为$ 0 $的情况, 一般取为较小的常数.其中累积梯度平方$ R_{t} $的形式为

$ \begin{equation} \begin{aligned} R_1& = \triangledown l(\beta_1)^{\bigodot2}, \\ R_{t}& = R_{t-1}+\triangledown l(\beta_t)^{\bigodot2}, \end{aligned} \end{equation} $

(2.9)

其中$ \bigodot $表示元素级别的运算, 对于向量$ x $, $ x^{\bigodot2} $即对向量里的每个元素进行平方.

由(2.9)式不难发现累积梯度平方$ R_t $是递增的, 这将造成学习率不断递减.因而引入RMsprop算法. RMSprop算法利用梯度平方的指数滑动平均过程替换()式的累加过程, 并引入滑动系数$ \phi $控制指数滑动平均的衰减速率, $ \phi \in [0, 1) $,

$ \begin{equation} \begin{aligned} \bar{R}_1& = (1-\phi)\triangledown l(\beta_1)^{\bigodot2}, \\ \bar{R}_t& = \phi \bar{R}_{t-1}+(1-\phi)\triangledown l(\beta_t)^{\bigodot2}. \end{aligned} \end{equation} $

(2.10)

由上可知, $ \bar{R}_t $是梯度平方$ \triangledown l(\beta)^{\bigodot2} $的指数滑动平均值. RMSprop算法利用指数滑动平均值$ \bar{R}_t $对迭代步长进行调整, 将(2.7)式调整为

$ \begin{equation} \beta_{t} = \beta_{t-1}- \frac{\alpha}{\sqrt{\bar{R}_t}+\tau} \cdot \triangledown l(\beta_t). \end{equation} $

(2.11)

尽管RMSprop算法在求解极值问题时已经达到较快的收敛速度, 但是为了获取更快的收敛速度, 引入Adam算法进行求解.在Adam算法下, 定义梯度$ \triangledown l(\beta) $的有偏一阶矩估计为

$ \begin{equation} \begin{aligned} {m}_1& = (1-\psi_1)\triangledown l(\beta_1), \\ {m}_t& = \psi_1m_{t-1}+(1-\psi_1)\triangledown l(\beta_t), \end{aligned} \end{equation} $

(2.12)

其中$ \psi_1 $为控制梯度$ \triangledown l(\beta) $指数滑动平均过程中衰减速率的超参数, $ \psi_1 \in [0, 1) $默认取值$ 0.9 $.

同时定义梯度$ \triangledown l(\beta) $的有偏二阶矩估计为

$ \begin{equation} \begin{aligned} {v}_1& = (1-\psi_2)\triangledown l(\beta_1)^{\bigodot2}, \\ {v}_t& = \psi_2v_{t-1}+(1-\psi_2)\triangledown l(\beta_t)^{\bigodot2}, \end{aligned} \end{equation} $

(2.13)

其中$ \psi_2 $为控制梯度平方$ \triangledown l(\beta)^{\bigodot2} $指数滑动平均过程中衰减速率的超参数, $ \psi_2 \in [0, 1) $默认取值0.999.

由(2.12)式与(2.13)式不难发现, 在$ \psi_1 $, $ \psi_2 $为默认取值的情况下, 在迭代初期, $ {m}_t $, $ {v}_t $是偏向于$ 0 $的, 因此需要进行修正

$ \begin{equation} \begin{aligned} \hat{m}_t& = {m_t}/{(1-\psi_1^t)}, \\ \hat{v}_t& = {v_t}/{(1-\psi_2^t)}. \end{aligned} \end{equation} $

(2.14)

进而, 将(2.7)式替换为

$ \begin{equation} \beta_{t} = \beta_{t-1}- \frac{\alpha}{\sqrt{\hat{v}_t}+\tau} \cdot \hat{m}_t. \end{equation} $

(2.15)

由此, 得出比例风险模型下, Adam算法的具体步骤如下

步骤$ 1 $ 给定初值$ \beta_0 $以及计算精度$ \epsilon $;

步骤$ 2 $ 计算梯度的一阶矩估计$ \hat{m}_t $, 其具体形式见(2.12), (2.14)式;

步骤$ 3 $ 计算梯度的二阶矩估计$ \hat{v}_t $, 其具体形式见(2.13), (2.14)式;

步骤$ 4 $ 基于如下迭代公式, 将第$ t-1 $步估计值$ \beta_{t-1} $更新为第$ t $步估计值$ \beta_{t} $,

$ \begin{equation} \beta_{t} = \beta_{t-1}- \frac{\alpha}{\sqrt{\hat{v}_t}+\tau} \cdot \hat{m}_t; \end{equation} $

(2.16)

步骤$ 5 $ 停止准则$ \parallel{\beta_{t}-\beta_{t-1}}\parallel<\epsilon $, 即:不满足时, 令$ t = t+1 $, 返回步骤$ 2 $; 满足时, 则输出估计值$ \beta_{t} $, 迭代结束.

在上一小节中, 介绍了Adam算法在比例风险模型下的求解步骤, 在给定初始学习率$ \alpha $, 平滑常数$ \tau $时, Adam算法在第$ t $步的迭代步长为

$ \begin{equation} \triangle\beta_{\text{Adam}} = \alpha \cdot \frac{\hat{m}_t}{\sqrt{\hat{v}_t}+\tau}. \end{equation} $

(2.17)

结合(2.13)式对$ \hat{v}_t $进行逐项展开

$ \begin{equation} \begin{aligned} \hat{v}_t& = \frac{1}{1-\psi_2^t}\left[ \psi_2v_{t-1}+(1-\psi_2)\triangledown l(\beta_t)^{\bigodot2}\right]\\ & = \frac{1}{1-\psi_2^t}\left[ \psi_2^2v_{t-2}+\psi_2(1-\psi_2)\triangledown l(\beta_{t-1})^{\bigodot2}+(1-\psi_2)\triangledown l(\beta_t)^{\bigodot2} \right]\\ & = \frac{1}{1-\psi_2^t}\left[ \psi_2^kv_{t-k}+\psi_2^{k-1}(1-\psi_2)\triangledown l(\beta_{t-k+1})^{\bigodot2}+\cdots+(1-\psi_2)\triangledown l(\beta_t)^{\bigodot2} \right], \end{aligned} \end{equation} $

(2.18)

其中$ k = 1, 2, \cdots, t $.由上式不难发现, 指数滑动平均对近期梯度平方赋予更高的权重, 而远期梯度平方的权重则以指数形式衰减, 从而有效控制了梯度平方的累积过程.然而, 这种梯度平方的短期记忆能力很可能造成算法收敛到次优的极值点.发现与文献[10]在ICLR $ 2018 $中正在审核的一篇文章中提出的问题不谋而合.

为了解决梯度平方的短期记忆问题, 引入记号$ P^{A}_t $表示Adam算法第$ t $步迭代的更新量

$ \begin{equation} P^{A}_t = \frac{\hat{m}_t}{\sqrt{\hat{v}_t}+\tau}. \end{equation} $

(2.19)

注意到更新量$ P^{A}_t $与Adam算法的迭代步长(2.17)式只差一个常数——学习率$ \alpha $.结合经典算法Momentum算法(见文献[4]中动量向前的思想, 动量向前的思想是模拟物理里动量的概念, 用积累之前的动量来替代真正的梯度, 记动量符号为$ M $ (Momentum)

$ \begin{equation} \begin{aligned} M_1& = 0, \\ M_t& = u\cdot {M_{t-1}}+\triangledown l(\beta), \end{aligned} \end{equation} $

(2.20)

其中$ u $为动量$ M $所占的权重.

记Adam算法的改进算法为MAdam算法(Modified Adam), 记MAdam算法的第$ t $步的迭代步长为$ P_t^M $.结合动量向前思想(2.20)式, 使用MAdam算法第$ t-1 $步的迭代步长与Adam算法第$ t $步的迭代更新量进行比较, 取二者之间的最大值, 与Adam算法第$ t $步的迭代更新量进行加权求和, 作为MAdam算法第$ t $步的迭代步长.其具体形式为

$ \begin{equation} \begin{aligned} P_1^M& = \eta \cdot P_1^A, \\ P_t^M& = \alpha \cdot {\hbox{Max}}(P_{t-1}^M, P_t^A)+\eta \cdot P_t^A, \end{aligned} \end{equation} $

(2.21)

其中$ \alpha $为学习率, $ \eta $为动量因子, 即Adam算法更新量$ P_t^A $的权重, $ \eta\in [0, 1) $.

结合(2.16)式, 将Adam算法的迭代过程修正为

$ \begin{equation} \beta_t = \beta_{t-1}-P_t^M. \end{equation} $

(2.22)

由于$ P_t^M $是$ P_{t-1}^M $与$ P_t^A $最大值与$ P_t^A $的加权求和, 因此不难发现在迭代过程中$ P_t^M $是以上二者最大值与$ P_t^A $的排列组合.因此分别分析最大值选取的两种情形下的收敛情况, 当$ P_t^A\geq P_{t-1}^M $时, MAdam算法退化为学习率为$ (\alpha+\eta) $的Adam算法

$ \begin{equation} \begin{aligned} P_1^M& = \eta \cdot P_1^A, \\ P_t^M& = (\alpha+\eta )\cdot P_t^A. \end{aligned} \end{equation} $

(2.23)

当$ P_t^A < P_{t-1}^M $时, MAdam算法则演变为衰减速率为$ \alpha $的Adam算法更新量的指数滑动平均过程

$ \begin{equation} \begin{aligned} P_1^M& = \eta \cdot P_1^A, \\ P_t^M& = \alpha \cdot P_{t-1}^M +\eta \cdot P_t^A. \end{aligned} \end{equation} $

(2.24)

由(2.24)式, 进行如下推导

$ \begin{equation} \begin{aligned} P_t^M& = \alpha \cdot P_{t-1}^M +\eta \cdot P_t^A, \\ P_t^M& = \sum\limits_{j = 1}^{t}\alpha^{t-j}\eta \cdot P_j^A. \end{aligned} \end{equation} $

(2.25)

结合(2.23)式与(2.25)式, 可知MAdam算法与Adam算法的关系为

$ \begin{equation} \begin{aligned} P_1^M& = \eta \cdot P_1^A, \\ P_t^M& = \begin{cases}{(\alpha+\eta )\cdot P_t^A}, & P_t^A \geq P_{t-1}^M, \\ {\sum\limits_{j = 1}^{t}\alpha^{t-j}\eta \cdot P_j^A}, & P_t^A < P_{t-1}^M.\end{cases} \end{aligned} \end{equation} $

(2.26)

由此, 得到比例风险模型下, MAdam算法的具体步骤如下

步骤$ 1 $ 给定初值$ \beta_0 $以及计算精度$ \epsilon $;

步骤$ 2 $ 计算梯度的一阶矩估计$ \hat{m}_t $, 其具体形式见(2.12), (2.14)式;

步骤$ 3 $ 计算梯度的二阶矩估计$ \hat{v}_t $, 其具体形式见(2.13), (2.14)式;

步骤$ 4 $ 计算Adam更新量:$ P_{t}^{A} = \frac{\hat{m}_t}{\sqrt{\hat{v}_t}+\tau} $;

步骤$ 5 $ 计算MAdam更新量:$ P_{t}^{M} = \alpha \cdot {\hbox{Max}}(P_{t}^{A}, P_{t-1}^{M})+\eta \cdot P_{t}^{A} $;

步骤$ 6 $ 基于如下迭代公式, 将第$ t-1 $步估计值$ \beta_{t-1} $更新为第$ t $步估计值$ \beta_{t} $:

$ \begin{equation} \beta_{t} = \beta_{t-1}- P_{t}^{M}; \end{equation} $

(2.27)

步骤$ 7 $ 停止准则$ \parallel{\beta_{t}-\beta_{t-1}}\parallel<\epsilon $, 即:不满足时, 令$ t = t+1 $, 返回步骤$ 2 $; 满足时, 则输出估计值$ \beta_{t} $, 迭代结束.

本节中, 通过数据模拟展示四种自适应优化算法替代Newton-Raphon算法的可行性, 展示在不同数据样本量以及不同删失率下, 各优化算法的数值计算时间.

在给定协变量$ Z_1 $, $ Z_2 $的条件下, 考虑失效时间$ \tilde{T} $的风险率函数服从如下形式的比例风险模型

$ \lambda(t|Z_1, Z_2) = \lambda_0(t)\text{exp}(\beta_1Z_1+\beta_2Z_2). $

为了比较各种参数设定下, 自适应优化算法与Newton-Raphon算法的参数估计情况, 设定参数真值为$ \beta_1 = 0.693 $, $ \beta_2 = -0.5 $.协变量$ Z_1 $从成功概率为$ 0.5 $的Bernouli分布中随机生成; 协变量$ Z_2 $从标准正态分布$ N(0, 1) $中随机生成.删失时间$ C $服从均匀分布$ U(0, c) $, 设定不同$ c $值, 使得$ 1000 $次模拟产生相应的删失率分别为$ \rho = 0.2 $, $ 0.5 $或$ 0.7 $.基准风险率$ \lambda_0(t) $分别设定为$ 1 $和$ 2t $, 由此可分别由风险率为$ \text{exp}{\left\{\beta_1Z_1+\beta_2Z_2\right\}} $的指数分布数据, 以及形状参数为$ 2 $, 尺度参数为$ 1/\sqrt{\text{exp}{\left\{\beta_1Z_1+\beta_2Z_2\right\}}} $的Weibull分布中随机生成失效时间$ \tilde{T} $.样本量$ n $设置为$ 100 $, $ 200 $, $ 250 $或$ 300 $.

模拟结果中, Newton-Raphon算法、Adam算法、MAdam算法、RMSprop算法以及Adagrad算法参数估计的估计均值(Mean)、参数估计的样本标准差(SD), 标准误差的估计均值(SE)和$ 95\% $置信区间覆盖率(CP)均由$ 1000 $次独立模拟结果获得.迭代误差$ \epsilon $设置为$ 10^{-5} $.算法的数值计算时间(Time)为$ 1000 $次模拟中平均每次计算所耗费的系统计算时间(单位:秒), 使用的计算系统配置为:CPU型号为Intel Xeon E5-2630 v2, 12核, CPU主频为$ 2.6 $ GHz, 内存大小为$ 64 $ GB.

Adam算法、MAdam算法、RMSprop算法以及Adagrad算法的初始学习率设置为$ \alpha = 0.01 $, 平滑常数设置为$ \tau = 10^{-8} $.RMSprop算法的滑动系数使用默认值$ \phi = 0.5 $. Adam算法的超参数使用默认值$ \psi_1 = 0.9 $, $ \psi_1 = 0.999 $, MAdam算法的超参数与Adam算法保持一致, 权重$ \eta = 0.05 $.数值模拟计算结果请见表 1– 表 4, 各算法计算时间, 请见图 1–图 6.

从表 1–表 4的数值计算结果可知, 在以上所有考虑的参数设定情况下, Adam算法、MAdam算法、RMSprop算法得到的极大似然估计均是无偏的, 随着样本量增大, 标准误差的估计均值SE和参数估计的样本标准差SD逐渐接近, 置信区间覆盖率CP也稳定在$ 0.95 $左右, 较为合理.除了Adagrad算法计算偏差过大以外, 其他三种算法得到的参数估计结果与Newton-Raphon算法估计结果一致, 这说明自适应优化算法替代Newton-Raphon算法的可行性.在实际应用中, 当Newton-Raphon算法出现无法求逆的情况下, 以上三种自适应优化算法是替代Newton-Raphon算法的可行方案之一.

进一步, 设置样本量$ n = 50 $, $ 100 $, $ 150 $, $ 200 $, $ 250 $, $ 300 $, 设定删失率分别为$ \rho = 0.2 $, $ 0.5 $, $ 0.7 $.在计算时间上, 自适应优化算法当中, MAdam算法的计算效率均高于其他三种自适应优化算法, 这说明对Adam算法的改进是可行的, 并有效地缩短了计算时间.以上几种自适应优化算法在收敛速度上与Newton-Raphon算法相比较慢.但是, 自适应优化算法能克服很多高维数据计算时的困难, 具有其计算优势和应用价值.

例$ 1 $ (癌症临床实验数据案例)首先分析一个癌症临床实验的实际案例, 数据来源:CluBear的公开数据(熊赳赳教育科技(北京)有限公司).该数据集合收集了$ 120 $个参与某临床试验的癌症病例, 目的是为了发现患者生存时间的相关影响因素, 并对比评估某种新治疗方案的疗效(同标准治疗方案对比).

该数据集随机抽取某临床试验的$ 120 $个病人, 因变量是病人的生存时间(单位:天).同时数据还提供该病人在试验结束时的生存状态.那些仍然生存的病人的真实生存时间未知, 因此他们的被记录的生存时间是删失的, 用右删失示性变量$ \delta = 0 $表示存活, $ \delta = 1 $表示死亡.解释变量有:不同治疗方案协变量$ Z_1 = 0 $表示标准治疗方案, $ Z_1 = 1 $表示新治疗方案; 癌细胞类型协变量$ Z_2 = 1, 2, 3, 4 $分别表示A、B、C、D四种类型的癌细胞; 主治医师临床打分$ Z_3 $表示医生对病人的综合打分, 连续型变量; 病人的年龄$ Z_4 $(单位:岁).

建立建立如下比例风险模型评估新治疗方案与生存时间之间的关系

$ \lambda(t|Z_1, Z_2, Z_3, Z_4) = \lambda_0(t)\text{exp}(\beta_1Z_1+\beta_2Z_2+\beta_3Z_3+\beta_4Z_4). $

模型参数估计的过程使用Newton-Raphon算法、Adam算法以及MAdam算法得出参数估计结果, 请见表 6.可以看出Adam算法以及MAdam算法与Newton-Raphon算法的参数估计结果是一致的.从数值估计结果来看, 新旧方案没有明显的差别, 协变量$ Z_1 $的检验$ p $值为$ 0.389 $.综上结果表明, 对该癌症案例具有影响的是癌细胞类型与主治医师临床打分, 没有充分的证据可以认定该癌症案例的治疗新方案能够有效的降低患者死亡的风险.

例$ 2 $(天津空气质量数据案例)本小节, 分析一组空气质量数据, 数据来源:CluBear的公开数据(熊赳赳教育科技(北京)有限公司).近年来, 城市雾霾天气在我国频繁出现, 空气质量关乎人民日常正常的生活与劳动, 因此大中型城市空气质量问题已经引起全社会高度关注.在一份亚洲开发银行与清华大学联合发布的名为《迈向环境可持续的未来中华人民共和国国家环境分析》(见文献[11])报告中指出, 尽管我国政府一直积极运用财政和行政手段治理大气环境污染, 但我国仍是世界上空气污染最严重国家之一.在我国空气污染最严重的城市排名中, 天津位列前茅.

本文获取了$ 1809 $条数据, 结构如下:因变量为空气质量指数(AQI); 协变量$ Z_1 $表示PM$ 2.5 $指数; $ Z_2 $表示空气质量的$ 7 $个水平:轻度污染、中度污染、重度污染、严重污染、良、优、无; $ Z_3 $表示PM$ 10 $指数; $ Z_4 $表示空气中SO$ _2 $含量; $ Z_5 $表示空气中Co含量; $ Z_6 $表示空气中NO$ _2 $含量; $ Z_7 $表示空气中O$ 3\_8 $h含量.建立如下比例风险模型来探究AQI影响因素

$ \lambda(t|Z) = \lambda_0(t)\text{exp}(\beta_1Z_1+\beta_2Z_2+\beta_3Z_3+\beta_4Z_4+\beta_5Z_5+\beta_6Z_6+\beta_7Z_7). $

用Newton-Raphon算法、Adam算法以及MAdam算法进行参数估计.其余变量均对AQI具有显著影响.

通过建立比例风险模型, 分析了以上各协变量对空气质量AQI的影响, 模型参数的数值计算结果请见表 7.发现, 三种算法的参数估计结果一致, 与Adam算法相比, MAdam算法的估计结果更接近Newton-Raphon算法, 并且计算速度更快.当待估参数维度增加时, Newton-Raphon算法涉及高维矩阵的求逆问题.因此, 相对于Newton-Raphon算法, 自适应优化算法在高维参数的估计问题上具有一定的计算优势.