基于R语言绘制BMI与事件风险关联的限制性立方样条风险比图

步骤说明与代码实现 #

先清理数据缺失值，拟合带RCS的Cox模型，再生成预测值完成绘图：

# 安装并加载所需包
install.packages(c("rms", "ggplot2", "survival", "dplyr"))
library(rms)
library(ggplot2)
library(survival)
library(dplyr)

# 用户提供的示例数据集
df <- data.frame(PatientID = c("0002" ,"0002", "0005", "0005" ,"0009" ,"0009" ,"0018", "0018" ,"0039" ,"0039" , "0043" ,"0043", "0046", "0046" ,"0048" ,"0048"),
                 time = c( 1961.810 , 929.466 , 978.166, 1005.820 , 925.752 , 969.469  ,943.398 ,  965.292 , 1996.404 ,  967.047 ,  NA , 893.428 , 921.606 , 976.192 , 929.590 , 950.493),
                 event = c(1 , 1 , 0 , 1 , 0 ,  1 , 0 , 0   ,  1 , 0 , 1 ,1 , 0 ,0  ,0 , 1),
                 BMI = c( 10.140 , 20.810 , 24.466 , 31.166, 26.469 , 40.398  ,20.034,  23.292 , 50.404 , 19.610 , 20.047, 37.517 , 36.428 , 19.606 , 20.590  ,29.493), stringsAsFactors = F)

# 1. 清理数据：移除time为NA的行
df_clean <- df %>% filter(!is.na(time))

# 2. 拟合带限制性立方样条的Cox比例风险模型
# 用分位数设置4个节点，避免极端值干扰，也可直接用rcs(BMI,4)自动生成节点
knots <- quantile(df_clean$BMI, c(0.05, 0.35, 0.65, 0.95))
cox_rcs <- cph(Surv(time, event) ~ rcs(BMI, knots = knots), data = df_clean, surv = T, x = T, y = T)

# 3. 生成覆盖BMI范围的预测序列
new_bmi <- seq(min(df_clean$BMI), max(df_clean$BMI), length.out = 100)
pred_data <- data.frame(BMI = new_bmi)

# 4. 预测HR及95%置信区间，以BMI中位数为参考值（HR=1）
ref_bmi <- median(df_clean$BMI)
pred <- predict(cox_rcs, newdata = pred_data, type = "lp", se.fit = T)
ref_lp <- predict(cox_rcs, newdata = data.frame(BMI = ref_bmi), type = "lp")
pred_data$HR <- exp(pred$fit - ref_lp)
pred_data$lower <- exp((pred$fit - ref_lp) - 1.96 * pred$se.fit)
pred_data$upper <- exp((pred$fit - ref_lp) + 1.96 * pred$se.fit)

# 5. 绘制限制性立方样条图
ggplot(pred_data, aes(x = BMI, y = HR)) +
  geom_line(color = "#1f77b4", linewidth = 1) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "#1f77b4", alpha = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray50") +
  labs(x = "BMI", y = "风险比（HR）", title = "BMI与事件关联的限制性立方样条HR图") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

关键说明 #

• 节点选择：用分位数节点可降低极端BMI值的影响，节点数通常设为3-5个，可根据数据分布调整。
• 参考值：默认用BMI中位数作为HR=1的基准，也可改为临床常用值（如25）。

步骤说明与代码实现 #

将BMI按指定阈值分组，拟合Cox回归后提取HR和置信区间绘图：

# 1. 创建BMI分组变量
df_clean$bmi_group <- cut(df_clean$BMI, 
                          breaks = c(-Inf, 10, 20, 30, 40, 50, Inf),
                          labels = c("≤10", "10-20", "20-30", "30-40", "40-50", ">50"))
# 设置"10-20"为参考组，可根据需求调整
df_clean$bmi_group <- relevel(df_clean$bmi_group, ref = "10-20")

# 2. 拟合单变量Cox回归
cox_group <- coxph(Surv(time, event) ~ bmi_group, data = df_clean)

# 3. 提取HR、95%CI及分组信息
library(broom)
cox_results <- tidy(cox_group, exponentiate = T, conf.int = T) %>%
  mutate(group = stringr::str_remove(term, "bmi_group")) %>%
  select(group, estimate, conf.low, conf.high)

# 4. 绘制带置信区间的XY图
ggplot(cox_results, aes(x = group, y = estimate)) +
  geom_point(size = 3, color = "#ff7f0e") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, color = "#ff7f0e") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray50") +
  scale_y_log10() + # HR为比值，对数刻度更直观展示对称的置信区间
  labs(x = "BMI分组", y = "风险比（HR）", title = "BMI分组与事件关联的HR及95%CI") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

关键说明 #

• 分组规则：可根据数据分布调整breaks和labels，确保覆盖所有样本（示例中"≤10"组无数据，绘图时会自动忽略）。
• 对数Y轴：HR的置信区间在对数尺度下呈对称分布，能更准确反映区间范围。

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