1

背景

我想绘制生存数据集随时间变化的风险比，包括其置信区间。作为示例，我将从survival包中获取一个简化的数据集：冒号数据集。

library(survival)
library(tidyverse)

# Colon survival dataset
data <- colon %>% 
  filter(etype == 2) %>% 
  select(c(id, rx, status, time)) %>% 
  filter(rx == "Obs" | rx == "Lev+5FU") %>% 
  mutate(rx = factor(rx))

数据集包含接受治疗的患者（即“Lev+5FU”）和未接受治疗的患者（即“Obs”）。生存曲线如下：

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% CI 以绘制它们ggplot？
  2. 是否有其他R软件包可以更方便地做同样的事情？
4

2 回答 2

2

注意：认识 Dion Groothof 的修正很重要。线和 CI 并不是真正的风险比。它们是随时间变化的 log-hazard-ratios的估计值和界限。您需要取幂才能获得 HR。

这些值在从返回的结果中cox.zph：

str(cox.zph(cox))
#----------------------
List of 7
 $ table    : num [1:2, 1:3] 1.188 1.188 1 1 0.276 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "rx" "GLOBAL"
  .. ..$ : chr [1:3] "chisq" "df" "p"
 $ x        : num [1:291] 0 0.00162 0.00323 0.00485 0.00646 ...
 $ time     : num [1:291] 23 34 45 52 79 113 125 127 138 141 ...
 $ y        : num [1:291, 1] 2.09 2.1 2.1 2.1 2.11 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:291] "23" "34" "45" "52" ...
  .. ..$ : chr "rx"
 $ var      : num [1, 1] 4.11
 $ transform: chr "km"
 $ call     : language cox.zph(fit = cox)
 - attr(*, "class")= chr "cox.zph"

要使用任何范式（base、lattice 或 ggplot2）绘制图，您可以使用timex 轴、x实线和“点”处的 y

 z <-  cox.zph(cox)
 ggdf <- data.frame( unclass(z)[c("time", "x","y")])
 ggplot(data=ggdf, aes(x=time, y=-x))+ 
        geom_line()+ ylim(range(z$y))+ 
        geom_point(aes(x=time,y=z$y) )

在此处输入图像描述

让 CI 看看getAnywhere(plot.cox.zph)

xx <- x$x
yy <- x$y
df <- max(df)
nvar <- ncol(yy)
pred.x <- seq(from = min(xx), to = max(xx), length = nsmo)
#------------
if (se) {
        bk <- backsolve(qmat$qr[1:df, 1:df], diag(df))
        xtx <- bk %*% t(bk)
        seval <- ((pmat %*% xtx) * pmat) %*% rep(1, df)
        temp <- 2 * sqrt(x$var[i, i] * seval)
        yup <- yhat + temp
        ylow <- yhat - temp
        yr <- range(yr, yup, ylow)
#---------------
if (se) {
            lines(pred.x, exp(yup), col = col[2], lty = lty[2], 
              lwd = lwd[2])
            lines(pred.x, exp(ylow), col = col[2], lty = lty[2], 
              lwd = lwd[2])
            }
于 2022-01-09T01:08:49.103 回答
1

该survminer软件包将为您执行此操作：

library(survminer)
ggcoxzph(cox.zph(cox))
于 2021-03-02T11:05:03.347 回答