4

如何绘制 cox 比例风险模型中连续协变量的代表值的生存曲线？具体来说，我想在 ggplot 中使用“survfit.cox”“survfit”对象来执行此操作。

这似乎是一个已经回答的问题，但我已经用“survfit”和“newdata”（以及许多其他搜索词）这两个词搜索了 SO 中的所有内容。这是迄今为止最接近回答我的问题的线程：Plot Kaplan-Meier for Cox regression

与该帖子的答案之一中提供的可重复示例保持一致：

url <- "http://socserv.mcmaster.ca/jfox/Books/Companion/data/Rossi.txt"
df <- read.table(url, header = TRUE)

library(dplyr)
library(ggplot2)
library(survival)
library(magrittr)
library(broom)

# Identifying the 25th and 75th percentiles for prio (continuous covariate)

summary(df$prio)

# Cox proportional hazards model with other covariates
# 'prio' is our explanatory variable of interest

m1 <- coxph(Surv(week, arrest) ~ 
                       fin + age + race + prio,
                     data = df)

# Creating new df to get survival predictions
# Want separate curves for the the different 'fin' and 'race'
# groups as well as the 25th and 75th percentile of prio

newdf <- df %$%
  expand.grid(fin = levels(fin), 
                    age = 30, 
                    race = levels(race), 
                    prio = c(1,4))

# Obtain the fitted survival curve, then tidy 
# into a dataframe that can be used in ggplot

survcurv <- survfit(m1, newdata = newdf) %>%
  tidy()

问题是，一旦我有了这个名为 的数据框survcurv，我就无法判断哪个“估计”变量属于哪个模式，因为没有保留任何原始变量。例如，哪个“估计”变量代表 30 岁的拟合曲线，race = 'other'，prio = '4'，fin = 'no'？

在我见过的所有其他示例中，通常将 survfit 对象放入通用plot()函数中并且不添加图例。我想使用 ggplot 并为每条预测曲线添加一个图例。

在我自己的数据集中，模型要复杂得多，曲线也比我在这里展示的要多得多，所以你可以想象看到 40 个不同的 'estimate.1'..'estimate.40' 变量很难理解什么是什么。

4

2 回答 2

4

感谢您提供一个措辞良好的问题和一个很好的例子。我有点惊讶，tidy这里在创造合理的输出方面做得相对较差。请参阅下文，了解我创建一些可绘制数据的尝试：

library(tidyr)
newdf$group <- as.character(1:nrow(newdf))

survcurv <- survfit(m1, newdata = newdf) %>%
  tidy() %>% 
  gather('key', 'value', -time, -n.risk, -n.event, -n.censor) %>% 
  mutate(group = substr(key, nchar(key), nchar(key)),
         key   = substr(key, 1, nchar(key) - 2)) %>% 
  left_join(newdf, 'group') %>% 
  spread(key, value)

并创建一个情节（也许你想geom_step改用，但不幸的是没有阶梯形丝带）：

ggplot(survcurv, aes(x = time, y = estimate, ymin = conf.low, ymax = conf.high,
                     col = race, fill = race)) +
  geom_line(size = 1) +
  geom_ribbon(alpha = 0.2, col = NA) +
  facet_grid(prio ~ fin)

在此处输入图像描述

于 2016-10-05T11:39:38.110 回答
3

尝试survcurv像这样定义你的：

survcurv <- 
  lapply(1:nrow(newdf),
         function(x, m1, newdata){
           cbind(newdata[x, ], survfit(m1, newdata[x, ]) %>% tidy)
         },
         m1, 
         newdf) %>%
  bind_rows()

这将包括所有预测值作为具有预测估计值的列。

于 2016-10-05T11:39:54.970 回答