R语言中如何用ggplot绘制含置信区间的生存数据时序风险比

背景 #

需要绘制生存数据集*随时间动态变化的风险比（Hazard Ratio）*及其配套95%置信区间，本次分析以survival包中简化的colon数据集作为示例，对应数据预处理代码如下：

library(survival)
library(tidyverse)

# 加载并预处理colon生存数据集
data <- colon %>% 
  filter(etype == 2) %>% 
  select(c(id, rx, status, time)) %>% 
  filter(rx == "Obs" | rx == "Lev+5FU") %>% 
  mutate(rx = factor(rx))

该数据集包含两个组别患者：未接受治疗的Obs组、接受Lev+5FU治疗的干预组。基础生存曲线的绘制代码如下：

fit <- survfit(Surv(time, status) ~ rx, data = data)
plot(fit)

对应绘制效果如下：

前期尝试 #

已通过cox.zph函数探索Cox模型的风险比随时间变化趋势，对应实现代码：

cox <- coxph(Surv(time, status) ~ rx, data = data)
plot(cox.zph(cox))

原生绘图效果如下：

目前期望将结果调整为ggplot风格的可视化形式，存在两个待解决问题：

1. 如何从cox.zph输出对象中提取各时间点的风险比数值、时间坐标以及对应的95%置信区间上下限，适配ggplot2的绘图数据格式要求？
2. 是否有其他R包可以更便捷地完成“带95%置信区间的时序风险比ggplot图”的绘制需求？

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