为 ggadjustedcurves 计算 SE 或 CI

Question

如何从 survminer 函数计算 ggadjustedcurves 的可变性指数（SE 或 CI）？我正在使用条件方法。有人有任何意见或资源吗？

Answer 1

简答：你不能。在 ggadjustedcurves 的文档中，未指定置信区间选项。现在，如果您尝试在 ggadjustedcurves 调用中使用 conf.int 参数，它不会出错，但它也不会工作。此外，在软件包的 github page 中已报告该问题，并且已请求但尚未添加此功能。

您可以通过添加 conf.int = TRUE 轻松计算和绘制来自同一包的 ggsurvplot 函数中的 CI。 ggsurvplot 是这样的：

library(survminer)
fit <- survfit(Surv(stop, event) ~ rx, data = bladder )
ggsurvplot(fit, data = bladder, conf.int=TRUE)

Answer 2

截至 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 函数定义的模型。建议的函数将通过所有层变量的所有值拆分曲线，这些变量在使用 strat() 函数的 cph 调用中定义（您还应该指定 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)

上面的代码给出了以下输出：

我知道代码很丑陋而且不是最有效的，但它对我来说是应该的，我希望它能对你或碰巧遇到这个问题的任何人有所帮助。

Answer 3

从 2022 年 4 月开始，您可以！不是 ggadjustedcurves 但有一个名为 adjustedCurves 的新包提供了各种计算：

https://github.com/RobinDenz1/adjustedCurves https://arxiv.org/abs/2203.10002

在我的测试下效果很好。

为 ggadjustedcurves 计算 SE 或 CI

Computing SE or CI for ggadjustedcurves

r

cox-regression

survival