Создание нового Stat с выживанием :: обьект объекта выживания (NA удален из данных в compute_group) - PullRequest
Новая
       51

Создание нового Stat с выживанием :: обьект объекта выживания (NA удален из данных в compute_group)

3 голосов
/ 04 мая 2020

Я хочу создать новую статистику, которая вычисляет выживаемость с интервалом цензуры с survival::survfit.formula. Но мне кажется, что в функции compute_group неправильный фрейм данных, и я изо всех сил пытаюсь найти причину этого.

Создание фрейма данных с точно таким же кодом «снаружи» и использование geom_path (который я хочу использовать для статистики) приводит к прекрасному результату (см. Ожидаемый результат). - кажется, что survfit.formula() создает NA внутри compute_group(), но я не понимаю, почему.

установка / добавление na.rm = TRUE/FALSE ничего не меняет.

Также использование Inf вместо NA для time2 не помогает. 

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

Failing 


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).

Создано в 2020-05-04 пакетом Представить (v0.3.0)

1 Ответ

3 голосов
/ 04 мая 2020

Отличный вопрос, Тьебо, спасибо за публикацию.

Как вы уже поняли, проблема в том, что значения NA удаляются из ваших данных, прежде чем они передаются в compute_group. Расширяющая виньетка ggplot не упоминает об этом, но ваши данные сначала передаются через функцию-член compute_layer вашего объекта ggproto. Поскольку вы не определили метод 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()

enter image description here

...