r - 如何生成协变量调整的 cox 生存/风险函数？

Question

我正在使用 survminer 包尝试为具有 5 个感兴趣子组的纵向学生级数据集生成生存和危险函数图。

我已经成功创建了一个模型，该模型显示了生存函数，无需使用 ggsurvplot 调整学生级别的协变量。

ggsurvplot(survfit(Surv(expectedgr, sped) ~ langstatus_new, data=mydata), pvalue=TRUE)

Output example

但是，我无法设法根据协变量调整这些曲线。我的目标是创建graphs like these 。正如您所看到的，这些是根据某些因子变量进行协变量调整的生存曲线。有人知道如何在R中获得这样的图表吗？

Answer 1

您希望从 Cox 模型中获取某些感兴趣协变量的某些值的生存概率，同时调整其他协变量。然而，由于我们没有对 Cox 模型中的生存​​时间分布做出任何假设，因此我们无法直接从中获得生存概率。我们首先必须估计基线风险函数，这通常是使用非参数 Breslow 估计器完成的。当 Cox 模型安装 coxph 时来自survival包中，我们可以通过调用 survfit() 来获得这样的概率功能。您可以咨询?survfit.coxph了解更多信息。

让我们看看如何使用 lung 来做到这一点数据集。

library(survival)

# select covariates of interest
df <- subset(lung, select = c(time, status, age, sex, ph.karno))

# assess whether there are any missing observations
apply(df, 2, \(x) sum(is.na(x))) # 1 in ph.karno

# listwise delete missing observations
df <- df[complete.cases(df), ]

# Cox model
fit <- coxph(Surv(time, status == 2) ~ age + sex + ph.karno, data = df)

## Note that I ignore the fact that ph.karno does not satisfy the PH assumption.

# specify for which combinations of values of age, sex, and 
# ph.karno we want to derive survival probabilies
ND1 <- with(df, expand.grid(
  age = median(age),
  sex = c(1,2),
  ph.karno = median(ph.karno)
))
ND2 <- with(df, expand.grid(
  age = median(age),
  sex = 1, # males
  ph.karno = round(create_intervals(n_groups = 3L))
))

# Obtain the expected survival times
sfit1 <- survfit(fit, newdata = ND1)
sfit2 <- survfit(fit, newdata = ND2)

函数背后的代码create_intervals()可以在 this post 中找到。我只是简单地替换了speed与 ph.karno在函数中。

输出sfit1包含 ND1 中指定的协变量组合的预期中位生存时间和相应的 95% 置信区间。 .

> sfit1
Call: survfit(formula = fit, newdata = ND)

    n events median 0.95LCL 0.95UCL
1 227    164    283     223     329
2 227    164    371     320     524

通过 times 获得特定随访时间的生存概率summary() 的参数方法。

# survival probabilities at 200 days of follow-up
summary(sfit1, times = 200)

输出再次包含预期生存概率，但现在经过 200 天的随访，其中 survival1对应ND1第一行的预期生存概率，即中位数为 age 的男性和女性患者中位数 ph.karno .

> summary(sfit1, times = 200)
Call: survfit(formula = fit, newdata = ND1)

 time n.risk n.event survival1 survival2
  200    144      71     0.625     0.751

与这两个概率相关的 95% 置信限可以从 summary() 中手动提取。 .

sum_sfit <- summary(sfit1, times = 200)
sum_sfit <- t(rbind(sum_sfit$surv, sum_sfit$lower, sum_sfit$upper))
colnames(sum_sfit) <- c("S_hat", "2.5 %", "97.5 %")
# ------------------------------------------------------
> sum_sfit
      S_hat     2.5 %    97.5 %
1 0.6250586 0.5541646 0.7050220
2 0.7513961 0.6842830 0.8250914

如果您想使用ggplot描述 ND1 中指定的值组合的预期生存概率(以及相应的 95% 置信区间)和ND2 ，我们首先需要制作data.frame以适当的格式包含所有信息。

# function which returns the output from a survfit.object
# in an appropriate format, which can be used in a call
# to ggplot()
df_fun <- \(surv_obj, newdata, factor) {
  len <- length(unique(newdata[[factor]]))
  out <- data.frame(
    time = rep(surv_obj[['time']], times = len),
    n.risk = rep(surv_obj[['n.risk']], times = len),
    n.event = rep(surv_obj[['n.event']], times = len),
    surv = stack(data.frame(surv_obj[['surv']]))[, 'values'],
    upper = stack(data.frame(surv_obj[['upper']]))[, 'values'],
    lower = stack(data.frame(surv_obj[['lower']]))[, 'values']
  )
  out[, 7] <- gl(len, length(surv_obj[['time']]))
  names(out)[7] <- 'factor'
  return(out)
}

# data for the first panel (A)
df_leftPanel <- df_fun(surv_obj = sfit1, newdata = ND1, factor = 'sex')

# data for the second panel (B)
df_rightPanel <- df_fun(surv_obj = sfit2, newdata = ND2, factor = 'ph.karno')

现在我们已经定义了 data.frame s，我们需要定义一个新函数来绘制 95% CI。我们为其指定通用名称 geom_stepribbon .

library(ggplot2)

# Function for geom_stepribbon
geom_stepribbon <- function(
  mapping     = NULL,
  data        = NULL,
  stat        = "identity",
  position    = "identity",
  na.rm       = FALSE,
  show.legend = NA,
  inherit.aes = TRUE, ...) {
  layer(
    data        = data,
    mapping     = mapping,
    stat        = stat,
    geom        = GeomStepribbon,
    position    = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params      = list(na.rm = na.rm, ... )
  )
}

GeomStepribbon <- ggproto(
  "GeomStepribbon", GeomRibbon,
  extra_params = c("na.rm"),
  draw_group = function(data, panel_scales, coord, na.rm = FALSE) {
    if (na.rm) data <- data[complete.cases(data[c("x", "ymin", "ymax")]), ]
    data   <- rbind(data, data)
    data   <- data[order(data$x), ]
    data$x <- c(data$x[2:nrow(data)], NA)
    data   <- data[complete.cases(data["x"]), ]
    GeomRibbon$draw_group(data, panel_scales, coord, na.rm = FALSE)
  }
)

最后，我们可以绘制 ND1 的预期生存概率和ND2 .

yl <- 'Expected Survival probability\n'
xl <- '\nTime (days)'

# left panel
my_colours <- c('blue4', 'darkorange')
adj_colour <- \(x) adjustcolor(x, alpha.f = 0.2)
my_colours <- c(
  my_colours, adj_colour(my_colours[1]), adj_colour(my_colours[2])
)
left_panel <- ggplot(df_leftPanel,
                     aes(x = time, colour = factor, fill = factor)) + 
  geom_step(aes(y = surv), size = 0.8) + 
  geom_stepribbon(aes(ymin = lower, ymax = upper), colour = NA) +
  scale_colour_manual(name = 'Sex',
                      values = c('1' = my_colours[1],
                                 '2' = my_colours[2]),
                      labels = c('1' = 'Males',
                                 '2' = 'Females')) +
  scale_fill_manual(name = 'Sex',
                    values = c('1' = my_colours[3],
                               '2' = my_colours[4]),
                    labels = c('1' = 'Males',
                               '2' = 'Females')) +
  ylab(yl) + xlab(xl) +
  theme(axis.text = element_text(size = 12),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = 'top')

# right panel
my_colours <- c('blue4', 'darkorange', '#00b0a4')
my_colours <- c(
  my_colours, adj_colour(my_colours[1]),
  adj_colour(my_colours[2]), adj_colour(my_colours[3])
)
right_panel <- ggplot(df_rightPanel,
                      aes(x = time, colour = factor, fill = factor)) + 
  geom_step(aes(y = surv), size = 0.8) +  
  geom_stepribbon(aes(ymin = lower, ymax = upper), colour = NA) +
  scale_colour_manual(name = 'Ph.karno',
                      values = c('1' = my_colours[1],
                                 '2' = my_colours[2],
                                 '3' = my_colours[3]),
                      labels = c('1' = 'Low',
                                 '2' = 'Middle',
                                 '3' = 'High')) +
  scale_fill_manual(name = 'Ph.karno',
                    values = c('1' = my_colours[4],
                               '2' = my_colours[5],
                               '3' = my_colours[6]),
                    labels = c('1' = 'Low',
                               '2' = 'Middle',
                               '3' = 'High')) +
  ylab(yl) + xlab(xl) +
  theme(axis.text = element_text(size = 12),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = 'top')

# composite plot
library(ggpubr)
ggarrange(left_panel, right_panel,
          ncol = 2, nrow = 1,
          labels = c('A', 'B'))

输出

解释

• 图 A 显示男性和女性患者的预期生存概率中位数 age中位数 ph.karno .
• B 组显示三名男性患者的预期生存概率，中位数 age与 ph.karno s 为 67(低)、83(中)和 100(高)。

这些生存曲线将始终满足 PH 假设，因为它们是从 Cox 模型导出的。

注意:使用function(x)而不是\(x)如果您使用 R <4.1.0

版本