Cox 回归：单因素与多因素

跟完这篇，您能把一份带随访时间、结局和多个协变量的数据，先用单因素 Cox 把候选变量筛一遍，再放进多因素模型算出校正后的风险比（HR）与置信区间，检验比例风险假设，画一张森林图，并把结果写成一段审稿人挑不出方法硬伤的结论。

Cox 回归在解决什么问题

KM 曲线加 log-rank 能告诉您两组生存有没有差，但它只给一个 p 值，答不了两个关键问题：差多少，以及排除了年龄、分期这些混杂之后还差不差。Cox 比例风险回归就是来补这两块的。

它的输出核心是风险比（hazard ratio, HR）：

HR > 1：该因素让事件（如死亡）发生得更快，是危险因素。
HR < 1：发生得更慢，是保护因素。
HR = 1：对风险没有影响。

举例，某基因高表达对比低表达 HR = 2.0，意思是「在任一时刻，高表达组正在发生事件的瞬时风险，是低表达组的 2 倍」。

有三点从一开始就得摆正，免得后面写错：

HR 是「风险速率之比」，不是「多少人会死」，也不是相对风险 RR 或比值比 OR。 它说的是发生的快慢，不直接等于绝对生存概率。
HR 是整个随访期的一个「平均」效应。 这个平均能不能代表全程，取决于下面要讲的比例风险假设。
Cox 半参数、不设定基线生存的具体形状，这正是它比参数模型更常用、更稳的原因。

数据要长什么样

和生存分析一样，最少三类列：随访时间、事件结局、以及您想纳入模型的协变量。

列名	含义	取值
`time`	从统一起点到事件、或到最后一次随访的时长	正数（全表统一单位：月 / 天 / 年）
`event`	随访结束时事件到底发生了没有	`1` = 发生，`0` = 删失
`age`	连续型协变量示例	数值
`sex` / `stage`	分类型协变量示例	因子（factor）

head(df)
#>   sample  time  event  age  sex   stage  gene
#> 1   S001   12.4     1   67    M     III  high
#> 2   S002   34.0     0   58    F      II  high     # 随访 34 个月仍存活，删失
#> 3   S003    5.1     1   72    M     III   low

分类变量务必先转成因子，并显式定好参照水平（reference level）——所有 HR 都是相对参照水平来解释的，参照选错，整张表都要重读一遍：

df$stage <- factor(df$stage, levels = c("I", "II", "III"))   # 以 I 期为参照
df$gene  <- factor(df$gene,  levels = c("low", "high"))       # 以低表达为参照

单因素 Cox：先把候选变量筛一遍

多因素模型不宜一股脑塞进所有变量。常规做法是先逐个变量单独跑 Cox，看哪些和结局有关联，再把有苗头的带进多因素模型。

单个变量的写法很短，summary() 里就有 HR、置信区间和 p 值：

library(survival)

fit <- coxph(Surv(time, event) ~ gene, data = df)
summary(fit)
#>            coef  exp(coef)  se(coef)     z   Pr(>|z|)
#> genehigh  0.693     2.000     0.180  3.85    0.00012
#>
#>          exp(coef)  exp(-coef)  lower .95  upper .95
#> genehigh     2.000       0.500      1.406      2.845

读这张输出：exp(coef) 就是 HR（此处 2.0），lower .95 / upper .95 是它的 95% 置信区间（1.41–2.85），Pr(>|z|) 是 p 值。置信区间不跨过 1，与 p<0.05 是一回事，都提示关联显著。

变量多时，套个循环逐个跑，把 HR、CI、p 汇成一张表：

vars <- c("age", "sex", "stage", "gene")

uni_tab <- lapply(vars, function(v) {
  f  <- as.formula(paste("Surv(time, event) ~", v))
  m  <- coxph(f, data = df)
  s  <- summary(m)
  data.frame(
    变量 = rownames(s$coefficients),
    HR   = round(s$coefficients[, "exp(coef)"], 2),
    下限 = round(s$conf.int[, "lower .95"], 2),
    上限 = round(s$conf.int[, "upper .95"], 2),
    p    = signif(s$coefficients[, "Pr(>|z|)"], 2)
  )
})
uni_tab <- do.call(rbind, uni_tab)

筛选阈值要说清一件事：进入多因素的门槛，常按惯例定在单因素 p<0.05，很多研究会放宽到 p<0.1，好让边缘但可能重要的变量有机会被校正。这是常用惯例、按研究和领域调，不是铁律。另外两点别忘：

临床上公认重要的混杂（如年龄、分期），即便单因素不显著，通常也应保留进多因素模型，因为它们的作用可能被其他变量掩盖。
别用纯自动逐步回归（stepwise）挑变量当唯一依据：它容易过拟合、结果不稳、还难复现。变量取舍最好有临床或前期证据支撑。

多因素 Cox：一起放进模型，得到校正后的 HR

把筛出来的变量用 + 一起写进公式，就是多因素模型。它给出的每个 HR，都是在同时固定其他变量的前提下，该变量单独的效应——这正是「校正混杂」的含义。

multi <- coxph(Surv(time, event) ~ age + stage + gene, data = df)
summary(multi)
#>            coef  exp(coef)  se(coef)     z   Pr(>|z|)
#> age       0.021     1.021     0.009  2.33    0.0198
#> stageII   0.480     1.616     0.240  2.00    0.0455
#> stageIII  0.990     2.691     0.250  3.96   7.5e-05
#> genehigh  0.560     1.751     0.185  3.03    0.0024
#>
#> Concordance = 0.71 (se = 0.02)

怎么读、怎么落到解释：

genehigh 的 HR 从单因素的 2.0 降到多因素的 1.75，说明原来那份差异里，有一部分其实是年龄、分期带来的；校正之后基因本身仍有独立作用。这种「校正前后 HR 变化」本身就是要报告、要讨论的信息。
连续变量看「每增加一个单位」的效应。 age 的 HR = 1.021，是每长 1 岁风险增加约 2%。若想按每 10 岁报告，可先把变量除以 10 再建模，解释更直观——换算尺度只影响表述，不改变模型本身。
分类变量逐水平对参照解释。 stageIII 的 HR = 2.69，指 III 期对比 I 期（参照）的风险；II 期对比 I 期为 1.62。
Concordance（C-index） 衡量模型区分预后好坏的能力，0.5 等于随机、越接近 1 越好，可作为整体表现的一个参考指标。

一条常被审稿人追问的经验法则：每纳入一个变量，至少要有约 10 个事件（events per variable, EPV ≥ 10）。事件数相对变量太少，模型会不稳、HR 偏离、置信区间虚宽。这是广为引用的经验惯例、非硬性规定，但纳入变量前值得先数一数 event == 1 有多少个。

比例风险假设是什么，怎么检验

Cox 之所以能用一个 HR 概括全程，靠的是一条前提——比例风险（proportional hazards, PH）假设：两组的风险之比在整个随访期里大致恒定，不随时间变。 一旦这个假设不成立，那个单一 HR 就是把「早期差很多、晚期不差」之类的时变效应硬压成了一个平均值，会误导读者。

标准检验用 Schoenfeld 残差，函数是 cox.zph()：

zph <- cox.zph(multi)
zph
#>          chisq  df      p
#> age       0.72   1  0.40
#> stage     5.90   1  0.015    # p 小，提示该变量可能违反 PH
#> gene      0.31   1  0.58
#> GLOBAL    6.80   3  0.079

plot(zph)      # 看每个变量的 Schoenfeld 残差随时间是否大致水平（水平即满足）

判读：某个变量或 GLOBAL 的 p<0.05，提示该项可能违反 PH 假设（p<0.05 仍是常用惯例、按数据判断，不是铁律）；plot() 里若残差曲线随时间明显上翘或下滑、而非大致水平，也是同样的信号。除此之外，KM 曲线交叉，往往也预示 PH 不成立。

假设被违反时，几条标准的处理路子（按情形选，不必全上）：

做法	适用	写法要点
分层（stratify）	该变量是混杂、但您不关心它的 HR	`+ strata(stage)`，让各层有各自的基线风险
加时间交互项	您关心该变量、且其效应确随时间变	在模型中引入变量与时间的交互（需按 `cox.zph` 结果与时间函数设定）
分时段报告 HR	效应有明显阶段性	分早期 / 晚期各报一个 HR，并说明为何分段

最要紧的是：如实报告您检验了 PH 假设、结果如何、违反时怎么处理。藏着不报，比违反本身更容易被挑。

画森林图

森林图是把多因素结果一图看尽的标准图：每个变量一行，点是 HR，横线是 95% 置信区间，竖直参考线在 HR = 1——线越靠右、且整条不压到 1，危险作用越确凿；压过 1 则不显著。

最省事的是 survminer::ggforest，直接吃 coxph 对象：

library(survminer)

ggforest(
  multi, data = df,
  main = "多因素 Cox 回归",
  fontsize = 0.9
)

想要更强的排版控制（自定义分组、合并单因素与多因素两列 HR），可用 forestplot 包手动搭。读图和写作时记住：看的是「线有没有压到竖线 1」，不是单看点估计离 1 多远——置信区间很宽、只是碰巧点估计偏大的变量，结论要留余地。

结果怎么写进论文

图和表都齐了，落到笔头。一段规范的 Cox 结果，通常把这几件事说全：

两张 HR 一起给。 论文里常见「Table 2」把单因素、多因素的 HR、95% CI、p 并排列出，让读者看清校正前后怎么变。
交代参照水平和变量尺度。 分类变量以哪一水平为参照、连续变量按每几个单位报告，都要写明，否则 HR 无法解释。
写清 PH 假设检验过。 一句「采用 Schoenfeld 残差检验比例风险假设」再加结果，是方法完整性的标配。
一句话范式（可直接改用）：

多因素 Cox 回归显示，在校正年龄与临床分期后，该基因高表达组的死亡风险仍显著高于低表达组（HR = 1.75，95% CI 1.22–2.52，p = 0.002）。

（句中数字仅为示例格式，请换成您自己的真实结果。）

同时把方法边界一并交代，别等审稿人来问：

HR 是时间平均效应，依赖 PH 假设成立。 若假设被违反且未处理，单个 HR 会误导。
只校正了模型里放进去的混杂，没测量到的混杂仍可能有偏；变量间高度共线时，各自 HR 也会不稳。
有竞争风险时（关心的结局会被其他事件「抢走」，如研究本病死亡却有人先死于他因），标准 Cox 会有偏，需竞争风险模型（如 Fine-Gray），本篇不覆盖。
单中心、小样本、EPV 偏低都会让结果偏乐观，如实写进局限。

在平台上跑一遍

以上从单因素筛变量、多因素建模、PH 假设检验到森林图，这整条链在百沐一下都不用写代码：上传带 time、event 和各协变量的数据表，用一句话说清目的——「对这份数据做 Cox 回归，先单因素筛一遍，把 p<0.1 的变量放进多因素模型，以低表达和 I 期为参照，检验比例风险假设，并画森林图」——平台会自动选方法、算 HR 与置信区间、检验 PH、出图，并给出可复现的 R 代码。想改参照水平、筛选阈值或变量尺度，直接在对话里说明即可。上手见智能统计使用指南，或直接开始：https://app.baimuyixia.com 。