Add possibility for left-truncation to ppc_km_overlay() and ppc_km_overlay_grouped()

NEWS.md

@@ -1,5 +1,6 @@
 # bayesplot (development version)
 
+* Add possibility for left-truncation to `ppc_km_overlay()` and `ppc_km_overlay_grouped()` by @Sakuski
 * Added `ppc_loo_pit_ecdf()` by @TeemuSailynoja
 
 # bayesplot 1.12.0

R/ppc-censoring.R

@@ -58,18 +58,33 @@
 #' \donttest{
 #' ppc_km_overlay_grouped(y, yrep[1:25, ], group = group, status_y = status_y)
 #' }
+#' # With left-truncation (delayed entry) times:
+#' condition <- y > mean(y) / 2
+#' left_truncation_y[condition] <- pmin(
+#'   runif(sum(condition), min = 0.6, max = 0.99) * y[condition],
+#'   min_vals[condition] - 0.001
+#' )
+#' \donttest{
+#' ppc_km_overlay(y, yrep[1:25, ], status_y = status_y,
+#'               left_truncation_y = left_truncation_y)
+#' }
 NULL
 
 #' @export
 #' @rdname PPC-censoring
 #' @param status_y The status indicator for the observations from `y`. This must
 #'   be a numeric vector of the same length as `y` with values in \{0, 1\} (0 =
 #'   right censored, 1 = event).
+#' @param left_truncation_y Optional parameter that specifies left-truncation
+#'   (delayed entry) times for the observations from `y`. This must
+#'   be a numeric vector of the same length as `y`. If `NULL` (default),
+#'   no left-truncation is assumed.
 ppc_km_overlay <- function(
   y,
   yrep,
   ...,
   status_y,
+  left_truncation_y = NULL,
   size = 0.25,
   alpha = 0.7
 ) {
@@ -79,8 +94,15 @@ ppc_km_overlay <- function(
   suggested_package("survival")
   suggested_package("ggfortify")
 
-  stopifnot(is.numeric(status_y))
-  stopifnot(all(status_y %in% c(0, 1)))
+  if (!is.numeric(status_y) || length(status_y) != length(y) || !all(status_y %in% c(0, 1))) {
+    stop("`status_y` must be a numeric vector of 0s and 1s the same length as `y`.")
+  }
+
+  if (!is.null(left_truncation_y)) {
+    if (!is.numeric(left_truncation_y) || length(left_truncation_y) != length(y)) {
+      stop("`left_truncation_y` must be a numeric vector of the same length as `y`.")
+    }
+  }
 
   data <- ppc_data(y, yrep, group = status_y)
 
@@ -96,7 +118,12 @@ ppc_km_overlay <- function(
                                  as.numeric(as.character(.data$group)),
                                  1))
 
-  sf_form <- survival::Surv(value, group) ~ rep_label
+  if (is.null(left_truncation_y)) {
+    sf_form <- survival::Surv(time = data$value, event = data$group) ~ rep_label
+  } else {
+    sf_form <- survival::Surv(time = left_truncation_y[data$y_id], time2 = data$value, event = data$group) ~ rep_label
+  }
+
   if (!is.null(add_group)) {
     data <- dplyr::inner_join(data,
                               tibble::tibble(y_id = seq_along(y),
@@ -164,6 +191,7 @@ ppc_km_overlay_grouped <- function(
   group,
   ...,
   status_y,
+  left_truncation_y = NULL,
   size = 0.25,
   alpha = 0.7
 ) {
@@ -175,6 +203,7 @@ ppc_km_overlay_grouped <- function(
     add_group = group,
     ...,
     status_y = status_y,
+    left_truncation_y = left_truncation_y,
     size = size,
     alpha = alpha
   )

tests/testthat/data-for-ppc-tests.R

@@ -3,6 +3,7 @@ y <- rnorm(100)
 yrep <- matrix(rnorm(2500), ncol = 100)
 group <- gl(4, 25, labels = LETTERS[1:4])
 status_y <- rep_len(0:1, length.out = length(y))
+left_truncation_y <- y - 10
 
 y2 <- rpois(30, 1)
 yrep2 <- matrix(rpois(30, 1), ncol = 30)
@@ -26,4 +27,22 @@ vdiff_loo_y <- rnorm(100, 30, 5)
 vdiff_loo_yrep <- matrix(rnorm(100 * 400, 30, 5), nrow = 400)
 vdiff_loo_lw <- vdiff_loo_yrep
 vdiff_loo_lw[] <- rnorm(100 * 400, -8, 2)
+
+
+vdiff_y3 <- rexp(50, rate = 0.2)
+vdiff_status_y3 <- rep_len(0:1, length.out = length(vdiff_y3))
+vdiff_group3 <- rep_len(c(1,2), length.out = 50)
+vdiff_left_truncation_y3 <- runif(length(vdiff_y3), min = 0, max = 0.6) * vdiff_y3
+
+simulate_truncated_exp <- function(n, rate, trunc_point) {
+  u <- runif(n)
+  return(trunc_point - log(u) / rate)
+}
+
+rate <- 0.2
+vdiff_yrep3 <- matrix(NA, nrow = 10, ncol = 50)
+for (i in 1:50) {
+  vdiff_yrep3[, i] <- simulate_truncated_exp(10, rate, vdiff_left_truncation_y3[i])
+}
+
 set.seed(seed = NULL)

tests/testthat/test-ppc-censoring.R

@@ -5,18 +5,24 @@ source(test_path("data-for-ppc-tests.R"))
 
 test_that("ppc_km_overlay returns a ggplot object", {
   skip_if_not_installed("ggfortify")
-  expect_gg(ppc_km_overlay(y, yrep, status_y = status_y, size = 0.5, alpha = 0.2))
+  expect_gg(ppc_km_overlay(y, yrep, status_y = status_y, left_truncation_y = left_truncation_y, size = 0.5, alpha = 0.2))
   expect_gg(ppc_km_overlay(y2, yrep2, status_y = status_y2))
 })
 
 test_that("ppc_km_overlay_grouped returns a ggplot object", {
   skip_if_not_installed("ggfortify")
   expect_gg(ppc_km_overlay_grouped(y, yrep, group,
-                                   status_y = status_y))
+                                   status_y = status_y,
+                                   left_truncation_y = left_truncation_y,
+                                   size = 0.5, alpha = 0.2))
   expect_gg(ppc_km_overlay_grouped(y, yrep, as.numeric(group),
-                                   status_y = status_y))
+                                   status_y = status_y,
+                                   left_truncation_y = left_truncation_y,
+                                   size = 0.5, alpha = 0.2))
   expect_gg(ppc_km_overlay_grouped(y, yrep, as.integer(group),
-                                   status_y = status_y))
+                                   status_y = status_y,
+                                   left_truncation_y = left_truncation_y,
+                                   size = 0.5, alpha = 0.2))
 
   expect_gg(ppc_km_overlay_grouped(y2, yrep2, group2,
                                    status_y = status_y2))
@@ -26,6 +32,34 @@ test_that("ppc_km_overlay_grouped returns a ggplot object", {
                                    status_y = status_y2))
 })
 
+test_that("ppc_km_overlay errors if bad status_y value", {
+  skip_if_not_installed("ggfortify")
+  expect_error(
+    ppc_km_overlay(y, yrep, status_y = FALSE),
+    "`status_y` must be a numeric vector of 0s and 1s the same length as `y`."
+  )
+  expect_error(
+    ppc_km_overlay(y, yrep, status_y = 1:10),
+    "`status_y` must be a numeric vector of 0s and 1s the same length as `y`."
+  )
+  expect_error(
+    ppc_km_overlay(y, yrep, status_y = rep(10, length(y))),
+    "`status_y` must be a numeric vector of 0s and 1s the same length as `y`."
+  )
+})
+
+test_that("ppc_km_overlay errors if bad left_truncation_y value", {
+  skip_if_not_installed("ggfortify")
+  expect_error(
+    ppc_km_overlay(y, yrep, status_y = status_y, left_truncation_y = "a"),
+    "`left_truncation_y` must be a numeric vector of the same length as `y`"
+  )
+  expect_error(
+    ppc_km_overlay(y, yrep, status_y = status_y, left_truncation_y = 1:10),
+    "`left_truncation_y` must be a numeric vector of the same length as `y`"
+  )
+})
+
 # Visual tests -----------------------------------------------------------------
 
 test_that("ppc_km_overlay renders correctly", {
@@ -44,6 +78,17 @@ test_that("ppc_km_overlay renders correctly", {
     size = 2,
     alpha = .2)
   vdiffr::expect_doppelganger("ppc_km_overlay (size, alpha)", p_custom)
+
+  p_base2 <- ppc_km_overlay(vdiff_y3, vdiff_yrep3, status_y = vdiff_status_y3)
+  vdiffr::expect_doppelganger("ppc_km_overlay (default 2)", p_base2)
+
+  p_custom2 <- ppc_km_overlay(
+    vdiff_y3,
+    vdiff_yrep3,
+    status_y = vdiff_status_y3,
+    left_truncation_y = vdiff_left_truncation_y3)
+  vdiffr::expect_doppelganger("ppc_km_overlay (left_truncation_y)",
+                              p_custom2)
 })
 
 test_that("ppc_km_overlay_grouped renders correctly", {
@@ -69,4 +114,21 @@ test_that("ppc_km_overlay_grouped renders correctly", {
     "ppc_km_overlay_grouped (size, alpha)",
     p_custom
   )
+
+  p_base2 <- ppc_km_overlay_grouped(vdiff_y3, vdiff_yrep3, vdiff_group3,
+                                   status_y = vdiff_status_y3)
+  vdiffr::expect_doppelganger("ppc_km_overlay_grouped (default 2)", p_base2)
+
+  p_custom2 <- ppc_km_overlay_grouped(
+    vdiff_y3,
+    vdiff_yrep3,
+    vdiff_group3,
+    status_y = vdiff_status_y3,
+    left_truncation_y = vdiff_left_truncation_y3
+  )
+
+  vdiffr::expect_doppelganger(
+    "ppc_km_overlay_grouped (left_truncation_y)",
+    p_custom2
+  )
 })