diff --git a/NEWS.md b/NEWS.md
index 83d34833..c4baf42a 100644
--- a/NEWS.md
+++ b/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
diff --git a/R/ppc-censoring.R b/R/ppc-censoring.R
index 2677cb58..7f39a78b 100644
--- a/R/ppc-censoring.R
+++ b/R/ppc-censoring.R
@@ -58,6 +58,18 @@
#' \donttest{
#' ppc_km_overlay_grouped(y, yrep[1:25, ], group = group, status_y = status_y)
#' }
+#' # With left-truncation (delayed entry) times:
+#' min_vals <- pmin(y, apply(yrep, 2, min))
+#' left_truncation_y <- rep(0, length(y))
+#' 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
@@ -65,11 +77,16 @@ NULL
#' @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 +96,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 +120,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 +193,7 @@ ppc_km_overlay_grouped <- function(
group,
...,
status_y,
+ left_truncation_y = NULL,
size = 0.25,
alpha = 0.7
) {
@@ -175,6 +205,7 @@ ppc_km_overlay_grouped <- function(
add_group = group,
...,
status_y = status_y,
+ left_truncation_y = left_truncation_y,
size = size,
alpha = alpha
)
diff --git a/man/PPC-censoring.Rd b/man/PPC-censoring.Rd
index 8e449ed9..736a15a9 100644
--- a/man/PPC-censoring.Rd
+++ b/man/PPC-censoring.Rd
@@ -6,9 +6,26 @@
\alias{ppc_km_overlay_grouped}
\title{PPC censoring}
\usage{
-ppc_km_overlay(y, yrep, ..., status_y, size = 0.25, alpha = 0.7)
+ppc_km_overlay(
+ y,
+ yrep,
+ ...,
+ status_y,
+ left_truncation_y = NULL,
+ size = 0.25,
+ alpha = 0.7
+)
-ppc_km_overlay_grouped(y, yrep, group, ..., status_y, size = 0.25, alpha = 0.7)
+ppc_km_overlay_grouped(
+ y,
+ yrep,
+ group,
+ ...,
+ status_y,
+ left_truncation_y = NULL,
+ size = 0.25,
+ alpha = 0.7
+)
}
\arguments{
\item{y}{A vector of observations. See \strong{Details}.}
@@ -27,6 +44,11 @@ additional advice specific to particular plots.}
be a numeric vector of the same length as \code{y} with values in \{0, 1\} (0 =
right censored, 1 = event).}
+\item{left_truncation_y}{Optional parameter that specifies left-truncation
+(delayed entry) times for the observations from \code{y}. This must
+be a numeric vector of the same length as \code{y}. If \code{NULL} (default),
+no left-truncation is assumed.}
+
\item{size, alpha}{Passed to the appropriate geom to control the appearance of
the \code{yrep} distributions.}
@@ -85,6 +107,18 @@ group <- example_group_data()
\donttest{
ppc_km_overlay_grouped(y, yrep[1:25, ], group = group, status_y = status_y)
}
+# With left-truncation (delayed entry) times:
+min_vals <- pmin(y, apply(yrep, 2, min))
+left_truncation_y <- rep(0, length(y))
+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)
+}
}
\references{
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari,
diff --git a/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-default-2.svg b/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-default-2.svg
new file mode 100644
index 00000000..792f234d
--- /dev/null
+++ b/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-default-2.svg
@@ -0,0 +1,71 @@
+
+
diff --git a/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-grouped-default-2.svg b/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-grouped-default-2.svg
new file mode 100644
index 00000000..a24e432f
--- /dev/null
+++ b/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-grouped-default-2.svg
@@ -0,0 +1,129 @@
+
+
diff --git a/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-grouped-left-truncation-y.svg b/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-grouped-left-truncation-y.svg
new file mode 100644
index 00000000..8b812eda
--- /dev/null
+++ b/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-grouped-left-truncation-y.svg
@@ -0,0 +1,129 @@
+
+
diff --git a/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-left-truncation-y.svg b/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-left-truncation-y.svg
new file mode 100644
index 00000000..9496f17d
--- /dev/null
+++ b/tests/testthat/_snaps/ppc-censoring/ppc-km-overlay-left-truncation-y.svg
@@ -0,0 +1,71 @@
+
+
diff --git a/tests/testthat/data-for-ppc-tests.R b/tests/testthat/data-for-ppc-tests.R
index 46c9bdec..8e2dfd4e 100644
--- a/tests/testthat/data-for-ppc-tests.R
+++ b/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)
diff --git a/tests/testthat/test-ppc-censoring.R b/tests/testthat/test-ppc-censoring.R
index f78639cb..611d2c3f 100644
--- a/tests/testthat/test-ppc-censoring.R
+++ b/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
+ )
})