加速失效时间模型的生存分析

什么是生存分析？

生存分析（回归）建模的是**事件发生的时间**。生存分析是一种特殊的回归分析，与传统回归任务的区别如下：

  • 标签值始终为正，因为事件发生所需的时间不可能是负数。

  • 标签值可能不完全已知，或者称为**截尾（censored）**，因为“测量时间需要时间”。

第二点至关重要，我们应该深入探讨。顾名思义，生存分析最早的应用之一是建模给定人群的死亡率。我们以NCCTG 肺癌数据集为例。前 8 列代表特征，最后一列“死亡时间”代表标签。

机构

年龄

性别

ph.ecog

ph.karno

pat.karno

meal.cal

wt.loss

死亡时间（天）

3

74

1

1

90

100

1175

不可用

306

3

68

1

0

90

90

1225

15

455

3

56

1

0

90

90

不可用

15

\([1010, +\infty)\)

5

57

1

1

90

60

1150

11

210

1

60

1

0

100

90

不可用

0

883

12

74

1

1

50

80

513

0

\([1022, +\infty)\)

7

68

2

2

70

60

384

10

310

仔细看第三位患者的标签。他的标签是一个范围，而不是一个单一的数字。第三位患者的标签被称为**截尾**，因为实验人员出于某种原因未能获得该标签的完整测量值。一种可能的情况是：患者在前 1010 天幸存下来，并在第 1011 天离开了诊所，因此未能直接观察到他的死亡。另一种可能性是：实验在他死亡被观察到之前被提前结束（因为实验不能永远进行）。无论如何，他的标签是\([1010, +\infty)\)，这意味着他的死亡时间可以是任何大于 1010 的数字，例如 2000、3000 或 10000。

截尾有四种类型

  • 未截尾：标签未截尾，给出一个单一的数字。

  • 右截尾：标签形式为\([a, +\infty)\)，其中\(a\)是下界。

  • 左截尾：标签形式为\([0, b]\)，其中\(b\)是上界。

  • 区间截尾：标签形式为\([a, b]\)，其中\(a\)和\(b\)分别是下界和上界。

右截尾是最常用的类型。

加速失效时间模型

加速失效时间 (AFT) 模型是生存分析中最常用的模型之一。该模型形式如下：

\[\ln{Y} = \langle \mathbf{w}, \mathbf{x} \rangle + \sigma Z\]

其中

  • \(\mathbf{x}\)是\(\mathbb{R}^d\)中的向量，代表特征。

  • \(\mathbf{w}\)是由\(d\)个系数组成的向量，每个系数对应一个特征。

  • \(\langle \cdot, \cdot \rangle\)是\(\mathbb{R}^d\)中的标准点积。

  • \(\ln{(\cdot)}\)是自然对数。

  • \(Y\)和\(Z\)是随机变量。

    • \(Y\)是输出标签。

    • \(Z\)是一个已知概率分布的随机变量。常见的选择包括正态分布、逻辑斯谛分布和极值分布。直观上，\(Z\)代表将预测值\(\langle \mathbf{w}, \mathbf{x} \rangle\)从真实的对数标签值\(\ln{Y}\)拉开的“噪声”。

  • \(\sigma\)是一个用于缩放\(Z\)大小的参数。

注意，该模型是线性回归模型\(Y = \langle \mathbf{w}, \mathbf{x} \rangle\)的广义形式。为了使 AFT 模型能够与梯度提升配合使用，我们对模型进行如下修改：

\[\ln{Y} = \mathcal{T}(\mathbf{x}) + \sigma Z\]

其中\(\mathcal{T}(\mathbf{x})\)表示给定输入\(\mathbf{x}\)时，决策树集成模型的输出。由于\(Z\)是一个随机变量，我们可以为表达式\(\ln{Y} = \mathcal{T}(\mathbf{x}) + \sigma Z\)定义一个似然函数。因此，XGBoost 的目标是通过拟合一个好的树集成模型\(\mathcal{T}(\mathbf{x})\)来最大化（对数）似然函数。

如何使用

第一步是将标签表示为范围的形式，以便每个数据点都有与之相关的两个数字，即标签的下界和上界。对于未截尾的标签，使用退化的区间形式\([a, a]\)。

截尾类型

区间形式

下界有限？

上界有限？

未截尾

\([a, a]\)

✔

✔

右截尾

\([a, +\infty)\)

✔

✘

左截尾

\([0, b]\)

✔

✔

区间截尾

\([a, b]\)

✔

✔

将下界数值收集到一个数组（称之为y_lower_bound）中，并将上界数值收集到另一个数组（称之为y_upper_bound）中。通过调用xgboost.DMatrix.set_float_info()，可以将范围标签与数据矩阵对象关联起来。

Python
import numpy as np
import xgboost as xgb

# 4-by-2 Data matrix
X = np.array([[1, -1], [-1, 1], [0, 1], [1, 0]])
dtrain = xgb.DMatrix(X)

# Associate ranged labels with the data matrix.
# This example shows each kind of censored labels.
#                         uncensored    right     left  interval
y_lower_bound = np.array([      2.0,     3.0,     0.0,     4.0])
y_upper_bound = np.array([      2.0, +np.inf,     4.0,     5.0])
dtrain.set_float_info('label_lower_bound', y_lower_bound)
dtrain.set_float_info('label_upper_bound', y_upper_bound)
R
library(xgboost)

# 4-by-2 Data matrix
X <- matrix(c(1., -1., -1., 1., 0., 1., 1., 0.),
            nrow=4, ncol=2, byrow=TRUE)
dtrain <- xgb.DMatrix(X)

# Associate ranged labels with the data matrix.
# This example shows each kind of censored labels.
#                   uncensored  right  left  interval
y_lower_bound <- c(        2.,    3.,   0.,       4.)
y_upper_bound <- c(        2.,  +Inf,   4.,       5.)
setinfo(dtrain, 'label_lower_bound', y_lower_bound)
setinfo(dtrain, 'label_upper_bound', y_upper_bound)

现在我们可以调用训练 API 了

Python
params = {'objective': 'survival:aft',
          'eval_metric': 'aft-nloglik',
          'aft_loss_distribution': 'normal',
          'aft_loss_distribution_scale': 1.20,
          'tree_method': 'hist', 'learning_rate': 0.05, 'max_depth': 2}
bst = xgb.train(params, dtrain, num_boost_round=5,
                evals=[(dtrain, 'train')])
R
params <- list(objective='survival:aft',
               eval_metric='aft-nloglik',
               aft_loss_distribution='normal',
               aft_loss_distribution_scale=1.20,
               tree_method='hist',
               learning_rate=0.05,
               max_depth=2)
watchlist <- list(train = dtrain)
bst <- xgb.train(params, dtrain, nrounds=5, watchlist)

我们将objective参数设置为survival:aft，将eval_metric设置为aft-nloglik，以便最大化 AFT 模型的对数似然。(XGBoost 实际上会最小化负对数似然，因此得名aft-nloglik。)

参数aft_loss_distribution对应于 AFT 模型中\(Z\)项的分布，而aft_loss_distribution_scale对应于缩放因子\(\sigma\)。

目前，您可以从以下三种概率分布中选择aft_loss_distribution：

aft_loss_distribution

概率密度函数 (PDF)

normal

\(\dfrac{\exp{(-z^2/2)}}{\sqrt{2\pi}}\)

logistic

\(\dfrac{e^z}{(1+e^z)^2}\)

extreme

\(e^z e^{-\exp{z}}\)

请注意，目前尚无法使用 scikit-learn 接口（例如 xgboost.XGBRegressor）设置范围标签。目前，您应该使用 xgboost.train 配合 xgboost.DMatrix 进行训练。有关 Python 示例的集合，请参阅 生存分析示例演练