如何在R中将平滑缩放Schoenfeld残差图Y轴转换为HR

完全可行，因为风险比（HR）就是回归系数Beta(t)的指数变换（$HR(t) = exp(\beta(t))$）。你之前的代码失败是因为test_ph$y是矩阵格式（每列对应一个协变量的Beta(t)估计值），直接传入plot()函数无法正确解析。下面是两种可行的代码方案：

方案1：基础R原生绘图 #

先提取并整理数据，逐个绘制每个协变量的HR随时间变化曲线：

# 先执行PH检验，transform="identity"让时间轴保持原始刻度
test_ph <- cox.zph(my_cox_model, transform="identity")

# 将Beta(t)转成HR，并整理成数据框
hr_df <- as.data.frame(exp(test_ph$y))
hr_df$time <- test_ph$time

# 按协变量数量调整绘图布局（示例为2行2列）
par(mfrow = c(2, 2))
# 循环绘制每个协变量的HR曲线
for (var_name in setdiff(colnames(hr_df), "time")) {
  plot(hr_df$time, hr_df[[var_name]],
       type = "l", lwd = 2,
       xlab = "时间", ylab = "风险比(HR)",
       main = paste(var_name, "的HR随时间变化"),
       ylim = c(0, max(hr_df[[var_name]]) * 1.1)) # 自适应Y轴范围
  abline(h = 1, col = "red", lty = 2) # 添加HR=1的参考线（对应Beta=0）
}
# 恢复默认绘图布局
par(mfrow = c(1, 1))

方案2：用ggplot2绘制更美观的分面/叠加图 #

如果需要更整洁的可视化效果，推荐用ggplot2：

library(ggplot2)
library(tidyr)

# 先按方案1的步骤生成hr_df
test_ph <- cox.zph(my_cox_model, transform="identity")
hr_df <- as.data.frame(exp(test_ph$y))
hr_df$time <- test_ph$time

# 将宽格式数据转为长格式，方便ggplot2处理
hr_long <- hr_df %>%
  pivot_longer(cols = -time, names_to = "协变量", values_to = "HR")

# 绘制分面图（每个协变量单独一张子图）
ggplot(hr_long, aes(x = time, y = HR, color = 协变量)) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = 1, color = "red", linetype = "dashed") +
  labs(x = "时间", y = "风险比(HR)", title = "各协变量HR随时间变化趋势") +
  facet_wrap(~协变量) +
  theme_bw()

# 或者绘制叠加图（所有协变量曲线放在同一张图）
ggplot(hr_long, aes(x = time, y = HR, color = 协变量)) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = 1, color = "red", linetype = "dashed") +
  labs(x = "时间", y = "风险比(HR)", title = "各协变量HR随时间变化趋势") +
  theme_bw()

关键说明 #

• HR=1对应原模型中Beta=0的情况，红色虚线用于快速判断：如果某协变量的HR曲线持续偏离1且波动明显，说明该协变量可能不满足比例风险假设。
• transform="identity"参数保证时间轴使用原始刻度，若需要对时间做对数变换，可改为transform="log"，此时HR曲线对应对数时间下的变化趋势。

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