在Python中实现Cox比例风险回归模型-牛站云

生存分析主要用于研究“某个事件发生之前需要多长时间”，例如患者从治疗开始到复发的时间、设备从投入运行到故障的时间，以及客户从注册到流失的时间。与普通回归模型相比，生存数据通常同时包含生存时间和删失信息，因此需要专门的统计方法进行处理。

Cox比例风险回归模型（Cox Proportional Hazards Model）是生存分析中应用非常广泛的一种方法。它既可以分析多个协变量对事件发生风险的影响，又不要求预先指定基准风险函数的具体分布，因此具有较强的实用性。

Python生态中可以使用 lifelines 等统计分析库快速完成Cox比例风险回归，包括数据准备、模型拟合、风险比计算、显著性检验以及比例风险假设检验等操作。

Cox比例风险模型的基本原理

Cox模型研究的是某一时刻个体发生事件的风险，即风险函数（Hazard Function）。

其经典形式为：


h(t|X) = h₀(t) × exp(β₁X₁ + β₂X₂ + ... + βₚXₚ)

其中：

h(t|X)：给定协变量条件下，在时间 t 发生事件的风险；
h₀(t)：基准风险函数；
X₁、X₂ ... Xₚ：解释变量或协变量；
β₁、β₂ ... βₚ：模型需要估计的回归系数。

Cox模型最大的特点之一是不需要明确指定基准风险函数 h₀(t) 的具体形式。因此，它属于半参数模型。

如果某个变量的回归系数为 β，那么该变量对应的风险比（Hazard Ratio，HR）为：


HR = exp(β)

HR可以直观地反映变量变化对事件风险的影响：

HR > 1：事件发生风险增加；
HR < 1：事件发生风险降低；
HR = 1：变量对风险没有明显影响。

例如某个治疗变量的HR为0.70，可以理解为在其他条件相同的情况下，该组的瞬时事件风险约为对照组的70%。

Cox模型中的删失数据

生存分析与普通回归的一个重要区别是存在删失（Censoring）。

例如研究患者从治疗开始到疾病复发的时间：

患者	时间	是否复发
A	10个月	1
B	15个月	0
C	8个月	1
D	20个月	0

这里通常使用一个事件状态变量表示结果：


1 = 事件发生
0 = 数据被删失

患者B虽然观察了15个月，但在观察结束时没有发生复发，因此我们只能知道其生存时间至少超过15个月，并不知道真正的复发时间。

Cox模型可以利用这类删失数据，而不需要简单地把这些样本直接删除。

Python安装lifelines

使用Python实现Cox比例风险回归，常用的库是 lifelines。

可以通过pip安装：

Bash
pip install lifelines

如果使用Anaconda，也可以执行：

Bash
conda install -c conda-forge lifelines

安装完成后进行导入：

Python
运行
import pandas as pd
from lifelines import CoxPHFitter

建议同时使用较新的Python环境和lifelines版本，以避免不同版本之间API差异带来的问题。

准备Cox回归数据

假设有一份患者生存数据：

Python
运行
import pandas as pd

data = pd.DataFrame({
    "age": [45, 52, 60, 39, 67, 55, 48, 72],
    "treatment": [1, 0, 1, 0, 1, 1, 0, 0],
    "stage": [1, 2, 3, 1, 3, 2, 2, 3],
    "time": [12, 8, 20, 15, 6, 18, 10, 5],
    "event": [1, 1, 0, 1, 1, 0, 1, 1]
})

print(data)

这里：

time：随访时间或生存时间；
event：事件是否发生；
age：年龄；
treatment：治疗组标记；
stage：疾病分期。

对于Cox模型而言，最重要的是明确哪个字段表示持续时间，哪个字段表示事件状态。

使用CoxPHFitter拟合模型

lifelines提供了 CoxPHFitter 类，可以直接进行Cox比例风险回归。

Python
运行
from lifelines import CoxPHFitter

cph = CoxPHFitter()

cph.fit(
    data,
    duration_col="time",
    event_col="event"
)

cph.print_summary()

其中：

Python
运行
duration_col="time"

指定生存时间字段。

Python
运行
event_col="event"

指定事件状态字段。

模型拟合完成后，print_summary()会输出回归系数、标准误、统计量、P值、置信区间以及风险比等重要信息。

查看回归系数

可以直接查看模型的回归系数：

Python
运行
print(cph.params_)

例如：


age          0.052
treatment   -0.410
stage        0.680

这些系数本身是在对数风险尺度上的结果，因此通常还需要结合HR进行解释。

例如：

Python
运行
import numpy as np

hr = np.exp(cph.params_)
print(hr)

如果年龄对应的系数为：


0.052

那么：


exp(0.052) ≈ 1.053

说明年龄每增加一个单位，模型估计的瞬时风险大约增加5.3%，前提是其他变量保持不变。

查看Hazard Ratio风险比

在实际分析中，HR通常比回归系数更加直观。

可以使用：

Python
运行
print(cph.hazard_ratios_)

也可以直接从模型结果中查看：

Python
运行
cph.print_summary()

假设得到：


             coef    exp(coef)
age          0.052     1.053
treatment   -0.410     0.664
stage        0.680     1.974

可以进行如下解释：

age的HR为1.053，年龄增加可能对应更高的事件风险；
treatment的HR为0.664，治疗组风险低于参考状态；
stage的HR为1.974，疾病分期每增加一级，事件风险可能明显升高。

需要注意，HR的解释必须结合变量的编码方式。尤其是分类变量，如果直接使用 0/1 编码，那么HR描述的是从0组变化到1组时风险的相对变化。

P值与统计显著性

Cox回归的结果通常需要结合P值判断变量是否具有统计学意义。

可以通过：

Python
运行
cph.summary

查看完整结果。

例如：

Python
运行
print(cph.summary[
    ["coef", "exp(coef)", "se(coef)", "p"]
])

如果某变量：


p < 0.05

通常可以认为该变量与生存结局之间存在统计学上的显著关联。

但不能简单地将：


p < 0.05

理解为“该变量一定具有实际意义”。

实际研究中还应该结合：

HR大小；
95%置信区间；
样本量；
临床或业务背景；
模型假设；
变量编码方式。

共同判断结果。

查看HR的置信区间

风险比最好不要只报告一个点估计，而应该同时报告置信区间。

可以查看：

Python
运行
print(cph.summary[
    ["exp(coef)", "exp(coef) lower 95%",
     "exp(coef) upper 95%", "p"]
])

例如：


treatment    HR=0.66    95% CI: 0.45-0.97

这意味着模型估计治疗变量对应的风险比为0.66，同时95%置信区间为0.45到0.97。

如果HR的95%置信区间跨过1，则通常意味着在对应显著性水平下，无法认为该变量的风险效应具有统计学显著性。

检验Cox模型的比例风险假设

Cox模型有一个非常重要的前提，即比例风险假设（Proportional Hazards Assumption）。

简单来说，该假设认为两个个体之间的风险比在研究期间应该保持相对稳定。

例如某治疗方案的HR为0.70，比例风险假设要求这种相对风险关系不会随着时间发生系统性的改变。

可以使用lifelines进行检查：

Python
运行
cph.check_assumptions(data, p_value_threshold=0.05)

该方法可以帮助发现可能违反比例风险假设的变量。

如果某个变量明显违反该假设，就不能直接忽略。可以考虑：

增加时间交互项；
对变量进行分层；
对变量进行适当变换；
使用时间变化协变量；
改用其他生存分析模型。

比例风险假设检查是Cox回归分析中非常重要的一步，不能只拟合模型并查看P值后就结束。

对分类变量进行处理

如果数据中存在多个类别的分类变量，不建议简单地把类别编号当成连续变量。

例如疾病类型：


1 = A型
2 = B型
3 = C型

直接放进Cox模型会隐含一种连续关系，即C型相对于B型的影响与B型相对于A型具有固定的线性差异，这通常并不合理。

更常见的方法是使用哑变量。

可以借助pandas：

Python
运行
data = pd.get_dummies(
    data,
    columns=["disease_type"],
    drop_first=True
)

例如原来的：


disease_type
A
B
C

可能转换为：


disease_type_B
disease_type_C

其中A作为参考组。

随后：

Python
运行
cph.fit(
    data,
    duration_col="time",
    event_col="event"
)

模型中的HR就可以解释为B组、C组相对于A组的风险比。

连续变量是否需要标准化

Cox模型并不要求所有连续变量必须标准化。

例如：

Python
运行
age
blood_pressure
cholesterol

可以直接进入模型。

不过，当变量尺度差异非常大，或者变量之间存在较强的数值尺度差异时，可以考虑标准化：

Python
运行
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()

data[["age", "blood_pressure"]] = scaler.fit_transform(
    data[["age", "blood_pressure"]]
)

标准化后的HR解释方式也会发生变化。例如年龄标准化后，HR表示年龄增加一个标准差所对应的风险变化，而不再是“每增加1岁”的风险变化。

因此，标准化并不只是为了让模型运行得更好，也会直接影响结果的解释方式。

处理缺失值

Cox模型建模前需要关注缺失数据。

可以先查看：

Python
运行
print(data.isnull().sum())

如果存在缺失值，需要根据实际业务场景选择处理方法。

简单情况下可以删除缺失样本：

Python
运行
data = data.dropna()

但对于样本量较小的数据集，直接删除可能造成明显的信息损失。

更复杂的研究通常会考虑：

均值或中位数填补；
多重插补；
基于模型的插补；
将缺失作为独立类别。

缺失值处理方式应该在分析报告中明确说明，而不是为了让代码正常运行而随意删除数据。

使用公式直接指定Cox模型

lifelines还支持通过公式选择变量。

例如：

Python
运行
cph = CoxPHFitter()

cph.fit(
    data,
    duration_col="time",
    event_col="event",
    formula="age + treatment + stage"
)

cph.print_summary()

这种方式对于变量较多的数据集更加清晰。

如果需要加入变量之间的交互项，也可以使用公式表达式，例如：

Python
运行
formula = "age + treatment + stage + age:treatment"

不过，加入交互项后模型解释会更加复杂，需要结合具体研究假设进行分析。

Cox模型的预测方法

模型拟合后，可以对新的样本进行风险预测。

例如：

Python
运行
new_data = pd.DataFrame({
    "age": [50, 65],
    "treatment": [1, 0],
    "stage": [2, 3]
})

risk = cph.predict_partial_hazard(new_data)

print(risk)

predict_partial_hazard()得到的是部分风险，也可以理解为相对于基准风险的相对风险指标。

如果希望获得生存概率，可以使用：

Python
运行
survival = cph.predict_survival_function(new_data)

print(survival)

得到的是不同时间点对应的生存概率估计。

还可以绘制生存曲线：

Python
运行
import matplotlib.pyplot as plt

cph.predict_survival_function(new_data).plot()

plt.xlabel("Time")
plt.ylabel("Survival probability")
plt.show()

通过曲线可以直观比较不同协变量组合对应的预测生存概率。

Cox回归与Kaplan-Meier曲线的区别

Cox回归和Kaplan-Meier（KM）分析都属于生存分析，但用途并不完全相同。

Kaplan-Meier更适合回答：

不同组的生存概率随时间如何变化？

例如：

Python
运行
from lifelines import KaplanMeierFitter

kmf = KaplanMeierFitter()

kmf.fit(
    data["time"],
    event_observed=data["event"]
)

kmf.plot_survival_function()

Cox回归则更适合回答：

在同时考虑年龄、治疗方式、疾病分期等变量后，这些因素与事件风险之间是什么关系？

因此，在实际项目中，两种方法往往可以结合使用：先通过KM曲线观察不同组的生存情况，再使用Cox回归控制混杂因素并量化变量对应的风险比。

单因素Cox回归

当变量较多时，可以先分别建立单因素Cox模型。

例如：

Python
运行
variables = ["age", "treatment", "stage"]

for variable in variables:
    cph = CoxPHFitter()
    cph.fit(
        data[["time", "event", variable]],
        duration_col="time",
        event_col="event"
    )

    print(variable)
    print(cph.summary[["coef", "exp(coef)", "p"]])

这样可以分别查看每个变量与生存结局之间的关联。

不过，单因素分析结果不能直接替代多因素分析。一个变量在单因素分析中显著，并不代表加入其他变量后仍然显著。

多因素Cox回归

多因素Cox回归可以同时纳入多个协变量：

Python
运行
cph = CoxPHFitter()

cph.fit(
    data,
    duration_col="time",
    event_col="event"
)

cph.print_summary()

例如同时考虑：


年龄
治疗方式
疾病分期
肿瘤大小
实验室指标

模型可以估计每个变量在控制其他变量后的独立关联。

这也是实际医学统计、生存预测和风险因素研究中最常见的Cox模型使用方式之一。

检查变量之间的共线性

如果两个或多个变量高度相关，可能导致Cox模型参数不稳定。

例如：


身高
体重
BMI

这些指标之间可能存在明显相关性。

可以使用相关系数矩阵进行初步检查：

Python
运行
print(data[["age", "stage"]].corr())

对于更复杂的数据，可以进一步计算VIF等指标。

如果发现严重共线性，可以根据研究目的：

删除高度相关变量；
合并变量；
重新定义指标；
使用降维方法；
采用正则化Cox模型。

不要为了追求更多变量而盲目把所有字段都放进模型。

Cox模型中的正则化

当变量数量较多，或者存在变量选择问题时，可以考虑正则化。

lifelines的 CoxPHFitter支持惩罚项，例如：

Python
运行
cph = CoxPHFitter(
    penalizer=0.1
)

cph.fit(
    data,
    duration_col="time",
    event_col="event"
)

正则化可以在一定程度上降低模型过拟合风险，尤其适用于变量较多而样本量有限的场景。

如果希望使用L1和L2混合惩罚，可以结合 l1_ratio：

Python
运行
cph = CoxPHFitter(
    penalizer=0.1,
    l1_ratio=0.5
)

其中：

l1_ratio=0：更接近L2正则化；
l1_ratio=1：L1正则化；
介于0和1之间：Elastic Net风格的组合。

具体参数需要通过合理的模型评估方法进行选择。

Cox模型常见错误

事件变量编码错误

事件列应该清晰表示事件是否发生。

通常：


1 = 发生事件
0 = 未发生事件或删失

如果把编码方向弄反，会导致整个模型解释出现问题。

生存时间存在负数或无效值

生存时间通常应该是非负数，并且需要符合具体研究设计。

可以检查：

Python
运行
print(data["time"].describe())
print((data["time"] < 0).sum())

把分类变量当作连续变量

例如把：


低风险 = 1
中风险 = 2
高风险 = 3

直接作为连续变量使用，可能产生不合理的线性假设。

需要根据变量本身的统计含义决定是否进行哑变量编码。

忽略比例风险假设

这是Cox回归中特别常见的问题。

即使模型能够成功运行，也不代表模型假设一定成立。

建议拟合模型后执行：

Python
运行
cph.check_assumptions(
    data,
    p_value_threshold=0.05
)

并结合统计结果和实际业务含义进行判断。

只关注P值

Cox分析不应该只看P值。

例如：


HR = 1.02
P < 0.001

虽然统计上可能显著，但实际影响程度可能很小。

反过来：


HR = 0.55
P = 0.08

也不能简单得出“没有任何影响”的结论，还需要考虑样本量、置信区间和统计效能。

一个完整的Cox回归示例

下面给出一个相对完整的Python实现：

Python
运行
import pandas as pd
from lifelines import CoxPHFitter

data = pd.DataFrame({
    "age": [45, 52, 60, 39, 67, 55, 48, 72,
            51, 63, 42, 58],
    "treatment": [1, 0, 1, 0, 1, 1, 0, 0,
                  1, 0, 1, 0],
    "stage": [1, 2, 3, 1, 3, 2, 2, 3,
              1, 3, 2, 2],
    "time": [12, 8, 20, 15, 6, 18, 10, 5,
             22, 7, 16, 11],
    "event": [1, 1, 0, 1, 1, 0, 1, 1,
              0, 1, 0, 1]
})

# 检查缺失值
print(data.isnull().sum())

# 创建Cox模型
cph = CoxPHFitter()

# 拟合模型
cph.fit(
    data,
    duration_col="time",
    event_col="event"
)

# 输出模型结果
cph.print_summary()

# 查看HR
print("Hazard Ratio:")
print(cph.hazard_ratios_)

# 检查比例风险假设
cph.check_assumptions(
    data,
    p_value_threshold=0.05
)

实际项目中，还应该根据数据规模、变量类型以及研究设计进一步完善数据预处理和模型验证。

如何正确解释Cox回归结果

假设模型得到：


Variable      HR      95% CI        P
age           1.04    1.01-1.07     0.008
treatment     0.62    0.43-0.90     0.011
stage         1.85    1.35-2.54     <0.001

可以按照以下思路解释。

对于年龄：

年龄每增加一个单位，事件发生风险约增加4%，其他变量保持不变。

对于治疗方式：

与参考组相比，治疗组的事件发生风险约为参考组的62%。

对于疾病分期：

疾病分期每增加一级，事件发生风险约增加85%。

这里的“风险增加”指的是瞬时风险（hazard），并不等同于“最终发生事件的概率增加85%”。这是解释Cox回归时需要特别注意的一点。

Cox模型结果的可视化

除了直接查看表格，还可以绘制模型的系数或风险比。

lifelines提供了方便的绘图方法：

Python
运行
cph.plot()

import matplotlib.pyplot as plt

plt.show()

也可以绘制HR及其置信区间，从而更直观地观察各变量的风险方向。

对于正式研究报告，通常建议使用森林图（Forest Plot）展示多变量Cox回归结果。森林图能够同时展示变量名称、HR、95%置信区间以及参考线，使读者快速判断风险因素和保护因素。

实际应用中的分析流程

一个较为规范的Python Cox比例风险回归分析流程可以概括为：


原始数据
   ↓
明确生存时间和事件状态
   ↓
检查缺失值和异常值
   ↓
处理分类变量
   ↓
探索变量分布和相关性
   ↓
单因素Cox分析
   ↓
建立多因素Cox模型
   ↓
检查比例风险假设
   ↓
评估共线性和模型稳定性
   ↓
计算HR及95%置信区间
   ↓
模型验证与结果可视化
   ↓
结合实际场景解释结果

需要特别强调的是，Cox回归不是“调用一个函数得到P值”这么简单。数据编码、删失机制、变量选择、比例风险假设和模型验证都会影响最终结论。

总结

Python中的 lifelines 为Cox比例风险回归提供了较完整的实现方案。使用 CoxPHFitter 可以方便地完成模型拟合、HR计算、置信区间分析、预测生存概率以及比例风险假设检验。

最基本的代码结构可以浓缩为：

Python
运行
from lifelines import CoxPHFitter

cph = CoxPHFitter()

cph.fit(
    data,
    duration_col="time",
    event_col="event"
)

cph.print_summary()