R нарисовать кривую выживания и рассчитать P-значение в определенные моменты - PullRequest
Новая
       24

R нарисовать кривую выживания и рассчитать P-значение в определенные моменты

0 голосов
/ 09 ноября 2018

Я пытаюсь выяснить, как создать кривую выживания и рассчитать P-значение для конкретной временной точки, а не всей кривой выживания.

Я использую методы surv и survfit из пакетов survminer, survival для создания объекта выживания и ggsurvplot для рисования кривой и ее p-значения. 

df_surv <- Surv(time = df$diff_in_days, event = df$survivalstat)
df_survfit <- survfit(dat_surv ~ Schedule, data = df)

ggsurvplot(
  df_survfit , 
  data = df,
  pval = TRUE
)

Теперь он вычисляет значение p по всей кривой 2500+ дней. Я также хотел бы рассчитать P-значение с точными интервалами. Допустим, я хотел бы знать вероятность выживания в / до 365 дней.

Я не могу просто отключить все записи, у которых время выживания превышает x (например, 365) дней, как показано ниже. Тогда вероятность выживания падает до 0%, поскольку субъекты, у которых это событие произошло позже 365, не учитываются.

В этом событии нет никого, кроме живого больше, чем x дней. 

df <- df[df$diff_in_days <= 365, ]

Как рассчитать P-значение в определенное время по общей кривой?

dput(head(df) моего кадра данных для воспроизводимого примера. 

structure(list(diff_in_days = structure(c(2160, 84, 273, 1245, 
2175, 114), class = "difftime", units = "days"), Schedule = c(1, 
1, 1, 2, 2, 2), survivalstat = c(0, 1, 1, 0, 1, 1)), row.names = c(12L, 
28L, 33L, 38L, 58L, 62L), class = "data.frame")

Мой фрейм данных

  • UID (каждая строка является новой записью)
  • Событие произошло нет / да (0,1)
  • Целочисленное количество дней до события (если событие еще не произошло, вычисляются дни от начала мониторинга до текущего (правая цензура))

EDIT:

, используя следующий код, чтобы установить для каждого события событие 0 через 365 дней. 

dat$survivalstat <- ifelse(dat$diff_in_days > 365, 0, dat$survivalstat)

Он вычисляет значение p, но по всей кривой. По истечении 365 дней он остается горизонтальным до конца в течение 2500+ дней (поскольку никаких событий не происходит), и все эти события после 365 дней все еще учитываются, поскольку они все еще находятся на кривой. (Я предполагаю, что даже несмотря на то, что все точки данных после 365 одинаковы, они все равно влияют на значение P?)

1 Ответ

0 голосов
/ 15 ноября 2018

Если вы хотите получить значение p в определенный момент времени, вы можете выполнить z-тест в определенный момент времени.В моем примере ниже я использовал набор данных легкого из пакета выживания.Чтобы лучше понять, подходит ли этот метод, я бы опубликовал этот вопрос на перекрестной валидации.

library(survival)
library(dplyr)
library(broom)
library(ggplot2)
fit1 <- survfit(Surv(time,status)~sex,data = lung)
          #turn into df
df <- broom::tidy(fit1) 

fit_df <- df  %>% 
          #group by strata
          group_by(strata) %>% 
          #get day  of interest or day before it
          filter(time <= 365) %>% 
          arrange(time) %>% 
          # pulls last date
          do(tail(.,1))

#calculate z score based on 2 sample test at that time point
z <- (fit_df$estimate[1]-fit_df$estimate[2]) /
      (sqrt( fit_df$std.error[1]^2+ fit_df$std.error[2]^2))
#get probability of z score
pz <- pnorm(abs(z))
#get p value
pvalue <- round(2 * (1-pz),2)



ggplot(data = df,  aes(x=time, y=estimate, group=strata, color= strata)) +
  geom_line(size = 1.5)+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2)+
  geom_vline(aes(xintercept=365))+
  geom_text(aes(x = 500,y=.8,label = paste0("p = " ,pvalue) ))+
  scale_y_continuous("Survival",
                     limits = c(0,1))+
  scale_x_continuous("Time")+
  scale_color_manual(" ", values = c("grey", "blue"))+
  scale_fill_discrete(guide = FALSE)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=14),
        axis.title.x = element_text(size =14),
        axis.text.y = element_text(size = 14),
        strip.text.x = element_text(size=14),
        axis.title.y = element_blank())+
  theme_bw()

enter image description here

Обновление - получение значения p до определенного момента времени с использованием лог-ранга 

#First censor and make follow time to the time point of interest 
lung2 <- lung %>% 
          mutate(time2 = ifelse(time >= 365, 365, time),
                 status2 = ifelse(time >= 365, 1,status))
#Compute log rank test using survdiff
sdf <- survdiff(Surv(time2,status2)~sex,data = lung2)
#extract p-value
p.val <- round(1 - pchisq(sdf$chisq, length(sdf$n) - 1),3)

Вggplot Код, приведенный выше, вы можете заменить pvalue на p.val, чтобы он показывал рейтинг в журнале.

...