计算广义调整曲线的SE或CI
计算广义调整曲线的SE或CI
提问于 2019-03-29 02:26:53
如何从survminer函数计算ggadjustedcurves的可变性指数(SE或CI)？我使用的是条件方法。有没有人有什么建议或资源？
回答 1
Stack Overflow用户
发布于 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)上面的代码给出了以下输出：
我知道代码很难看，也不是最高效的，但它的工作方式对我来说是应该的，我希望它能对你或其他碰巧遇到这个问题的人有所帮助。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/55404550
复制相关文章
相似问题