如何在R中基于Cox回归结果绘制交互项对风险比的影响图

前置准备 #

• 需提前安装加载survival（Cox模型拟合）、ggplot2（可视化）包，也可选择安装ggeffects包简化边际效应提取流程

# 加载依赖包
library(survival)
library(ggplot2)
# 可选依赖
# install.packages("ggeffects")
library(ggeffects)

拟合模型 #

与你给出的示例模型保持一致：

fit <- coxph(Surv(time, status) ~ ph.karno*sex + ph.ecog, data=lung)

生成绘图数据并出图 #

注意：计算时会将控制变量ph.ecog固定为样本均值，排除其他变量干扰，仅展示ph.karno与sex的交互对风险比的影响

方法1：手动计算预测值（无需额外依赖包） #

# 构造覆盖ph.karno全部取值范围、sex取两个水平的预测数据集
new_data <- expand.grid(
  ph.karno = seq(min(lung$ph.karno, na.rm = TRUE), max(lung$ph.karno, na.rm = TRUE), by = 1),
  sex = c(1, 2),
  ph.ecog = mean(lung$ph.ecog, na.rm = TRUE)
)
# 预测线性项及标准误
pred_res <- predict(fit, newdata = new_data, se.fit = TRUE)
# 计算风险比及95%置信区间
new_data$hr <- exp(pred_res$fit)
new_data$hr_low <- exp(pred_res$fit - 1.96 * pred_res$se.fit)
new_data$hr_high <- exp(pred_res$fit + 1.96 * pred_res$se.fit)
# 转换sex为因子方便标注
new_data$sex <- factor(new_data$sex, levels = c(1,2), labels = c("男性", "女性"))

# 绘图
ggplot(new_data, aes(x = ph.karno, y = hr, color = sex, fill = sex)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = hr_low, ymax = hr_high), alpha = 0.2, color = NA) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray40") +
  labs(x = "ph.karno评分", y = "风险比(HR)", color = "性别", fill = "性别") +
  theme_bw()

方法2：用ggeffects包快速实现 #

# 直接提取交互项边际效应，自动固定其他变量为均值
pred_eff <- ggpredict(fit, terms = c("ph.karno [all]", "sex"))
# 绘图
plot(pred_eff) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray40") +
  labs(x = "ph.karno评分", y = "风险比(HR)", color = "性别", fill = "性别") +
  theme_bw()

两种方法输出结果完全一致，图中不同颜色的线对应不同性别，阴影部分为HR的95%置信区间，灰色虚线为HR=1的无效应参考线，你可以直接调整ggplot2的相关参数自定义图的样式。

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