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

Question

我想弄清楚如何生成生存曲线并计算特定时间点而不是整个生存曲线的 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 天/至多 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（每一行都是一个新条目）
事件发生no/yes (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 检验。在下面的示例中，我使用了生存包中的肺部数据集。为了更好地帮助查看此方法是否合适，我将 post 这个问题进行交叉验证。

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，以便它显示对数排名分数。

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

R draw survival curve and calculate P-value at specific times

r

dataframe

survival-analysis

更新 - 使用对数秩获得特定时间点的 p 值