为ggflexsurv绘图添加阴影置信区间时遇'surv'未找到错误

问题背景 #

survminer包的ggflexsurv仅支持虚线形式的置信区间，无法直接使用阴影填充样式。修改源码得到ggflexsurvplotmod函数后，运行测试代码时出现错误：

Error in geom_ribbon():
Problem while computing aesthetics.
i Error occurred in the 5th layer.
Caused by error:
object 'surv' not found

错误原因分析 #

1. 图层美学继承问题：ggsurvplot_core生成的图层默认包含y=surv的美学映射，geom_ribbon未明确指定x变量，且fill=strata未放入aes()中，导致stat_stepribbon试图在summ数据框中寻找不存在的surv变量。
2. 辅助函数语法错误：.is_all_covariate_factor函数中存在语法错误（x$datasapply缺少分隔逗号），is_factor_or_character函数中误写is.facet应为is.factor。

修复后的完整代码 #

ggflexsurvplotmod <- function(fit, 
                               data = NULL,
                               fun = c("survival", "cumhaz"),
                               summary.flexsurv = NULL,
                               size = 1, 
                               conf.int = FALSE,
                               conf.int.flex = conf.int, 
                               conf.int.km = FALSE,
                               legend.labs = NULL,...)
  
{
  
  if (!requireNamespace("flexsurv", quietly = TRUE)){
    stop("flexsurv package needed for this function to work. Please install it.")}
  
  if(!inherits(fit, "flexsurvreg")){
    stop("Can't handle an object of class ", class(fit))
  }
  fun <- match.arg(fun)
  
  data <- .get_data(fit, data = data, complain = FALSE)
  
  summ <- .summary_flexsurv(fit, type = fun,
                            summary.flexsurv = summary.flexsurv)
  
  .strata <- summ$strata
  
  n.strata <- .strata %>% levels() %>% length()
  
  fit.ext <- .extract.survfit(fit)
  
  surv.obj <- fit.ext$surv
  
  surv.vars <- fit.ext$variables
  
  .formula <- fit.ext$formula
  
  isfac <- .is_all_covariate_factor(fit)
  
  if(!all(isfac)){
    .formula <- .build_formula(surv.obj, "1")
    n.strata <- 1
  }
  
  if(n.strata == 1 & missing(conf.int)) conf.int <- TRUE
  
  # Fit KM survival curves
  x <- do.call(survival::survfit, 
             list(formula = .formula, data = data))

  fun <- if(fun == "survival") NULL else fun

  ggsurv <- ggsurvplot_core(x, 
                          data = data, 
                          size = 0.5,
                          fun = fun, 
                          conf.int = conf.int.km,
                          legend.labs = legend.labs, ...)
  
  # Overlay the fitted models
  if(!is.null(legend.labs)){
    if(n.strata != length(legend.labs))
      stop("The length of legend.labs should be ", n.strata )
    
    summ$strata <- factor(summ$strata,
                          levels = levels(.strata),
                          labels = legend.labs)
  }
  
  time <- est <- strata <- lcl <- ucl <- NULL

  ggsurv$plot <- 
    ggsurv$plot +
    geom_line(aes(time, 
                  est, 
                  color = strata),
              data = summ, 
              size = size)
  
  # 阶梯型置信区间带的自定义统计量
  stairstepn <- function( data, direction="hv", yvars="y" ){
    direction <- match.arg(direction, c("hv", "vh") )
    data <- as.data.frame(data)[ order( data$x ), ]
    n <- nrow( data )
  
    if ( direction == "vh" ) {
      xs <- rep( 1:n, each = 2 )[ -2 * n ]
      ys <- c( 1, rep( 2:n, each = 2 ) )
    } else {
      ys <- rep( 1:n, each = 2 )[ -2 * n ]
      xs <- c( 1, rep( 2:n, each = 2))
    }
  
    data.frame(
      x = data$x[ xs ],
      data[ ys, yvars, drop=FALSE ],
      data[ xs, setdiff( names( data ), c( "x", yvars ) ), drop=FALSE ]
    ) 
  }
  
  stat_stepribbon <- function(mapping = NULL, 
                              data = NULL, 
                              geom = "ribbon", 
                              position = "identity", 
                              inherit.aes = TRUE) {
    ggplot2::layer(
      stat = Stepribbon, 
      mapping = mapping, 
      data = data, 
      geom = geom,
      position = position, 
      inherit.aes = inherit.aes
    )
  }
  
  StatStepribbon <- ggproto("stepribbon", 
                            Stat,
                            compute_group = function(., data, scales, 
                                       direction = "hv",
                                       yvars = c( "ymin", "ymax" ), ...){
                                stairstepn(data = data,
                                            direction = direction,
                                            yvars = yvars )
                              },
                            required_aes = c( "x", "ymin", "ymax" ))
  
  # 添加阴影置信区间
  if(conf.int.flex){
    ggsurv$plot <- ggsurv$plot +
      stat_stepribbon(
        aes(x = time, ymin = lcl, ymax = ucl, fill = strata),
        data = summ,
        alpha = 0.3,
        inherit.aes = FALSE  # 禁用继承，避免引入surv变量
      )
  }
  
  ggsurv
  
}

.summary_flexsurv <- function(fit, type = "survival", summary.flexsurv = NULL){
  
  summ <- summary.flexsurv
  
  if(is.null(summary.flexsurv)){
    summ <- summary(fit, type = type)}
  
  if(length(summ) == 1){
    summ <- summary(fit)[[1]] %>%
      dplyr::mutate(strata = "All")
  } else {
    .strata <- names(summ)
    summ <- purrr::pmap(list(.strata, summ), function(.s, .summ){
      dplyr::mutate(.summ, strata = .s )
    })
    summ <- dplyr::bind_rows(summ)
    summ$strata <- factor(summ$strata, levels = .strata)
  }
  
  summ
  
}

# 检查所有协变量是否为因子或字符型
.is_all_covariate_factor <- function(fit){
  x <- fit
  mf <- stats::model.frame(x)
  Xraw <- mf[,attr(mf, "covnames.orig"), drop=FALSE]
  # 修复语法错误：添加逗号分隔
  all(sapply(Xraw, is_factor_or_character))
}

is_factor_or_character <- function(x){
  # 修复函数名错误：is.facet改为is.factor
  is.factor(x) | is.character(x)
}

# 将修改后的函数替换survminer包中的原函数
assignInNamespace("ggflexsurvplot", ggflexsurvplotmod, ns = "survminer")

关键修改点 #

1. 修复辅助函数错误：
   • 修正.is_all_covariate_factor中的语法错误，添加逗号分隔x$data和sapply
   • 将is.facet改为正确的is.factor
2. 调整置信区间图层：
   • 使用stat_stepribbon并明确指定x=time，禁用美学继承（inherit.aes=FALSE），避免继承原图层的y=surv映射
   • 将fill=strata放入aes()中，确保按分层填充颜色
3. 修正函数逻辑：将fun <- match.arg(fun)移到错误判断之后，避免逻辑顺序错误

测试验证 #

运行以下代码验证修复效果：

library(flexsurv)
library(survminer)
library(ggplot2)

# 拟合模型
code_example <- flexsurvreg(formula = Surv(time, status) ~ sex, data = lung, dist = "gompertz")

# 绘制带阴影置信区间的生存曲线
ggflexsurvplotmod(code_example, conf.int = TRUE)

内容的提问来源于stack exchange，提问作者Ashley