生存分析 Step 2：亚组分析与交互效应
survival-analysis-step2-zh.Rmd这一步要解决什么问题
这篇教程是 leo.ukb 生存分析部分的第二步。主 Cox
分析之后，最常见的 下一个问题通常是：
exposure 和结局的关联，会不会在不同人群中表现得不一样？如果不一样， 是在哪个尺度上不一样？
这个问题实际上会对应三类相关但不同的输出：
- 正式的乘性交互检验
- 用来展示各层 HR 模式的 subgroup 表
- 用于 binary-binary 对比的加性交互结果
Step 2 的重点不是把这三者混在一起，而是学会什么时候该用哪一个、结果 该怎么读。
先看数据和变量角色
这里继续使用
survival::lung，但这次我们会把数据框整理成更适合教学的
结构，而不只是为了跑模型。
lung_df <- dplyr::transmute(
survival::lung,
outcome = as.integer(status == 2),
outcome_censor = time / 365.25,
age = age,
age_high = factor(
ifelse(age >= 65, "Older", "Younger"),
levels = c("Younger", "Older")
),
sex = factor(sex, levels = c(1, 2), labels = c("Male", "Female")),
ecog_group3 = factor(
dplyr::case_when(
ph.ecog == 0 ~ "ECOG0",
ph.ecog == 1 ~ "ECOG1",
ph.ecog >= 2 ~ "ECOG2plus",
TRUE ~ NA_character_
),
levels = c("ECOG0", "ECOG1", "ECOG2plus")
),
ecog_binary = factor(
ifelse(ph.ecog >= 2, "Poor", "Good"),
levels = c("Good", "Poor")
)
)正式建模前，先像真实研究一样看变量角色和分布。
dplyr::glimpse(lung_df)
#> Rows: 228
#> Columns: 7
#> $ 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,…
#> $ age_high <fct> Older, Older, Younger, Younger, Younger, Older, Older, …
#> $ sex <fct> Male, Male, Male, Male, Male, Male, Female, Female, Mal…
#> $ ecog_group3 <fct> ECOG1, ECOG0, ECOG0, ECOG1, ECOG0, ECOG1, ECOG2plus, EC…
#> $ ecog_binary <fct> Good, Good, Good, Good, Good, Good, Poor, Poor, Good, P…| Variable | Type | Teaching role |
|---|---|---|
outcome |
二分类事件指标 | 结局事件 |
outcome_censor |
按年计的随访时间 | 结局时间 |
age |
连续 exposure | 乘性交互示例中的主 exposure |
sex |
二分类 modifier | 乘性交互中的二分类 effect modifier |
ecog_group3 |
ECOG 体能状态（三组） | subgroup 展示中的多分类 modifier |
age_high |
二分类年龄分组 | 加性交互教学用的年龄二分类变量 |
ecog_binary |
ECOG 体能状态（二分） | 加性交互教学用的二分化 ECOG 体能状态变量 |
交互分析之前，最好还要看一下各组的事件数和 binary-binary 联合分布。
| ecog_group3 | outcome | n |
|---|---|---|
| ECOG0 | 0 | 26 |
| ECOG0 | 1 | 37 |
| ECOG1 | 0 | 31 |
| ECOG1 | 1 | 82 |
| ECOG2plus | 0 | 6 |
| ECOG2plus | 1 | 45 |
| NA | 1 | 1 |
| age_high | ecog_binary | n |
|---|---|---|
| Younger | Good | 107 |
| Younger | Poor | 20 |
| Younger | NA | 1 |
| Older | Good | 69 |
| Older | Poor | 31 |
第二张表尤其重要，因为加性交互需要 binary-binary 结构，而且四个联合组 都要有样本。
建模前先决定你要回答哪个尺度上的问题
最容易避免混乱的方法，是先决定科学问题，再决定函数。
decision_table <- data.frame(
Goal = c(
"正式乘性交互检验",
"subgroup 展示",
"加性交互"
),
Exposure_types = c(
"连续、二分类、多分类都可以",
"连续、二分类、多分类都可以",
"当前建议二分类"
),
Modifier_types = c(
"连续、二分类、多分类都可以",
"二分类或多分类",
"当前建议二分类"
),
Main_function = c(
"leo_cox_interaction()",
"leo_cox_subgroup()",
"leo_cox_add_interaction()"
),
Main_question = c(
"效应是否在 HR 尺度上发生改变？",
"每一层内部的 HR 模式是什么？",
"联合效应是否超过相加预期？"
)
)
knitr::kable(
decision_table,
caption = "Step 2 里三类工具分别适合回答什么问题"
)| Goal | Exposure_types | Modifier_types | Main_function | Main_question |
|---|---|---|---|---|
| 正式乘性交互检验 | 连续、二分类、多分类都可以 | 连续、二分类、多分类都可以 | leo_cox_interaction() | 效应是否在 HR 尺度上发生改变？ |
| subgroup 展示 | 连续、二分类、多分类都可以 | 二分类或多分类 | leo_cox_subgroup() | 每一层内部的 HR 模式是什么？ |
| 加性交互 | 当前建议二分类 | 当前建议二分类 | leo_cox_add_interaction() | 联合效应是否超过相加预期？ |
可以把这张表记成一句话：
- 乘性交互最通用
- subgroup 表要求 modifier 是分类变量，因为它需要显式分层
- 加性交互最严格，当前最适合 binary-binary 对比
如果 modifier 是连续变量，leo_cox_interaction()
仍然可以做正式检验； 但 leo_cox_subgroup()
通常要等你先定义一个有临床意义的分组之后才适合用。
乘性交互：最通用的正式检验
在 Cox 模型里，乘性交互问的是：
exposure 和结局的关联，会不会随着 modifier 不同而在 HR 尺度上改变？
示例 1：连续 exposure × 二分类 modifier
这里先问：age 和结局的关联，会不会因 sex
不同而改变？
res_inter_sex <- leo_cox_interaction(
df = lung_df,
y_out = c("outcome", "outcome_censor"),
x_exp = "age",
x_inter = "sex",
x_cov = list(
"Crude" = NULL,
"Model A" = c("ecog_group3")
),
verbose = FALSE
)
res_inter_sex$result
#> Model Exposure Interaction N Case N Control N Person-years Interaction DF
#> 1 Crude age sex 227 164 63 190.3409 1
#> 2 Model A age sex 227 164 63 190.3409 1
#> P for interaction Exposure class Interaction class
#> 1 0.602 Continuous Binary
#> 2 0.180 Continuous Binary这张表要这样读：
-
Exposure class = Continuous表明age被当成连续 exposure -
Interaction class = Binary表明sex是二分类 modifier -
P for interaction是正式的异质性检验
在这个教学示例里，校正后的
P for interaction = 0.18，说明没有清楚证据
支持 age 的效应因 sex 不同而在乘法 HR
尺度上发生变化。
这并不等于“什么都没有”，而是表示：
当前数据没有足够强的证据支持 sex 对 age 效应存在乘性交互。
示例 2：连续 exposure × 多分类 modifier
同一个函数也可以处理多分类 modifier。
res_inter_ecog <- leo_cox_interaction(
df = lung_df,
y_out = c("outcome", "outcome_censor"),
x_exp = "age",
x_inter = "ecog_group3",
x_cov = list("Crude" = NULL),
verbose = FALSE
)
res_inter_ecog$result
#> Model Exposure Interaction N Case N Control N Person-years Interaction DF
#> 1 Crude age ecog_group3 227 164 63 190.3409 2
#> P for interaction Exposure class Interaction class
#> 1 0.037 Continuous Categorical (3 levels)这里 modifier
是三分类变量，P for interaction = 0.037。
这表示：
age和结局的关联，在不同 ECOG 组之间似乎并不完全一样， 至少在乘法 HR 尺度上存在异质性证据。
既然 modifier 是分类变量，下一步就很自然进入 subgroup 展示。
subgroup 分析：在正式检验之后展示 HR 模式
leo_cox_subgroup() 不是用来替代 interaction test
的。
它更像是正式检验之后的展示层，用来回答：
每一层内部的 HR 到底长什么样？
res_sub_ecog <- leo_cox_subgroup(
df = lung_df,
y_out = c("outcome", "outcome_censor"),
x_exp = "age",
x_subgroup = "ecog_group3",
x_cov = list("Crude" = NULL),
verbose = FALSE
)
res_sub_ecog$result
#> Subgroup Level N Model Exposure Outcome Exposure level Case N
#> 1 ecog_group3 ECOG0 63 Crude age outcome Per unit increase 37
#> 2 ecog_group3 ECOG1 113 Crude age outcome Per unit increase 82
#> 3 ecog_group3 ECOG2plus 51 Crude age outcome Per unit increase 45
#> Control N Person-years HR 95% CI P value P for interaction
#> 1 26 60.69268 1.059 1.013, 1.107 0.011 0.037
#> 2 31 97.28131 1.004 0.980, 1.029 0.731 0.037
#> 3 6 32.36687 0.994 0.958, 1.030 0.733 0.037这张表告诉我们的描述性模式是：
- 在
ECOG0组中，age每增加 1 岁，HR 为1.059 - 在
ECOG1组中，HR 接近于1.004 - 在
ECOG2plus组中，HR 也接近于0.994 - 同时表里附的
P for interaction = 0.037提醒我们： 这些 subgroup HR 要放在正式 interaction test 的背景下解读
所以更稳妥的论文表述应当是：
age与结局的关联在不同 ECOG 组之间存在差异（P for interaction = 0.037）。 在 subgroup 分析中，该关联在ECOG0组更强，而在ECOG1和ECOG2plus中更接近于零效应。
这比“一个 subgroup 显著、另一个不显著”要科学得多。
加性交互：回答的是另一类问题
乘性交互和加性交互不是一回事。
加性交互问的是：
两个 binary 因素一起出现时，它们的联合效应是否超过相加预期？
这也是为什么当前教学工作流里，加性交互最适合 binary-binary。
这里我们用：
-
ecog_binary：ECOG 体能状态较差 vs 较好 -
age_high：Older vs Younger
注意：这里把 age 二分类，只是为了教学加性交互。
真实研究中，这类 cutpoint 应该在分析前基于临床意义预先定义，而不是为了
让函数能跑才临时分组。
res_add <- leo_cox_add_interaction(
df = lung_df,
y_out = c("outcome", "outcome_censor"),
x_exp = "ecog_binary",
x_inter = "age_high",
x_cov = list("Crude" = NULL),
verbose = FALSE
)
res_add$result
#> Model Exposure Interaction N Case N Control N Person-years
#> 1 Crude ecog_binary age_high 227 164 63 190.3409
#> Reference group Recode applied Multiplicative interaction
#> 1 ecog_binary=Good, age_high=Younger FALSE 0.454
#> Multiplicative 95% CI P for interaction RERI RERI 95% CI AP
#> 1 0.226, 0.911 0.029 -1.475 -3.199, 0.249 -0.739
#> AP 95% CI S S 95% CI
#> 1 -1.723, 0.245 0.403 0.154, 1.056这张表可以分层来读。
先看结构列：
-
Reference group：哪一个联合组被当成参考组 -
Recode applied：函数是否为了更清楚的参考结构做了内部重编码
再看核心指标：
-
Multiplicative interaction：乘法尺度上的交互 -
RERI：relative excess risk due to interaction -
AP：attributable proportion due to interaction -
S：synergy index
这些指标的零值或无交互值并不一样：
-
Multiplicative interaction的无交互值是1 -
RERI的无交互值是0 -
AP的无交互值是0 -
S的无交互值是1
在这个教学示例里：
-
Multiplicative interaction = 0.454，且P for interaction = 0.029 -
RERI = -1.475，AP = -0.739，点估计都指向负向加性交互 - 但
RERI、AP、S的置信区间都跨过了各自的无交互值
这里最容易误读的是第一行。Multiplicative interaction
不是联合组本身的 HR，而是：
HR11 / (HR10 × HR01)
所以它小于 1 的意思不是“联合组一定是保护性的”或“联合组 HR 一定小于 1”， 而是：
联合效应低于按两个单独效应直接相乘时的预期。
举个最直观的例子：
- 如果
ecog_binary单独 HR 是2 -
age_high单独 HR 是3 - 那么在“没有乘法交互”时，联合组预期 HR 是
2 × 3 = 6 - 如果乘法交互项是
0.454，联合组 HR 大约就是2 × 3 × 0.454 = 2.724
也就是说，联合组风险仍然可能高于参考组，只是没有高到“纯乘法预期”那么强。 这更接近乘法尺度上的负向交互或拮抗，而不是“联合组 HR < 1”。
所以更谨慎的解释应当是：
该 binary-binary 对比在乘法尺度上显示出交互证据，但在加法尺度上， 虽然点估计提示可能存在负向交互，置信区间仍较宽，尚不足以支持清楚的 加性交互结论。
这正好说明了一个教学重点：
同一组数据在乘法尺度和加法尺度上的答案，可能并不完全一致。
常见结果组合应该怎么解读
实务上可以按下面 4 种情况来理解：
-
乘性交互显著，加性交互不显著
说明在 HR 尺度上有异质性证据，但加法尺度证据仍不稳定。
-
加性交互显著，乘性交互不显著
说明联合效应超过相加预期，但未见清楚的乘法尺度偏离。
-
两种尺度都显著
说明在两种尺度上都存在明确的交互证据。
-
两种尺度都不显著
说明当前数据没有清楚支持效应修饰。
但无论哪种情况，都要记住一句话：
“一个 subgroup 显著、另一个 subgroup 不显著” 本身不能证明存在交互。
在论文里应该怎么写
下面给出 3 个和上面输出对应的写法模板。
对于没有明确乘性交互的情况：
调整后未观察到
age与结局的关联在不同sex间存在显著差异 （P for interaction = 0.18）。
对于存在乘性交互并辅以 subgroup 表的情况：
age与结局的关联在不同 ECOG 组之间存在差异 （P for interaction = 0.037）。在 subgroup 分析中，ECOG0组中age每增加 1 岁的 HR 为1.059（95% CI 1.013, 1.107），而在ECOG1与ECOG2plus中接近零效应。
对于 binary-binary 的加性交互：
对于 ECOG 体能状态较差与高年龄的 binary 对比，乘法尺度交互项估计为
0.454（95% CI 0.226, 0.911），提示联合效应低于按两个单独 HR 直接 相乘时的预期；这并不表示联合组 HR 必然低于 1。加法尺度指标仍较不稳定 （例如RERI = -1.475,95% CI -3.199, 0.249）。
这种写法的好处是：
- 先交代是哪一个尺度
- 再交代方向和大小
- 最后再决定能否下比较强的结论