如何绘制含tt()项的coxph模型的风险比（HR）

针对你用tt()构建的含时间依赖系数的Cox模型，以下是绘制karno变量时变HR曲线的具体步骤：

步骤1：拟合模型并提取核心参数 #

先运行你提供的模型代码，再提取系数和方差协方差矩阵，同时生成需要评估的时间序列：

library(survival)
library(ggplot2)

# 拟合时变系数Cox模型
fit.tt3 <- coxph(Surv(time, status) ~ trt + prior + karno + tt(karno), data = veteran,
                 tt = function(x, t, ...) x * nsk(t, knots = c(5, 20, 50, 150, 300),
                                                  Boundary.knots = FALSE))

# 提取模型系数与协方差矩阵
coefs <- coef(fit.tt3)
vc <- vcov(fit.tt3)

# 生成时间点序列：覆盖数据中所有观察时间范围
t_seq <- seq(0, max(veteran$time), by = 1)

步骤2：计算各时间点的HR及置信区间 #

编写函数计算每个时间点对应的HR和95%置信区间，批量处理所有时间点：

# 定义单时间点HR计算函数
calc_hr <- function(t) {
  # 计算时间t对应的样条项取值
  spline_val <- nsk(t, knots = c(5, 20, 50, 150, 300), Boundary.knots = FALSE)
  # 计算karno的总时变系数：基础系数 + tt项系数*样条值
  beta_total <- coefs["karno"] + coefs["tt(karno)"] * spline_val
  # 计算总系数的方差（考虑协方差）
  var_total <- vc["karno", "karno"] + vc["tt(karno)", "tt(karno)"] + 2 * vc["karno", "tt(karno)"]
  # 转换为HR并计算置信区间
  hr <- exp(beta_total)
  hr_lower <- exp(beta_total - 1.96 * sqrt(var_total))
  hr_upper <- exp(beta_total + 1.96 * sqrt(var_total))
  return(data.frame(time = t, hr = hr, lower = hr_lower, upper = hr_upper))
}

# 批量计算所有时间点的HR数据
hr_df <- do.call(rbind, lapply(t_seq, calc_hr))

步骤3：绘制时变HR曲线 #

用ggplot2绘制带置信区间的曲线，直观展示时间依赖性：

ggplot(hr_df, aes(x = time, y = hr)) +
  geom_line(color = "darkblue", linewidth = 1) +
  # 绘制95%置信区间带
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "lightblue", alpha = 0.3) +
  # 添加HR=1的参考线（代表无风险差异）
  geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
  labs(x = "时间", y = "Karno评分的风险比(HR)", 
       title = "Karno评分的时变风险比随时间变化",
       subtitle = "含95%置信区间") +
  theme_minimal()

图表解释 #

• 红色虚线代表HR=1：曲线在虚线以上时，Karno评分每增加1单位对应更高的死亡风险；在虚线以下则对应更低的死亡风险
• 浅蓝色置信区间带展示了HR估计的不确定性，区间越宽说明该时间点的估计精度越低
• 可根据科研需求调整时间序列的步长或范围，修改t_seq的参数即可

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