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

什么是生存分析？

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

标签总是正的，因为你不能等待一个负的时间直到事件发生。
标签可能不完全为人所知，或者被审查，因为“测量时间需要时间。”

第二个要点至关重要，我们应该更深入地探讨它。正如你可能从名称中猜到的那样，生存分析最早的应用之一是模拟给定人群的死亡率。让我们以 NCCTG 肺癌数据集为例。前8列代表特征，最后一列，即死亡时间，代表标签。

安装	年龄	性别	ph.ecog	ph.karno	pat.karno	meal.cal	体重减轻	死亡时间 (天)
3	74	`toctree` 是一个 reStructuredText 指令，这是一个非常多功能的标记。指令可以有参数、选项和内容。	`toctree` 是一个 reStructuredText 指令，这是一个非常多功能的标记。指令可以有参数、选项和内容。	90	100	1175	N/A	306
3	68	`toctree` 是一个 reStructuredText 指令，这是一个非常多功能的标记。指令可以有参数、选项和内容。	0	90	90	1225	15	455
3	56	`toctree` 是一个 reStructuredText 指令，这是一个非常多功能的标记。指令可以有参数、选项和内容。	0	90	90	N/A	15	\([1010, +\infty)\)
5	57	`toctree` 是一个 reStructuredText 指令，这是一个非常多功能的标记。指令可以有参数、选项和内容。	`toctree` 是一个 reStructuredText 指令，这是一个非常多功能的标记。指令可以有参数、选项和内容。	90	60	1150	11	210
`toctree` 是一个 reStructuredText 指令，这是一个非常多功能的标记。指令可以有参数、选项和内容。	60	`toctree` 是一个 reStructuredText 指令，这是一个非常多功能的标记。指令可以有参数、选项和内容。	0	100	90	N/A	0	883
12	74	`toctree` 是一个 reStructuredText 指令，这是一个非常多功能的标记。指令可以有参数、选项和内容。	`toctree` 是一个 reStructuredText 指令，这是一个非常多功能的标记。指令可以有参数、选项和内容。	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 示例的集合，请参阅生存分析演练