r - 使用 survival::survfit 对象创建新 Stat 失败(NA 从 compute_group 中的数据中删除)

Question

我想创建一个新的统计数据，它使用 survival::survfit.formula 计算区间删失生存。但是我似乎在 compute_group 函数中得到了一个错误的数据框，我很难找到它的原因。

使用完全相同的“外部”代码创建数据框并使用 geom_path(我想将其用于统计)，结果很好(参见预期结果)。 - 似乎 survfit.formula() 正在 compute_group() 中创建 NA，但我不明白为什么。

设置/添加 na.rm = TRUE/FALSE 不会改变任何东西。

对于 time2 使用 Inf 而不是 NA 也没有帮助。

library(ggplot2)
library(survival)
set.seed(42)
testdf <- data.frame(time = sample(30, replace = TRUE), time2 = c(20, 10, 10, 30, rep(NA, 26)))

fit_icens <-
  survival::survfit.formula(
    survival::Surv(time = time, time2 = time2, type = "interval2") ~ 1,
    data = testdf
  )

预期结果

path <- data.frame(time = fit_icens$time, time2= fit_icens$surv)

ggplot(path, aes(x = time, y = time2)) +
  geom_path() +
  coord_cartesian(ylim = c(0, 1))

失败

StatIcen <- ggplot2::ggproto("StatIcen", Stat,
  required_aes = c("time", "time2"),
  compute_group = function(data, scales) {
    fit_icens <-
      survival::survfit.formula(
        survival::Surv(time = data$time, time2 = data$time2, type = "interval2") ~ 1,
        data = data
      )
    path <- data.frame(x = fit_icens$time, y = fit_icens$surv)
    path
  }
)

stat_icen <- function(mapping = NULL, data = NULL, geom = "path",
                      position = "identity", show.legend = NA,
                      inherit.aes = TRUE, ...) {
  layer(
    stat = StatIcen, data = data, mapping = mapping, geom = geom,
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(...)
  )
}

ggplot(testdf, aes(time = time, time2 = time2)) +
  stat_icen()
#> Warning: Removed 26 rows containing non-finite values (stat_icen).

^{由 reprex package 创建于 2020-05-04 (v0.3.0)}

Answer 1

好问题 Tjebo，感谢发帖。

正如您已经发现的那样，问题在于 NA 值在传递给 compute_group 之前从您的数据中剥离。 Extending ggplot 小插图没有提到这一点，但您的数据首先通过 ggproto 对象的 compute_layer 成员函数传递。由于您尚未定义 compute_layer 方法，因此您的 StatIcen 类继承了 ggplot2::Stat 类的方法。

如果您在 ggplot2::Stat$compute_layer 中查看此方法的源代码，您将看到这是您的 NA 值被剥离的地方，使用remove_missing 函数，删除数据框中任何指定列中缺失值的行。据推测，如果 NA 值出现在您的 time 列中，您仍然希望删除它们，但如果它们出现在 time2 中则不需要。

所以我在这里所做的就是从 Stat$compute_layer 复制源代码并稍微调整 remove_missing 调用，然后使其成为 StatIcen 的成员:

StatIcen <- ggplot2::ggproto("StatIcen", Stat,
  required_aes = c("time", "time2"),
  compute_group = function(data, scales){
    fit_icens <- survival::survfit.formula(
      survival::Surv(time = data$time,  time2 = data$time2, 
                     type = "interval2") ~ 1, data = data)
    data.frame(x = fit_icens$time, y = fit_icens$surv)
  },
  compute_layer = function (self, data, params, layout) 
  {
    ggplot2:::check_required_aesthetics(self$required_aes, c(names(data), 
        names(params)), snake_class(self))
    data <- remove_missing(data, params$na.rm, "time", 
                           ggplot2:::snake_class(self), finite = TRUE)
    params <- params[intersect(names(params), self$parameters())]
    args <- c(list(data = quote(data), scales = quote(scales)), params)
    ggplot2:::dapply(data, "PANEL", function(data) {
        scales <- layout$get_scales(data$PANEL[1])
        tryCatch(do.call(self$compute_panel, args), 
                 error = function(e) {
            warning("Computation failed in `", 
                    ggplot2:::snake_class(self), 
                    "()`:\n", e$message, call. = FALSE)
            ggplot2:::new_data_frame()
        })
    })
  }
)

现在我们得到:

ggplot(testdf, aes(time = time, time2 = time2)) + stat_icen()