如何从cox.zph()提取真实时间尺度的平滑曲线坐标？

我希望通过对比Cox模型生成的时变风险比（HR）与cox.zph()生成的平滑曲线，评估模型中设定的时间-协变量交互项对潜在非比例风险的拟合效果。我为trt变量设置了时变效应：将所有个体的时间按事件时间拆分，用得到的time变量（计数过程格式数据）与trt构建交互项（我知道也可使用tt()公式，但此处因需绘制HR曲线而采用该步骤）。

随后我想将其与cox.zph()生成的平滑曲线对比，但原绘图的y轴压缩严重，因此想提取平滑值重新绘图。查阅cox.zph()帮助文档得知，可通过plot = FALSE返回曲线的x、y值列表，但提取的x（时间值）范围被限制在0-1区间，请问如何获取真实的观测时间？

library(survival)
library(tidyverse)
# Load veteran data and convert trt to numeric, with a reference cat = 0
vdata1 <-  veteran
vdata1$trt <-  as.numeric(vdata1$trt) - 1
vdata1$id <-  seq(1:dim(vdata1)[1])
vdata1 <-  vdata1 |> 
  select(id, everything())

# Fit basic model with treatment as only predictor
vfit1 <- coxph(Surv(time, status) ~ trt, vdata1)
summary(vfit1)

# Now refit model with trt:time interaction
# Split time at every event and create vector of unique event times
event_times <- sort(unique(with(vdata1, time[status == 1])))
# Create new df in CP form with splits at every event time
vdata2 <- survSplit(Surv(time, status) ~., vdata1, cut = event_times)
# Model
vfit2 <- coxph(Surv(tstart, time, status) ~ trt + trt:nsk(time, df = 3), vdata2)
summary(vfit2)
# Create newdat df
tdata <- expand.grid(trt = 1, time = seq(1, 1000, length = 100))
# Predict HR and plot
yhat <- predict(vfit2, newdata = tdata, se.fit = TRUE, type = "risk", reference = "zero")
tdata$fit <-  yhat$fit
ggplot(tdata, aes(x = time, y = fit)) + 
  geom_line(col = 'blue') +
  scale_x_continuous(limits = c(0, 1000), breaks = seq(0, 1000, by = 50)) +
  scale_y_continuous(limits = c(0, 2), breaks = seq(0, 2, by = 0.25)) +
  xlab("Observation Time") + ylab("HR") + ggtitle("Cox Model - Spline trt:time interaction") +
  theme_bw(base_size = 15) 

# Compare HR plot for trt to that from cox.zph
# Check prop. hazards
zph <- cox.zph(vfit1)
plot(zph, hr = T)
# Extract smoother coordinates
zph <-  plot(zph, plot = F)
coords <-  data.frame(x = zph$x, y = exp(zph$y[,1]))
# Plot
ggplot(coords, aes(x = x, y = y)) + 
  geom_line(col = 'blue') +
  scale_x_continuous(limits = c(0, 1000), breaks = seq(0, 1000, by = 50)) +
  scale_y_continuous(limits = c(0, 2), breaks = seq(0, 2, by = 0.25)) +
  xlab("Observation Time") + ylab("HR") + ggtitle("Smoother from cox.zph") +
  theme_bw(base_size = 15) 

plot(cox.zph_obj, plot=FALSE)返回的x值是标准化后的时间（范围0-1），这是cox.zph绘图函数默认的处理方式，目的是让平滑拟合更稳定。要还原为真实观测时间，只需做线性映射：

1. 从原数据集或模型中获取真实时间的最大值（或完整范围）
2. 将标准化x值映射到真实时间区间

修改后的核心代码段 #

# Compare HR plot for trt to that from cox.zph
# Check prop. hazards
zph <- cox.zph(vfit1)
plot(zph, hr = T)
# Extract smoother coordinates and convert x to real time
zph_plot <-  plot(zph, plot = F)
# 获取真实时间的最大值
max_real_time <- max(vdata1$time)
# 将标准化x值转换为真实时间（若原时间从0开始，直接乘最大值即可）
coords <-  data.frame(
  x = zph_plot$x * max_real_time, 
  y = exp(zph_plot$y[,1])
)
# Plot
ggplot(coords, aes(x = x, y = y)) + 
  geom_line(col = 'blue') +
  scale_x_continuous(limits = c(0, 1000), breaks = seq(0, 1000, by = 50)) +
  scale_y_continuous(limits = c(0, 2), breaks = seq(0, 2, by = 0.25)) +
  xlab("Observation Time") + ylab("HR") + ggtitle("Smoother from cox.zph") +
  theme_bw(base_size = 15) 

补充说明 #

如果原数据的时间不是从0开始，使用更通用的线性变换公式：

min_real_time <- min(vdata1$time)
max_real_time <- max(vdata1$time)
coords$x <- min_real_time + zph_plot$x * (max_real_time - min_real_time)

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