生存分析必修：Cox模型残差检验与PH假设验证

残差检验：Cox比例风险模型的假设检验条件

做生存分析的同学应该都听过一句话：跑一个 Cox 模型很简单，但模型的两个前提假设有没有满足，才是真正见功力的时候。

这个系列讲到第三十讲，前面我们把 Cox 模型的基本原理、HR 解读、森林图绘制都过了一遍。今天专门聊一个最容易被忽视、但审稿人和统计分析报告中几乎必问的环节——残差检验。很多人拿到数据，coxph()一行代码跑完，看看系数显著就收工了，结果被审稿人一句 "Have you tested the proportional hazards assumption?" 直接问住。这篇文章就把 Cox 比例风险模型的两大假设、四类残差、R 语言完整检验流程一次讲透。

先说结论：Cox 模型的核心假设是比例风险假设（Proportional Hazards，PH 假设）和对数线性假设（Log-linearity）。检验这两条，靠的就是 Schoenfeld 残差、Martingale 残差、Deviance 残差和 Score 残差。下面我们从"为什么要检验"开始，一步步把检验逻辑和实操代码全部打通。

1. 为什么要关注 Cox 模型的假设条件

1.1 Cox 模型的假设：比例风险到底在说什么

我见过太多人把 Cox 模型当成黑箱，其实它的数学形式非常简洁：

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

这里 h₀(t) 是基线风险函数，它不限制具体分布，所以 Cox 模型叫半参数模型。问题出在 exp(βX) 这一项——它意味着任意两个个体 (A 和 B) 的风险比是：

HR = h(t|X_A) / h(t|X_B) = exp[β(X_A - X_B)]

关键来了：这个比值里没有t。也就是说，无论随访到第 1 天还是第 1000 天，两个个体的风险比始终保持恒定。这就是"比例风险"四个字的由来——各组风险函数在整个时间轴上成固定比例，可以画成一组永远不相交的曲线。

这个假设不满足会怎样？最典型的情况是治疗效果的"时间衰减"。比如某项手术干预，术后 30 天内风险确实显著低于保守治疗，但半年后这种优势消失了——这时 HR 是随时间变化的，如果用单一 HR 来概括，就会出现"平均化错误"，高估或低估真实效应。审稿人问 PH 假设，本质上就是在问：你给出的这个 HR，能不能代表整个随访期？

1.2 第二个假设：对数线性容易被人忽略

第二个常被忽略的假设是对数线性假设：连续型自变量每增加一个单位，对数风险（log hazard）呈线性变化。换句话讲，X 每增加 1 个单位，HR 恒定为 exp(β)。

这个假设不满足的后果同样严重。举个例子，年龄与疾病风险通常不是纯粹的线性关系——年轻段风险平缓，老年段风险陡增。如果你强行把年龄作为连续变量塞进模型，得到的结果是"平均效应"，可能既低估了老年人的风险，又高估了中年人的风险。

很多人以为只有线性回归才需要检验线性假设，其实 Cox 模型作为广义线性模型的一个变体，同样默认对数风险与协变量之间是线性关系。这就是为什么在 Cox 模型里，我们也需要用残差去"审问"这个设定。

1.3 假设不满足时的后果：从统计功效到解释偏差

如果你跳过了假设检验，后果远不止"被审稿人刁难"这么简单：

• HR 变成一个无意义的加权平均：真实风险比随时间变化时，单一 HR 不能描述任何特定时间点的效应。
• 检验功效下降：违反 PH 假设时，Cox 模型对组间差异的检验功效会明显降低，本来该显著的差异可能无法检测出来。
• 预测偏差：基于错误假设建立的模型用于个体风险预测时，偏差会随随访时间增大。
• 参数解释困难：分类协变量（比如治疗组与对照组）违反了 PH 假设时，分组因素与时间之间存在交互作用，β 本身就失去了"总体效应"的含义。

所以，不管是发表论文还是做临床预测模型，假设检验这一关必须过。

2. 四类残差：谁负责检验什么

既然要检验假设，就得有工具。Cox 模型的残差和普通线性回归的残差（观测值减预测值）不是一个概念，它没那么直观。好消息是，R 的survival包帮我们把四种残差都封装好了，但要选对类型、看懂结果，你得先明白它们各自"管什么"。

2.1 Schoenfeld 残差：PH 假设的专用探针

Schoenfeld 残差的核心逻辑是：对每一个发生事件的时间点，计算"在这个时间点发生事件的个体的协变量观测值"与"此时所有仍处于风险集的个体的协变量期望值"之差。

听起来复杂，但你只需要理解一个直觉：如果 PH 假设成立，那么 Schoenfeld 残差与时间没有系统性关联，残差点应围绕零线随机波动。如果残差随时间呈现上升、下降或先升后降的趋势，说明该变量的效应不是固定的——比例风险假设可能被违反。

更妙的是，对 Schoenfeld 残差做缩放处理后，其平滑曲线可以直接解释为"该变量的系数 β 随时间的变化轨迹"。也就是说，你能直观看到：性别效应在随访早期很大，到后期是否趋于消失。

2.2 Martingale 残差：用于检查函数形式

Martingale 残差的定义和我们熟悉的残差差异很大。它的取值范围是 (-∞, 1]，可以理解为"删失指示变量"与"模型预测累积风险"之间的差值。对删失个体，残差通常是负值；对发生事件的个体，正值说明实际事件比预测来得更早。

它的主要用途是检查协变量的函数形式。比如你把年龄按连续变量放入模型，如果 Martingale 残差随年龄变化的平滑曲线明显偏离零线或者呈非线性形状，说明年龄的线性设定不合理，可能需要考虑分段、二次项或样条变换。这个思路和线性回归中"残差对自变量作图"异曲同工。

2.3 Deviance 残差：找异常值和强影响点

Deviance 残差是对 Martingale 残差做标准化变换得到的，好处是让分布更接近正态，从而更容易识别异常值。通常我们约定：当 Deviance 残差的绝对值过大（比如大于 3 或 4）时，该样本可能对模型结果产生很大影响，值得拿出来单独检查。

实践中，Deviance 残差常用来排查数据录入错误、极端预后个体等"可疑样本"。我会在后面演示怎么用阈值筛选这些点。

2.4 Score 残差（dfbeta）：评估每个观测对系数的杠杆

Score 残差还有一个更好记的名字：dfbeta 残差。它直接衡量"如果删除第 i 个观测，模型中的各个回归系数 β 会变化多少"。变化量大，说明这个观测是强影响点，也许正在"拽着"你的 HR 往某个方向偏。

四种残差各有分工，下一个表格可以直接抄走：

3. R 语言实战：PH 假设的完整检验流程

3.1 环境准备与数据加载

我用 R 自带数据集lung做演示。这个数据来自北美肺癌研究组，记录了 228 名晚期肺癌患者的生存时间、生存状态和若干临床特征。虽然数据年代久远，但胜在拿过来就能用，很适合教学演示。

# 加载所需的包 library(survival) # Cox 模型及残差检验核心包 library(survminer) # 可视化辅助包 library(dplyr) # 数据操作 library(ggplot2) # 高级绘图 # 载入内置数据集并查看结构 data("lung") str(lung)

在lung数据集中，time表示生存天数，status中 1 表示删失、2 表示死亡，age为年龄，sex为性别（1=男，2=女），ph.ecog为 ECOG 体力评分（0-5，分数越高体能越差）。演示前先做一个简单处理：把status转成 0/1 格式（0=删失，1=事件），并给sex加上标签，方便看结果。

lung <- lung %>% mutate( status_bin = ifelse(status == 2, 1, 0), sex_label = factor(sex, levels = c(1, 2), labels = c("Male", "Female")) )

3.2 建立基础 Cox 模型

我们先建立一个包含年龄、性别和体能评分的 Cox 模型，这也是论文中最常见的组合：

# 建立 Cox 比例风险模型 fit_cox <- coxph(Surv(time, status_bin) ~ age + sex_label + ph.ecog, data = lung) # 查看模型摘要 summary(fit_cox)

输出的exp(coef)一列就是 HR。比如ph.ecog的 HR 通常大于 1，表示 ECOG 评分越高死亡风险越大。性别变量的 HR 小于 1，表示女性（标签为 Female 的组）相对男性风险更低。输出中还有一个 Likelihood ratio test 的 p 值，用于整体模型检验。

3.3 用cox.zph()一行代码做 PH 假设检验

下面进入正题。R 中检验 PH 假设的标准做法是cox.zph()函数，它利用的正是缩放 Schoenfeld 残差：

# 检验 PH 假设 ph_test <- cox.zph(fit_cox) # 输出检验结果 print(ph_test)

输出结果大致长这样：

这里的原假设是"该变量的风险比不随时间变化（即满足 PH 假设）"。如果 p 值小于 0.05，就需要警惕：该变量可能不满足比例风险假设。GLOBAL 则是整体检验，对所有变量做联合检验。

用上面这组结果来说，三个变量的 p 值都大于 0.05，GLOBAL 的 p 值也远大于 0.05，说明在这个数据上，PH 假设是成立的。不过我要强调一点：p 值只是参考，图形判断更加重要。

3.4 画图：直接从图形看风险是否成比例

光看 p 值容易踩坑，尤其是样本量大的时候，微小偏离也会让 p 值很显著。所以我强烈建议每次做cox.zph()都出图看一眼。

# 绘制 Schoenfeld 残差图 plot(ph_test)

也可以用ggcoxzph()让图更美观：

library(survminer) ggcoxzph(ph_test)

图形中每一行对应一个协变量。图中的实线是缩放 Schoenfeld 残差的平滑拟合线，虚线是置信带。如果中间的实线大致水平，并且拟合线的变化幅度始终在置信带内部，说明该变量满足 PH 假设。如果平滑曲线呈现明显的上升/下降趋势或者穿过置信带，说明该变量的效应随时间发生系统性变化。

如果某个变量不满足 PH 假设，你还能从图中看出变化方向：曲线上升意味着效应随时间增强，下降则意味着效应随时间减弱（比如治疗获益逐渐消失）。

3.5 按时间分组检验：处理 p 值不显著但图形可疑的情况

有一个场景很常见：p 值大于 0.05，但图形上看起来曲线在早期有波动。这时候可以用cox.zph()的transform参数做补充分析。默认情况下，cox.zph()使用 Kaplan-Meier 变换后的时间尺度，也可以换成其他变换：

# 用对数时间尺度重新检验 cox.zph(fit_cox, transform = "log") # 用原始时间尺度检验 cox.zph(fit_cox, transform = "identity")

为什么要换时间尺度？因为不同变换对时间轴不同区间的敏感性不同。默认的 KM 变换对删除/事件分布比较均匀的时间段更敏感，而对数变换则更关注早期与晚期的对比。实际操作中，我会把默认结果和图形结合起来判断，而不是机械地看某一个 p 值。

3.6 对分层变量和分类变量的 PH 检验

如果模型中包含多分类变量（比如血型 A/B/AB/O），应该先将它转换为因子再纳入模型，此时cox.zph()会自动对该因子做整体的检验：

# 因子变量的 PH 检验 lung$ph.ecog_factor <- factor(lung$ph.ecog) fit_factor <- coxph(Surv(time, status_bin) ~ age + sex_label + ph.ecog_factor, data = lung) cox.zph(fit_factor)

注意，这种情况下 Schoenfeld 残差图会有多个面板，每个面板对应因子变量中的一个哑变量。解读方式和连续变量一致，只要有一个哑变量的趋势明显，就说明这个分类变量的比例风险假设可能存在问题。

4. 对数线性假设与异常点诊断实操

4.1 利用 Martingale 残差检查连续变量的函数形式

PH 假设通过后，接下来要检查对数线性假设。这里用 Martingale 残差：

# 拟合一个只含截距的模型（或目标模型），提取 Martingale 残差 res_martingale <- residuals(fit_cox, type = "martingale") # 把残差和年龄画在一起，看是否存在非线性趋势 plot_data <- data.frame(age = lung$age, martingale = res_martingale) ggplot(plot_data, aes(x = age, y = martingale)) + geom_point(alpha = 0.5) + geom_smooth(method = "loess", se = TRUE, color = "red") + geom_hline(yintercept = 0, linetype = "dashed") + labs(title = "Martingale Residuals vs. Age", x = "Age", y = "Martingale Residual")

如果红色平滑线基本在零线附近水平波动，说明年龄的线性设定可以接受。如果曲线呈明显的倒 U 型或 S 型，则要考虑年龄的变换形式。

另一种更正式的检查方式是 Box-Tidwell 检验思路：把连续变量的交互项纳入模型，看交互项是否显著。比如检验年龄是否满足对数线性，可以构造age * log(time)交互项，但这在模型解释上比较复杂，实践中我更推荐直接看 Martingale 残差图 + 用样条函数验证。

4.2 用 Deviance 残差识别异常样本

检查完函数形式，顺手用 Deviance 残差扫一遍异常值：

# 计算 Deviance 残差 res_deviance <- residuals(fit_cox, type = "deviance") # 看看哪些样本的残差异常大 outliers <- which(abs(res_deviance) > 3) print(outliers) # 查看这些样本的原始记录 lung[outliers, c("time", "status_bin", "age", "sex", "ph.ecog")]

Deviance 残差的阈值没有绝对标准，但绝对值大于 3 是一个常用的参考线。找到这些异常点之后，不要急着删除，而是要认真核查：是数据录入错误？是极端预后的真实病例？还是某些协变量组合本身就罕见？如果属于最后一种，建议在论文的敏感性分析部分说明"剔除异常点后结果方向一致"，以证明结论的稳健性。

4.3 dfbeta 残差评估强影响点

dfbeta 可以直接量化每个观测对模型的杠杆作用：

# 计算 dfbeta 残差 dfbeta_res <- residuals(fit_cox, type = "dfbeta") # 查看前几行 head(dfbeta_res) # 用索引图画每个变量的 dfbeta par(mfrow = c(2, 2)) for (i in 1:ncol(dfbeta_res)) { plot(dfbeta_res[, i], type = "h", main = colnames(dfbeta_res)[i], ylab = "dfbeta", xlab = "Observation Index") abline(h = 0, lty = 2) }

如果某条竖线明显比其它观测高出很多，说明这个观测对相应回归系数的估计影响很大。此时可以做一个快速敏感性分析：删掉这个观测重新拟合模型，看 HR 和 p 值是否发生质的改变。如果结果差别很大，说明结论对这个样本过于敏感，需要在讨论中如实说明。

5. PH 假设不满足时的三大应对策略

如果检验结果确实显示某个变量违反了 PH 假设，不要慌，这在实际数据中非常常见。我有三个常用方案，按实施难度从低到高介绍。

5.1 方案一：分层处理

最简单的办法是对违反 PH 假设的变量做分层，也就是允许不同层拥有不同的基线风险函数：

# 直接用 strata() 包裹该变量 fit_strata <- coxph(Surv(time, status_bin) ~ age + sex_label + strata(ph.ecog), data = lung) summary(fit_strata)

分层后，ph.ecog不再估计回归系数，但每一层都有自己的 h₀(t)，相当于各层之间互不干扰，风险比自然不再受 PH 假设约束。代价是：你无法直接回答"ph.ecog 每增加一级风险增加多少"这个问题。如果分层变量是你最关心的暴露因素，分层策略就不太适合了。

5.2 方案二：时间交互项

保持变量在模型中，但允许其效应随时间变化。具体做法是在coxph()中通过tt()指定时间变换函数：

# 假设 ph.ecog 违反 PH 假设，加入其对数时间交互项 fit_tt <- coxph( Surv(time, status_bin) ~ age + sex_label + ph.ecog + tt(ph.ecog), data = lung, tt = function(x, t, ...) x * log(t) ) summary(fit_tt)

这时ph.ecog的 HR = exp(β₁ + β₂ * log(t))，是时间的函数。你可以在文中报告：随访早期 ph.ecog 每增加一级的 HR 是多少，随访到某个时间点后又变成多少。这种方案解释起来稍微复杂，但在临床研究中特别受欢迎，因为它揭示的是"效应随时间如何演变"。

常见的tt变换函数有x * log(t)、x * t、x * sqrt(t)，具体选哪个可以用AIC()比较不同模型的拟合优度。注意，tt()中的t必须是模型内部认可的时间变量，不可自己定义，这点很多人容易写错。

5.3 方案三：扩展为时变系数模型或拆分时间区间

如果时间交互项无法充分捕捉时变效应，可以考虑更灵活的时变系数模型。R 中有几个方向：

• survival包的重叠版本支持tt()更复杂的函数。
• timereg包的timecox()可以估计非参数时变系数，画出的系数曲线非常直观。
• 也可以把时间轴按分位数切分成几段，分段拟合 Cox 模型，如果各段 HR 接近，说明偏差不大：

# 用 survSplit 按时间切分数据 lung_split <- survSplit(Surv(time, status_bin) ~ ., data = lung, cut = c(180, 365), episode = "tgroup") # 拟合分段模型，允许 ph.ecog 在不同时间段有不同效应 fit_split <- coxph(Surv(tstart, time, status_bin) ~ age + sex_label + ph.ecog + ph.ecog:tgroup, data = lung_split) summary(fit_split)

这里的截断点（180、365 天）需要结合临床意义选择，而不是随意取。分段模型的优点是解释直观，缺点是切分点的选择带有主观性，建议做敏感性分析。

6. 常见问题与避坑实录

6.1 常见问题速查表

把实际分析中最常踩的坑整理成表，排查问题时可以直接对号入座：

6.2 不要只盯着 p 值

这是我最想强调的一点。cox.zph()的 p 值受样本量影响极大。样本量达到几千甚至上万时，哪怕是微小的风险比漂移也会得到 p < 0.05；反过来，样本量只有几十例时，即使风险比在时间轴上波动很大，检验也可能不显著。

我自己的经验是：p 值 + 图形 + HR 实际变化幅度三者一起看。比如 HR 从 1.5 缓慢变到 1.6，即使 p 值显著，临床意义可能有限；如果 HR 从 2.0 直线降到 1.0，p 值不显著也要警惕，很可能是功效不足。计算不同时间段的 HR 可以用时间交互项模型输出预测值来实现。

6.3 多重比较的顾虑

模型里协变量越多，cox.zph()检验的变量也越多，出现假阳性的概率随之上升。当变量很多时，可以考虑用 Bonferroni 校正，比如 α 从 0.05 调整为 0.05/k（k 为检验的变量个数）。更实际的做法是优先关注主要暴露因素和临床上已知会随时间变化的因素，把 PH 检验当作有针对性的诊断，而不是机械地做一遍筛查。

6.4 删失比例极高时怎么办

如果数据删失比例超过 80%，事件数很少，残差检验的功效会非常有限。此时 Schoenfeld 残差图上的点可能非常稀疏，平滑曲线也不可靠。这种情况下，我建议缩减模型，只保留核心变量；或者在方法学部分坦诚说明"由于事件数有限，PH 假设检验功效不足，结果应谨慎解读"。

6.5 一个完整的分析流程模板

最后给大家一个可以直接套用的完整流程：

1. 检查生存数据格式（time + status 0/1）。
2. 用coxph()建立初步模型。
3. 用cox.zph()对模型做整体 PH 检验，并出图。
4. 对连续变量做 Martingale 残差图，检查线性设定。
5. 用 Deviance 残差和 dfbeta 检查异常点与强影响点。
6. 如果 PH 检验通过，报告模型参数；如果不通过，依次尝试分层、时间交互、分段模型，并用 AIC 或对数似然比比较模型。
7. 对最终模型重新做cox.zph()，确认问题已缓解。

# 最后检查一下最优模型是否满足 PH 假设 cox.zph(fit_final)

我在实际项目中见过太多"一把梭"式的 Cox 分析——跑完模型直接写报告，结果被评审或上级追问"PH 假设验证了吗"的时候才开始补救。提前做好残差检验，不仅能让你的统计结果经得起推敲，更重要的是能帮你发现数据中隐藏的时变效应和异常样本，这些信息本身往往比一个简单的 HR 更有临床价值。把今天这套流程跑熟，生存分析这块基本功就算真正过关了。