如何从survminer函数计算ggadjustedcurves
的可变性指数(SE或CI)?我使用的是条件方法。有没有人有什么建议或资源?
发布于 2019-12-03 23:00:00
截至2019年底,仍无法使用ggadjustedcurves
函数绘制置信区间。然而,可以使用rms
软件包中的survest
函数从cox模型及其相应的置信区间估计调整后的生存概率。
使用这个,ggplot2
,reshape2
和pammtools
我写了一个函数,它能够用置信区间绘制调整后的生存曲线估计,按任意数量的变量分层:
library(rms)
library(ggplot2)
library(reshape2)
library(pammtools)
adjusted_surv_curves <- function(model_cph, ..., ci=TRUE, conf.int=0.95,
from=0, to=max(surv_prob$time), step=0.5,
xlab="Time", ylab="Survival Probability", title="",
labels=names(strata), legend.title="Strata",
plot.ylim=c(0, 1), theme=theme_bw(), size=1,
palette="Set1", ci.alpha=0.2, return.data=F) {
# check input
stopifnot(class(model_cph)==c("cph", "rms", "coxph"))
# check if packages are loaded
if(!(all(c("ggplot2", "reshape2", "rms", "pammtools") %in% (.packages())))){
stop("One or more of the following packages are not attached: ggplot2, reshape2, rms, pammtools")
}
# to get strata, this needs to be called without
# specified times
surv_prob <- survest(model_cph)
strata <- surv_prob$strata
# calc adjusted survival probabilities
surv_prob <- survest(model_cph, times = seq(from, to, by=step),
conf.int = conf.int)
# extract estimates from survest object
plotd_surv <- data.frame(time=surv_prob$time)
plotd_lower <- data.frame(time=surv_prob$time)
plotd_upper <- data.frame(time=surv_prob$time)
for (i in 1:length(strata)) {
plotd_surv[,ncol(plotd_surv)+1] <- surv_prob$surv[i,]
plotd_lower[,ncol(plotd_lower)+1] <- surv_prob$lower[i,]
plotd_upper[,ncol(plotd_upper)+1] <- surv_prob$upper[i,]
}
# put together in new data frame
# melt dataframes
plotd_surv <- melt(plotd_surv, id.vars = "time")
colnames(plotd_surv)[3] <- "est"
plotd_lower <- melt(plotd_lower, id.vars = "time")
colnames(plotd_lower)[3] <- "lower"
plotd_upper <- melt(plotd_upper, id.vars = "time")
colnames(plotd_upper)[3] <- "upper"
# merge to one
plotdata <- merge(plotd_surv, plotd_lower, by=c("time", "variable"))
plotdata <- merge(plotdata, plotd_upper, by=c("time", "variable"))
# return data instead, if specified
if (return.data) {
return(plotdata)
}
# plot curves
p <- ggplot(data=plotdata, aes(x=time)) +
geom_step(aes(y=est, color=variable), size=size) +
theme + ylim(plot.ylim) +
scale_colour_brewer(palette = palette, name=legend.title,
labels=labels) +
xlab(xlab) + ylab(ylab) + ggtitle(title)
# add confidence interval if not specified otherwise
if (ci) {
p <- p + pammtools::geom_stepribbon(
aes(ymin=lower, ymax=upper, fill=variable),
alpha=ci.alpha) +
scale_fill_brewer(palette = palette, name=legend.title,
labels=labels)
}
# add additional ggplot parameter, if specified
add_params <- list(...)
if (length(add_params)!=0) {
for (i in 1:length(add_params)) {
p <- p + add_params[[i]]
}
}
# return plot
return(p)
}
作为输入,由cph
函数定义的模型是必需的。建议的函数将根据在cph
调用中使用strat()
函数定义的所有地层变量的所有值分割曲线(还应指定x=TRUE
和y=TRUE
以进行正确估计)。to
、from
和step
参数可用于指定应计算生存概率的间隔和频率。绘图本身的各种参数可以随意修改,因为它还接受ggplot
调用的任何其他参数。
注意:生存概率及其对应的置信区间的估计值与ggadjustedcurves
函数的估计值不同,因为它们的计算方法不同。有关如何计算置信区间的更多信息,可以参考官方文档:https://rdrr.io/cran/rms/man/survest.cph.html
例如,假设已经定义了函数:
library(rms)
library(ggplot2)
library(reshape2)
library(pammtools)
# load example dataset
data("lung")
# define model
model <- cph(with(data=lung, Surv(time, status) ~ strat(sex) + age + ph.ecog), x=TRUE, y=TRUE)
# plot adjusted survival curves by sex
adjusted_surv_curves(model)
上面的代码给出了以下输出:
我知道代码很难看,也不是最高效的,但它的工作方式对我来说是应该的,我希望它能对你或其他碰巧遇到这个问题的人有所帮助。
https://stackoverflow.com/questions/55404550
复制相似问题