R 绘制生存曲线并计算特定时间的 P 值

Question

我试图弄清楚如何生成生存曲线并计算特定时间点的 P 值，而不是整个生存曲线的 P 值。

我使用包survminer、survival 中的surv 和survfit 方法来创建生存对象，并使用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
)

现在它计算 2500 多天的整个曲线的 p 值。我还想以精确的间隔计算 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)
• 距事件发生的整数天数（如果尚未发生，则计算从监控开始到当前的天数（右删失））

编辑：

使用以下代码将每个人的事件发生设置为 365 天后的 0。

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

它确实计算了 p 值，但仍在整个曲线上。 365 天后，它保持水平直到 2500 多天（因为没有事件发生）结束，并且 365 天后的那些事件仍然被考虑在内，因为它们仍在曲线中。 （我假设即使 365 之后的所有数据点都相同，它们仍然会影响 P 值？）

Answer 1

如果您想要特定时间点的 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()

更新 - 使用对数秩获取 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，这样它就会显示对数排名分数。