SciPy - Cox Proportional Hazards Model 生存分析

Cox Proportional Hazards Model 是一种流行的用于生存分析的统计方法。它有助于估计各种变量对事件（如故障或死亡）发生所需时间的影响。该模型在处理截尾数据时特别有价值，在这种数据中，有些个体在研究结束时可能尚未经历该事件。

虽然 SciPy 没有内置的 Cox Proportional Hazards model，但基于 SciPy 的 lifelines 库为生存分析提供了全面支持，包括 Cox Proportional Hazards model。

使用 lifelines 在 Python 中实现 Cox 比例风险模型的步骤

以下是使用 lifelines 在 Python 中实现 Cox Proportional Hazards 的步骤 −

安装 Lifelines 库

首先，我们需要在命令提示符中执行以下命令来安装 lifelines 库（如果之前未安装） −

pip install lifelines

导入库

之后，我们需要导入所有必要的库 −

import pandas as pd
import numpy as np
from lifelines import CoxPHFitter

准备数据

对于 Cox Proportional Hazards 模型，我们的数据集应至少包含以下两个组成部分 −

• Duration, 即事件或删失发生的时间。
• Event/Censoring Indicator, 如果事件发生则为 1，如果观测被删失则为 0。

以下是实现 Cox Proportional Hazards Model 的示例数据集 −

# 示例数据集
data = {
    'age': [60, 65, 70, 80, 85],
    'sex': [1, 0, 1, 1, 0],  # 1 表示男性，0 表示女性
    'duration': [5, 6, 7, 8, 9],  # 年数作为 Duration
    'event': [1, 0, 1, 1, 0]  # 1 表示事件发生（例如死亡），0 表示删失
}
# 创建 DataFrame
df = pd.DataFrame(data)

拟合 Cox 比例风险模型

在这里，我们将显示模型摘要，包括每个预测变量的估计系数、标准误、z 分数和 p 值，这有助于我们理解它们对 hazard ratio 的影响。

# 实例化 Cox Proportional Hazards 模型
cph = CoxPHFitter()

# 使用数据集拟合模型
cph.fit(df, duration_col='duration', event_col='event')

# 打印模型摘要
cph.print_summary()

进行预测

模型拟合完成后，我们可以为新数据预测生存函数或计算累积 hazard。

# 为新个体预测生存函数
new_data = pd.DataFrame({
    'age': [75],
    'sex': [1],
})

# 为新个体预测生存函数
survival_function = cph.predict_survival_function(new_data)
print(survival_function)

示例

以下示例模拟了一个包含不同风险因素（如年龄、性别、吸烟状况）的数据集，并分析它们对生存时间的影响 −

import pandas as pd
from lifelines import CoxPHFitter

# 模拟数据集：包含年龄、性别和吸烟状况的生存数据
data = {
    'age': [60, 62, 65, 70, 72, 75, 80, 85],      # 个体年龄
    'sex': [1, 0, 1, 1, 0, 1, 0, 1],              # 1 表示男性，0 表示女性
    'smoking_status': [1, 0, 1, 0, 1, 0, 1, 0],   # 1 表示吸烟者，0 表示非吸烟者
    'duration': [5, 6, 7, 8, 9, 5, 6, 7],         # 生存时间（年）
    'event': [1, 1, 0, 1, 1, 0, 1, 1]             # 1 表示事件发生（例如死亡），0 表示删失
}

# 创建 DataFrame
df = pd.DataFrame(data)

# 实例化 Cox Proportional Hazards 模型
cph = CoxPHFitter()

# 使用数据集拟合模型
cph.fit(df, duration_col='duration', event_col='event')

# 打印模型摘要以分析结果
cph.print_summary()

# 为新个体（例如 70 岁吸烟者）预测生存函数
new_data = pd.DataFrame({
    'age': [70],
    'sex': [1],          # 男性
    'smoking_status': [1],  # 吸烟者
})

# 为新个体预测生存函数
survival_function = cph.predict_survival_function(new_data)

# 打印新个体的生存函数
print(survival_function)

以上示例的输出如下 −

< lifelines.CoxPHFitter: fitted with 8 total observations, 2 right-censored observations >
             duration col = 'duration'
                event col = 'event'
      baseline estimation = breslow
   number of observations = 8
number of events observed = 6
   partial log-likelihood = -7.39
         time fit was run = 2025-01-31 12:06:31 UTC

---
                coef exp(coef)  se(coef)  coef lower 95%  coef upper 95% exp(coef) lower 95% exp(coef) upper 95%
covariate
age            -0.02      0.98      0.06           -0.14            0.11                0.87                1.12
sex            -0.23      0.79      1.06           -2.31            1.84                0.10                6.32
smoking_status -0.55      0.58      1.03           -2.57            1.47                0.08                4.33

                cmp to     z    p  -log2(p)
covariate
age               0.00 -0.26 0.80      0.33
sex               0.00 -0.22 0.83      0.28
smoking_status    0.00 -0.54 0.59      0.76
---
Concordance = 0.47
Partial AIC = 20.79
log-likelihood ratio test = 0.33 on 3 df
-log2(p) of ll-ratio test = 0.07

            0
5.0  0.918415
6.0  0.734865
7.0  0.610552
8.0  0.435392
9.0  0.191870