Расчет прогнозируемых кривых выживаемости для непрерывных ковариаций в ggplot

Question

Как я могу построить кривые выживания для репрезентативных значений непрерывного ковариата в модели пропорциональных рисков Cox? В частности, я хотел бы сделать это в ggplot, используя объект «survfit.cox» «survfit».

Это может показаться вопросом, на который уже был дан ответ, но я просмотрел все в SO с терминами «survfit» и «newdata» (плюс множество других поисковых терминов). Это нить, которая ближе всего к ответу на мой вопрос до сих пор: 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() 

Проблема в том, что когда-то у меня есть этот dataframe называется survcurv, я не могу определить, какая из переменных «оценки» принадлежит к какому шаблону, потому что ни одна из исходных переменных не сохраняется. Например, какая из переменных «оценки» представляет собой установленную кривую для 30-летнего возраста, race = 'other', prio = '4', fin = 'no'?

Во всех других примерах, которые я видел, обычно один объект помещается в общую функцию plot() и не добавляет легенды. Я хочу использовать ggplot и добавить легенду для каждой из предсказанных кривых.

В моем собственном наборе данных модель намного сложнее, и здесь есть намного больше кривых, чем показано здесь, так как вы можете себе представить, что видят 40 различных оценок «оценка.1» .. «оценка.40». трудно понять, что к чему.

Answer 1

Спасибо за предоставление хорошо сформулированного вопроса и хороший пример. Я немного удивлен тем, что 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) 

Answer 2

Попробуйте определения ваших survcurv как это:

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

Это будет включать все значения предиктора в виде столбцов с прогнозируемыми оценками.