生存分析 Step 1：主 Cox 分析
survival-analysis-step1-zh.Rmd这一步要解决什么问题
这篇教程是 leo.ukb
生存分析部分的第一步。正式进入亚组分析、交互效应或
中介分析之前，我们首先要回答最基础的问题：
这个 cohort 里，exposure 和 incident outcome 之间有没有主关联？
这就是主 Cox 分析的任务。
在真实研究里，把流程拆成几个阶段会更清楚：
- 先建立 exposure-outcome 的主关联
- 再问这个关联会不会在不同人群里不一样
- 再问这个关联有没有可能部分通过某条 mediator 路径传递
所以 Step 1 是后面所有分析的地基。
先从研究问题和变量角色开始
下面继续用 survival::lung 做教学示例，但我们先按真实
cohort 分析的思路 定义好每列变量的角色。
lung_df <- dplyr::transmute(
survival::lung,
outcome = as.integer(status == 2),
outcome_censor = time / 365.25,
age = age,
sex = factor(sex, levels = c(1, 2), labels = c("Male", "Female")),
ecog_group = factor(
ph.ecog,
levels = 0:3,
labels = c("ECOG0", "ECOG1", "ECOG2", "ECOG3")
)
)先把变量角色说清楚，而不是直接开始跑模型。
| 变量 | 类型 | Step 1 角色 |
|---|---|---|
outcome |
二分类事件指示变量 |
incident outcome 是否发生 |
outcome_censor |
以年为单位的随访时间 |
incident outcome 随访时间 |
age |
连续 exposure | 第一个示例里的主 exposure |
sex |
二分类协变量 | 基础人口学调整变量 |
ecog_group |
ECOG 体能状态（4 组） | ECOG 体能状态调整变量 |
建模前先检查分析数据
主 Cox 分析要先让人信服，最简单的办法就是在建模前先把 cohort 结构、事件负担 和随访情况看清楚。
dplyr::glimpse(lung_df)
#> Rows: 228
#> Columns: 5
#> $ outcome <int> 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ outcome_censor <dbl> 0.8377823, 1.2457221, 2.7652293, 0.5749487, 2.4175222, …
#> $ age <dbl> 74, 68, 56, 57, 60, 74, 68, 71, 53, 61, 57, 68, 68, 60,…
#> $ sex <fct> Male, Male, Male, Male, Male, Male, Female, Female, Mal…
#> $ ecog_group <fct> ECOG1, ECOG0, ECOG0, ECOG1, ECOG0, ECOG1, ECOG2, ECOG2,…| outcome | outcome_censor | age | sex | ecog_group |
|---|---|---|---|---|
| 1 | 0.8377823 | 74 | Male | ECOG1 |
| 1 | 1.2457221 | 68 | Male | ECOG0 |
| 0 | 2.7652293 | 56 | Male | ECOG0 |
| 1 | 0.5749487 | 57 | Male | ECOG1 |
| 1 | 2.4175222 | 60 | Male | ECOG0 |
| 0 | 2.7980835 | 74 | Male | ECOG1 |
event_summary <- data.frame(
指标 = c("Participants", "Events", "Event rate"),
数值 = c(
nrow(lung_df),
sum(lung_df$outcome),
sprintf("%.1f%%", 100 * mean(lung_df$outcome))
)
)
knitr::kable(
event_summary,
caption = "总体事件负担"
)| 指标 | 数值 |
|---|---|
| Participants | 228 |
| Events | 165 |
| Event rate | 72.4% |
| sex | outcome | n |
|---|---|---|
| Male | 0 | 26 |
| Male | 1 | 112 |
| Female | 0 | 37 |
| Female | 1 | 53 |
followup_summary <- summary(lung_df$outcome_censor)
knitr::kable(
data.frame(
统计量 = names(followup_summary),
数值 = as.numeric(followup_summary)
),
caption = "以年为单位的随访时间概况"
)| 统计量 | 数值 |
|---|---|
| Min. | 0.0136893 |
| 1st Qu. | 0.4565366 |
| Median | 0.6995209 |
| Mean | 0.8356809 |
| 3rd Qu. | 1.0855578 |
| Max. | 2.7980835 |
这些检查能先回答几个最基础的问题：
- 事件率大概是多少
- follow-up 的量级是否合理
- 不同协变量层次里有没有明显过稀的单元
建模前先决定 exposure 的尺度
真正开始拟合之前，先决定 exposure 是按连续变量解释，还是按分类变量解释。 这会直接决定 hazard ratio 的含义。
scale_table <- data.frame(
Exposure_设定 = c("连续 exposure", "分类 exposure"),
主要科学问题 = c(
"每增加 1 个单位，hazard 改变多少？",
"每个 level 相对于参考组的 hazard 如何？"
),
本教程示例 = c("age", "ecog_group"),
主要设置 = c("默认", "x_exp_type = 'categorical'")
)
knitr::kable(
scale_table,
caption = "主 Cox 模型中 exposure 尺度的选择"
)| Exposure_设定 | 主要科学问题 | 本教程示例 | 主要设置 |
|---|---|---|---|
| 连续 exposure | 每增加 1 个单位，hazard 改变多少？ | age | 默认 |
| 分类 exposure | 每个 level 相对于参考组的 hazard 如何？ | ecog_group | x_exp_type = ‘categorical’ |
实践中可以记住一句话：
- 如果问题是“每增加 1 单位会怎样”，就保留连续形式
- 如果问题是“某一类和参考组相比如何”，就按分类 exposure 跑
建模前先定义好调整层次
在 cohort
研究里，模型名字最好对应科学上的调整层次，而不是只为了写代码方便。
leo_cox() 接受命名
list，这些名字会直接出现在最终结果表里。
cox_models <- list(
"Crude" = NULL,
"Model A" = c("sex"),
"Model B" = c("sex", "ecog_group")
)
cox_models
#> $Crude
#> NULL
#>
#> $`Model A`
#> [1] "sex"
#>
#> $`Model B`
#> [1] "sex" "ecog_group"这会让结果天然更像论文表，而不是后面再手工整理一遍。
主 Cox 示例 1：连续 exposure
我们先看 age 和结局之间的主关联。
res_cox <- leo_cox(
df = lung_df,
y_out = c("outcome", "outcome_censor"),
x_exp = "age",
x_cov = cox_models,
verbose = FALSE
)
res_cox
#> Exposure Outcome Case N Control N Person-years Class Crude HR
#> 1 age outcome 164 63 190.3409 Continuous 1.019
#> Crude 95% CI Crude P value Model A HR Model A 95% CI Model A P value
#> 1 1.001, 1.038 0.040 1.017 0.999, 1.036 0.061
#> Model B HR Model B 95% CI Model B P value
#> 1 1.011 0.993, 1.029 0.246主 Cox 结果应该怎么读
这张表已经很接近论文主表：
-
Case N、Control N、Person-years告诉你用了多少信息 -
Class告诉你 exposure 是按连续还是分类处理的 - 每个模型块里直接给
HR、95% CI和P value
对于连续 exposure，hazard ratio 的含义是：
exposure 每增加 1 个单位，对应的 hazard ratio 是多少
在这里就是年龄每增加 1 岁的 effect。
如果需要 tidy 表用于后续编程，可以设置 simplify = "tidy"
或者 用 simplify = FALSE 配合
leo_cox_format()：
res_cox_full <- leo_cox(
df = lung_df,
y_out = c("outcome", "outcome_censor"),
x_exp = "age",
x_cov = cox_models,
simplify = FALSE,
verbose = FALSE
)
age_tidy <- leo_cox_format(res_cox_full, style = "tidy")
age_model_b <- dplyr::filter(age_tidy, Model == "Model B")
age_model_b
#> Model Exposure Outcome Level Case N Control N Person-years HR
#> 1 Model B age outcome <NA> 164 63 190.3409 1.011
#> 95% CI P value Class
#> 1 0.993, 1.029 0.246 Continuous这个教学示例里，Model B 的估计值是
HR = 1.011，意味着在调整 sex 和 ECOG 后，年龄每增加 1
岁，incident hazard 大约增加 1.1%。但它的置信区间跨过了 1，
说明这个调整后的估计仍然不够精确，证据并不强。
这正是主 Cox 结果最标准的阅读顺序：
- 先看方向和效应大小
- 再看不确定性
- 最后再做科学上保守的解释
什么时候切换到 tidy 输出
宽表适合直接读和直接汇报。
如果你接下来还要继续做 dplyr 处理、拼接多组结果或画图，tidy
表会更方便。
dplyr::select(
age_tidy,
Model, Exposure, Outcome, HR, `95% CI`, `P value`
)
#> Model Exposure Outcome HR 95% CI P value
#> 1 Crude age outcome 1.019 1.001, 1.038 0.040
#> 2 Model A age outcome 1.017 0.999, 1.036 0.061
#> 3 Model B age outcome 1.011 0.993, 1.029 0.246可以把规则记成一句话：
-
simplify = "wide"（默认）或simplify = "tidy"直接拿结果表 -
simplify = FALSE拿完整对象，可以继续用leo_cox_format()或查看拟合结果
也可以设置 p_fmt = "scientific" 来显示精确 p 值，而不是
<0.001。
主 Cox 示例 2：分类 exposure
同一套工作流也适用于分类 exposure。
res_cox_cat <- leo_cox(
df = lung_df,
y_out = c("outcome", "outcome_censor"),
x_exp = "ecog_group",
x_cov = list("Crude" = NULL, "Model A" = c("sex")),
x_exp_type = "categorical",
verbose = FALSE
)
res_cox_cat
#> Exposure Outcome Case N Control N Person-years Class
#> 1 ECOG0 (Ref) outcome 37 26 60.6926762 Categorical (4 levels)
#> 2 ECOG1 outcome 82 31 97.2813142 Categorical (4 levels)
#> 3 ECOG2 outcome 44 6 32.0438056 Categorical (4 levels)
#> 4 ECOG3 outcome 1 0 0.3230664 Categorical (4 levels)
#> Crude HR Crude 95% CI Crude P value Model A HR Model A 95% CI
#> 1 1.000 1.000, 1.000 1.000 1.000 1.000, 1.000
#> 2 1.446 0.980, 2.134 0.063 1.519 1.028, 2.246
#> 3 2.500 1.610, 3.883 <0.001 2.579 1.660, 4.007
#> 4 9.097 1.218, 67.937 0.031 7.756 1.037, 58.039
#> Model A P value
#> 1 1.000
#> 2 0.036
#> 3 <0.001
#> 4 0.046这张表和年龄模型的读法不一样：
-
ECOG0是参考组 - 其他每一行都表示某个 level 相对
ECOG0的比较 - 所以这里不再是“每增加 1 单位”，而是 “level vs reference”
ecog_tidy <- leo_cox(
df = lung_df,
y_out = c("outcome", "outcome_censor"),
x_exp = "ecog_group",
x_cov = list("Crude" = NULL, "Model A" = c("sex")),
x_exp_type = "categorical",
simplify = "tidy",
verbose = FALSE
)
dplyr::select(
ecog_tidy,
Model, Exposure, Level, HR, `95% CI`, `P value`
)
#> Model Exposure Level HR 95% CI P value
#> 1 Crude ECOG0 (Ref) ECOG0 1.000 1.000, 1.000 1.000
#> 2 Crude ECOG1 ECOG1 1.446 0.980, 2.134 0.063
#> 3 Crude ECOG2 ECOG2 2.500 1.610, 3.883 <0.001
#> 4 Crude ECOG3 ECOG3 9.097 1.218, 67.937 0.031
#> 5 Model A ECOG0 (Ref) ECOG0 1.000 1.000, 1.000 1.000
#> 6 Model A ECOG1 ECOG1 1.519 1.028, 2.246 0.036
#> 7 Model A ECOG2 ECOG2 2.579 1.660, 4.007 <0.001
#> 8 Model A ECOG3 ECOG3 7.756 1.037, 58.039 0.046在这个示例里，Model A 提示较差的 ECOG 体能状态与更高的 incident hazard 相关：
-
ECOG1vsECOG0:HR = 1.519 -
ECOG2vsECOG0:HR = 2.579 -
ECOG3vsECOG0:HR = 7.756
但最后一行也提醒我们一个真实研究里必须警惕的问题：ECOG3
这一层在这个小示例中只有 1
个事件，因此它的置信区间极宽。稀疏单元不是表格上的小问题，而是实质性的解释风险。
Step 1 常见结果模式应该怎么解读
主 Cox 结果通常按下面这个顺序理解：
- 如果 HR 明显偏离 1，且置信区间不跨 1，那么主关联方向和统计证据都比较清楚
- 如果 HR 偏离 1，但置信区间跨 1，说明数据提示可能存在关联，但估计还不够精确
- 如果点估计接近 1，通常提示效应量可能较小；但统计证据强弱仍要结合置信区间精度、样本量和具体研究背景一起判断
- 对分类 exposure，一定要同时看效应量和每个 level 的样本/事件负担；稀疏组很容易产生看起来很大的不稳定 HR
所以 Step 1 不是“把 Cox 跑出来”就结束，而是要学会分辨：
- 信号
- 不确定性
- 稀疏数据带来的假象
论文里 Step 1 结果应该怎么写
对于连续 exposure，可以写成：
在 fully adjusted model 中，年龄每增加 1 岁，incident hazard 估计增加约 1.1% （
HR = 1.011,95% CI 0.993 to 1.029），但其置信区间包含 1。
对于分类 exposure，可以写成：
与
ECOG0相比，较差的 ECOG 体能状态与更高的 incident hazard 相关， 其中ECOG2vsECOG0的对比最强（HR = 2.579,95% CI 1.660 to 4.007）， 该估计在 sex 调整后仍然存在。
这类写法的共同点是：
- 先写清 exposure contrast
- 报
HR、95% CI和模型层次 - 结论保持保守，不过度外推