如何生成协变量调整的 cox survival/hazard 函数？
How to generate covariate-adjusted cox survival/hazard functions?
我正在使用 survminer 包尝试为具有 5 个感兴趣子组的纵向学生级数据集生成生存和危险函数图。
我已经成功创建了一个模型，该模型显示了生存函数 ，而没有 使用 ggsurvplot.
调整学生水平的协变量
ggsurvplot(survfit(Surv(expectedgr, sped) ~ langstatus_new, data=mydata), pvalue=TRUE)
Output example
但是，我无法针对协变量调整这些曲线。我的目标是创造graphs like these。如您所见，这些是根据某些因子变量进行协变量调整的生存曲线。有谁知道如何在 R 中获得这样的图表？
您想从 Cox 模型中获取某些感兴趣协变量的某些值的生存概率，同时针对其他协变量进行调整。但是，由于我们没有对 Cox 模型中的生存时间分布做出任何假设，因此我们无法直接从中获得生存概率。我们首先必须估计基线风险函数，这通常是使用非参数 Breslow 估计器完成的。当 Cox 模型与 survival 包中的 coxph 相匹配时，我们可以通过调用 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() 背后的代码可以在 中找到。我只是简单地将函数中的 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
特定随访时间的生存概率可通过 summary() 方法的 times 参数获得。
# 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
可以从 summary().
中手动提取与这两个概率相关的 95% 置信限度
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 和 ND2 中指定的值组合的预期生存概率（以及相应的 95% 置信区间） , 我们首先需要制作 data.frames 以适当的格式包含所有信息。
# 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.frames，我们需要定义一个新函数来绘制 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 分别为 67（低）、83（中）和 100（高）。
这些生存曲线将始终满足 PH 假设，因为它们源自 Cox 模型。
注意：如果您使用 R <4.1.0
版本，请使用 function(x) 而不是 \(x)
虽然正确，但我相信 Dion Groothof 的回答中描述的方法不是通常感兴趣的方法。通常，研究人员感兴趣的是可视化变量 adjusted 对混杂因素的因果效应。简单地显示一个单一协变量组合的预测生存曲线在这里并不能真正起到作用。我建议阅读 confounder-adjusted 生存曲线。例如，参见 https://arxiv.org/abs/2203.10002。
可以使用 adjustedCurves 包在 R 中计算这些类型的曲线：https://github.com/RobinDenz1/adjustedCurves
在您的示例中，可以使用以下代码：
library(survival)
library(devtools)
# install adjustedCurves from github, load it
devtools::install_github("/RobinDenz1/adjustedCurves")
library(adjustedCurves)
# "event" needs to be binary
lung$status <- lung$status - 1
# "variable" needs to be a factor
lung$ph.ecog <- factor(lung$ph.ecog)
fit <- coxph(Surv(time, status) ~ ph.ecog + age + sex, data=lung,
x=TRUE)
# calculate and plot curves
adj <- adjustedsurv(data=lung, variable="ph.ecog", ev_time="time",
event="status", method="direct",
outcome_model=fit, conf_int=TRUE)
plot(adj)
生成以下输出：
这些生存曲线针对 age 和 sex 的影响进行了调整。有关此调整如何工作的更多信息，请参阅 adjustedCurves 包的文档或我上面引用的文章。
我正在使用 survminer 包尝试为具有 5 个感兴趣子组的纵向学生级数据集生成生存和危险函数图。
我已经成功创建了一个模型，该模型显示了生存函数 ，而没有 使用 ggsurvplot.
ggsurvplot(survfit(Surv(expectedgr, sped) ~ langstatus_new, data=mydata), pvalue=TRUE)
Output example
但是，我无法针对协变量调整这些曲线。我的目标是创造graphs like these。如您所见，这些是根据某些因子变量进行协变量调整的生存曲线。有谁知道如何在 R 中获得这样的图表？
您想从 Cox 模型中获取某些感兴趣协变量的某些值的生存概率，同时针对其他协变量进行调整。但是，由于我们没有对 Cox 模型中的生存时间分布做出任何假设，因此我们无法直接从中获得生存概率。我们首先必须估计基线风险函数，这通常是使用非参数 Breslow 估计器完成的。当 Cox 模型与 survival 包中的 coxph 相匹配时，我们可以通过调用 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() 背后的代码可以在 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
特定随访时间的生存概率可通过 summary() 方法的 times 参数获得。
# 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
可以从 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 和 ND2 中指定的值组合的预期生存概率（以及相应的 95% 置信区间） , 我们首先需要制作 data.frames 以适当的格式包含所有信息。
# 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.frames，我们需要定义一个新函数来绘制 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分别为 67（低）、83（中）和 100（高）。
这些生存曲线将始终满足 PH 假设，因为它们源自 Cox 模型。
注意：如果您使用 R <4.1.0
版本，请使用function(x) 而不是 \(x)
虽然正确，但我相信 Dion Groothof 的回答中描述的方法不是通常感兴趣的方法。通常，研究人员感兴趣的是可视化变量 adjusted 对混杂因素的因果效应。简单地显示一个单一协变量组合的预测生存曲线在这里并不能真正起到作用。我建议阅读 confounder-adjusted 生存曲线。例如，参见 https://arxiv.org/abs/2203.10002。
可以使用 adjustedCurves 包在 R 中计算这些类型的曲线：https://github.com/RobinDenz1/adjustedCurves
在您的示例中，可以使用以下代码：
library(survival)
library(devtools)
# install adjustedCurves from github, load it
devtools::install_github("/RobinDenz1/adjustedCurves")
library(adjustedCurves)
# "event" needs to be binary
lung$status <- lung$status - 1
# "variable" needs to be a factor
lung$ph.ecog <- factor(lung$ph.ecog)
fit <- coxph(Surv(time, status) ~ ph.ecog + age + sex, data=lung,
x=TRUE)
# calculate and plot curves
adj <- adjustedsurv(data=lung, variable="ph.ecog", ev_time="time",
event="status", method="direct",
outcome_model=fit, conf_int=TRUE)
plot(adj)
生成以下输出：
这些生存曲线针对 age 和 sex 的影响进行了调整。有关此调整如何工作的更多信息，请参阅 adjustedCurves 包的文档或我上面引用的文章。