R语言Cox回归亚组HR计算及森林图绘制技术问询

问题1：计算不同年龄亚组的给药效应HR #

你可以选择两种常用方案实现，推荐优先使用交互项方案，可同时检验亚组效应的差异显著性：

方案1：亚组拆分独立建模 #

直接按年龄分层拆分数据集，分别拟合Cox回归得到对应亚组的给药HR：

library(survival)
# 年龄≤40岁亚组（age40=0）
mod_le40 <- coxph(Surv(time, status) ~ drug, data = df[df$age40 == 0, ])
# 提取HR、95%置信区间、p值
res_le40 <- summary(mod_le40)$coefficients[, c("exp(coef)", "lower .95", "upper .95", "Pr(>|z|)")]

# 年龄>40岁亚组（age40=1）
mod_gt40 <- coxph(Surv(time, status) ~ drug, data = df[df$age40 == 1, ])
res_gt40 <- summary(mod_gt40)$coefficients[, c("exp(coef)", "lower .95", "upper .95", "Pr(>|z|)")]

方案2：交互项模型 #

拟合带drug和亚组变量交互项的Cox模型，可同时得到亚组效应和组间差异显著性：

# 拟合交互项模型
mod_inter <- coxph(Surv(time, status) ~ drug * age40, data = df)
sum_inter <- summary(mod_inter)

# ≤40岁亚组给药HR：直接取drug主效应的指数化系数
hr_le40 <- exp(coef(mod_inter)["drug"])
# >40岁亚组给药HR：主效应系数+交互项系数后指数化
hr_gt40 <- exp(coef(mod_inter)["drug"] + coef(mod_inter)["drug:age40"])
# 亚组效应差异p值：取交互项的显著性p值
p_diff <- sum_inter$coefficients["drug:age40", "Pr(>|z|)"]

如果交互项p值<0.05，说明给药的保护/风险效应在两个年龄亚组中存在显著统计学差异。 #

问题2：多亚组效应森林图绘制 #

可以使用survminer或ggplot2实现，以下是可直接运行的全流程代码，支持同时展示age40、高血压、吸烟等多个二分类亚组的给药效应：

步骤1：安装加载依赖包 #

install.packages(c("survival", "ggplot2", "dplyr"), ask = FALSE)
library(survival)
library(ggplot2)
library(dplyr)

步骤2：批量计算所有亚组的效应值 #

# 1. 补充你需要的亚组变量（示例为高血压hyp、吸烟smoke，已有age40可直接用）
set.seed(123)
df <- df %>%
  mutate(
    hyp = rbinom(nrow(.), 1, 0.3), # 高血压：1=是，0=否
    smoke = rbinom(nrow(.), 1, 0.4) # 吸烟：1=是，0=否
  )

# 2. 定义所有要分析的亚组变量名
sub_vars <- c("age40", "hyp", "smoke")

# 3. 批量计算各亚组的HR、95%CI、p值
sub_res <- lapply(sub_vars, function(v){
  mod <- coxph(formula = as.formula(paste0("Surv(time, status) ~ drug * ", v)), data = df)
  coefs <- coef(mod)
  se <- sqrt(diag(vcov(mod)))
  cov_inter <- vcov(mod)[paste0("drug:",v), "drug"]
  # 亚组水平为0的结果
  hr0 <- exp(coefs["drug"])
  low0 <- exp(coefs["drug"] - 1.96*se["drug"])
  high0 <- exp(coefs["drug"] + 1.96*se["drug"])
  p0 <- 2*pnorm(-abs(coefs["drug"]/se["drug"]))
  # 亚组水平为1的结果
  hr1 <- exp(coefs["drug"] + coefs[paste0("drug:",v)])
  se1 <- sqrt(se["drug"]^2 + se[paste0("drug:",v)]^2 + 2*cov_inter)
  low1 <- exp((coefs["drug"] + coefs[paste0("drug:",v)]) - 1.96*se1)
  high1 <- exp((coefs["drug"] + coefs[paste0("drug:",v)]) + 1.96*se1)
  p1 <- 2*pnorm(-abs((coefs["drug"] + coefs[paste0("drug:",v)])/se1))
  # 交互项p值
  p_inter <- 2*pnorm(-abs(coefs[paste0("drug:",v)]/se[paste0("drug:",v)]))
  # 整理输出
  bind_rows(
    data.frame(var = v, level = paste0(v, "=0"), hr = hr0, low = low0, high = high0, p = p0, p_inter = p_inter),
    data.frame(var = v, level = paste0(v, "=1"), hr = hr1, low = low1, high = high1, p = p1, p_inter = p_inter)
  )
}) %>% bind_rows()

步骤3：绘制森林图 #

ggplot(sub_res, aes(x = hr, y = reorder(level, desc(level)))) +
  # HR=1参考线
  geom_vline(xintercept = 1, linetype = "dashed", color = "gray60", linewidth = 0.8) +
  # 置信区间条
  geom_errorbarh(aes(xmin = low, xmax = high), height = 0.2, linewidth = 0.8, color = "#2c3e50") +
  # HR点
  geom_point(size = 3, color = "#2980b9") +
  # 坐标轴设置，x轴用对数刻度符合HR展示惯例
  scale_x_log10(limits = c(0.2, 5), breaks = c(0.2, 0.5, 1, 2, 5)) +
  # 标签设置
  labs(x = "风险比（HR，给药组vs未给药组）", y = "亚组", title = "不同亚组给药效应森林图") +
  # 按亚组变量分面
  facet_wrap(~var, scales = "free_y", ncol = 1, strip.position = "left") +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_blank(),
    strip.text.y.left = element_text(angle = 0, hjust = 1, size = 10)
  )

得到的森林图中，每个横条对应一个亚组的HR 95%置信区间，若横条不与HR=1的参考线相交，说明该亚组的给药效应具有统计学显著性。

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