3

如何计算ggadjustedcurvessurvminer 函数的变异指数(SE 或 CI)?我正在使用条件方法。任何人有任何输入或资源?

4

2 回答 2

3

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

可以ggsurvplot通过添加从同一包轻松计算和绘制 CI函数conf.int = TRUE。它ggsurvplot是这样的:

library(survminer)
fit <- survfit(Surv(stop, event) ~ rx, data = bladder )
ggsurvplot(fit, data = bladder, conf.int=TRUE)
于 2019-03-28T22:05:03.523 回答
2

截至 2019 年底,仍然无法使用该ggadjustedcurves函数绘制置信区间。survest然而,可以使用包中的函数从cox 模型及其相应的置信区间估计调整后的生存概率rms

使用它,ggplot2我编写了一个函数reshape2pammtools该函数能够绘制带有置信区间的调整后的生存曲线估计值,由任意数量的变量分层:

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=TRUEy=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)

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

adjusted_surv_curve

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

于 2019-12-03T14:59:59.893 回答