R/ppc_km_nada.R

In bentrueman/bgamcar1: Fit Bayesian GAMs With CAR(1) Errors To Censored Data

#' Generate ECDFs for posterior predictive checks with left-censored data
#'
#' @param input A dataframe for which to generate model predictions.
#' @param object A `brms` model object.
#' @param draw_ids Draw IDs from model object.
#' @param car1 Logical. Add CAR(1) errors?
#' @param seed Set a random seed.
#' @param ... Arguments passed on to `add_pred_draws_car1()`.
#'
#' @return A dataframe containing ECDFs generated by `NADA::cenfit()` for observations and model predictions.
#' @importFrom NADA cenfit flip Cen summary
#' @importFrom tibble tibble
#' @importFrom rlang .data
#' @importFrom dplyr filter %>% group_by ungroup mutate bind_rows select pull
#' @importFrom withr with_seed
#' @importFrom tidybayes add_predicted_draws
#' @importFrom purrr map
#' @importFrom tidyr nest unnest
#' @export
#'
#' @examples
#' library("brms")
#' seed <- 1
#' data <- read.csv(paste0(system.file("extdata", package = "bgamcar1"), "/data.csv"))
#' fit <- fit_stan_model(
#'    paste0(system.file("extdata", package = "bgamcar1"), "/test"),
#'    seed,
#'    bf(y | cens(ycens, y2 = y2) ~ 1),
#'    data,
#'    prior(normal(0, 1), class = Intercept),
#'    car1 = FALSE,
#'    save_warmup = FALSE,
#'    chains = 3
#'  )
#'  ppc_km_nada(data, fit, draw_ids = 34, seed = seed, car1 = FALSE)
ppc_km_nada <- function(input,
                        object,
                        draw_ids = 1:200,
                        car1 = TRUE,
                        seed,
                        ...) {
  cenfit_summ <- NULL
  data <- NULL

  varnames <- extract_resp(object)

  if (is.na(varnames$cens)) stop("Model formula does not include censoring")

  cenfit_in <- tibble(
    y = pull(input, .data[[varnames$resp]]),
    ci = pull(input, .data[[varnames$cens]]) == "left"
  ) %>%
    filter(!is.na(.data$y))

  cenfit_brms <- NADA::cenfit(
    obs = cenfit_in$y,
    censored = cenfit_in$ci
  )

  preds <- if (car1) {
    withr::with_seed(
      seed,
      {
        add_pred_draws_car1(input, object, draw_ids = draw_ids, type = "prediction", ...)
      }
    )
  } else {
    tidybayes::add_predicted_draws(input, object, draw_ids = draw_ids, seed = seed)
  }

  pp <- preds %>%
    ungroup() %>%
    group_by(.data$.draw) %>%
    nest() %>%
    ungroup() %>%
    mutate(
      cenfit = map(
        .data$data,
        ~ with(
          .x,
          NADA::cenfit(.prediction, censored = rep(FALSE, length(.prediction)))
        )
      ),
      cenfit_summ = map(cenfit, NADA::summary)
    ) %>%
    unnest(cenfit_summ) %>%
    select(-c(data, cenfit))

  bind_rows(
    "Posterior draws" = pp,
    "Observations" = summary(cenfit_brms),
    .id = "type"
  )
}