R中survival包coxph分层变量做模型诊断报错问题排查

架构老林

2026-8-29

问题描述

数据集包含3组仅暴露于单一媒介的受试对象，初始构建的分层Cox比例风险模型代码如下：

# 暴露变量treatment已预先转换为因子类型
fullModel <- coxph(Surv(time, status) ~ strata(treatment), data = d)

执行比例风险假设诊断（cox.zph()）时触发报错：

test.assump <- cox.zph(fullModel)
Error in cox.zph(fullModel) : 
there are no score residuals for a Null model

移除公式里的strata()包裹后，诊断可正常运行，输出结果为：

chisq df    p
treatment  1.29  2 0.52
GLOBAL     1.29  2 0.52

可复现报错的示例代码：

data <- list(time=c(4,3,1,1,2,2,3,2,4,1,3,4), 
             status=c(1,1,1,0,1,1,0,1,1,0,0,1),  
             treatment=c(0,0,0,0,1,1,1,1,2,2,2,2))

cox.test <- coxph(Surv(time, status) ~ strata(treatment), data = data)

test.coxas <- cox.zph(cox.test)
ggcoxzph(test.coxas)

ggcoxdiagnostics(test.coxas, type = "schoenfeld",
                 linear.predictions = F)

核心疑问：是否可以先移除strata()参数完成模型诊断，之后再加回该参数，用ggsurvplot()绘制不同暴露组的生存曲线？上述操作存在什么错误？

解答

报错的核心原因非常直接：你写的~ strata(treatment)本质是没有任何协变量的空模型。strata()的作用是通知模型对不同分层单独拟合基线风险，不会把括号里的变量作为协变量估计回归系数。模型里连一个待估计的协变量效应都没有，自然没法计算得分残差、做比例风险假设检验，cox.zph()抛出空模型报错完全符合预期。

你设想的「先删掉strata()做诊断，再加回strata()画图」的操作逻辑完全错误，这是两个设定完全不同的模型，结果不能互用：

去掉strata()的~ treatment模型，前提假设是3个暴露组的基线风险完全一致，仅存在成比例的风险比差异。此时cox.zph()输出的检验结果，对应的是「treatment作为协变量时，其效应是否满足比例风险假设」，这个结果只适配把treatment作为协变量的非分层模型，和分层模型没有任何关系。
带strata(treatment)的分层模型，本身就允许3个组的基线风险完全不同，模型根本不会估计treatment的效应值——你都没估计这个变量的效应，自然不存在「这个变量的效应是否违反比例风险假设」的说法，也就没法对它做PH检验。

如果你的研究目的是比较3个暴露组的生存差异、估计暴露的效应值，从一开始就不应该用strata()包裹treatment：直接把treatment作为协变量放入模型即可，此时做cox.zph()诊断是完全正确的，后续用ggsurvplot()绘制分组生存曲线也直接用这个非分层模型就行，不需要额外加分层项。
只有当你需要控制的混杂变量不满足PH假设、且你不需要估计这个混杂变量的效应时，才应该把这类变量放到strata()里。你关心的核心暴露变量treatment绝对不能直接作为分层项，否则模型根本不会输出暴露的效应估计值，完全失去组间比较的意义。
如果你的模型里确实纳入了其他需要分层控制的、违反PH假设的变量，做cox.zph()诊断时只需要检验模型里实际纳入的协变量项即可，分层项本身不需要、也无法通过cox.zph()完成比例风险检验。

内容的提问来源于stack exchange，提问作者asellus