first commit Aalen-Johansen

R/aalen_johansen.R

@@ -0,0 +1,114 @@
+#' Aalen Johansen Estimator
+#'
+#' This function calculates the Aalen Johansen estimator of
+#' adverse events  observed in `[0, tau]`.
+#' Please also refer to  \insertCite{stegherr_meta_analytic_2021;textual}{savvyr}.
+#'
+#' Estimating and comparing adverse event probabilities in the presence of varying follow-up times and competing events
+#' Regina Stegherr, Claudia Schmoor, Michael Lübbert, Tim Friede, Jan Beyersmann
+#' First published: 18 May 2021 https://doi.org/10.1002/pst.2130
+#'
+#' @typed data: data.frame
+#'  with columns including
+#'  - `time_to_event`: Time to the first AE, death or soft competing event.
+#'  - `type_of_event`: 0 for censored, 1 for AE, 2 for death, 3 for soft competing event.
+#'
+#' @typed tau: number
+#'  milestone at which Aalen-Johansen is computed.
+#'
+#' @typedreturn vector
+#'  with the following entries:
+#'
+#' - `ae_prob`: Estimated probability of AE.
+#' - `ae_prob_var`: Variance of that estimate.
+#' - `ce_prob`: Estimated probability of competing events.
+#' - `ce_prob_var`: Variance of competing events.
+#'
+#' @export
+#'
+#' @references
+#' \insertRef{stegherr_meta_analytic_2021}{savvyr}
+#'
+#' @examples
+#' set.seed(123)
+#' dat <- generate_data(n = 5, cens = c(2, 5), haz_ae = 2, haz_death = 3, haz_soft = 5)
+#' one_minus_kaplan_meier(dat, tau = 4)
+aalen_johansen <- function(data,
+                           ce,
+                           tau) {
+  assert_data_frame(data, any.missing = FALSE, min.rows = 1, min.cols = 2)
+  assert_numeric(data$time_to_event, lower = 0, finite = TRUE)
+  assert_integerish(data$type_of_event, any.missing = FALSE)
+  assert_subset(data$type_of_event, c(0, 1, 2, 3))
+  assert_number(tau, finite = TRUE)
+  assert_true(tau > 0)
+  assert_subset(ce, c(2, 3))
+
+  data$type_of_event_accounted <- ifelse(ce == 2 & data$type_of_event == 3, 0,
+    ifelse(ce == 3 & data$type_of_event == 3, 2, data$type_of_event)
+  )
+
+  time <- data$time_to_event
+  type2 <- data$type_of_event_accounted
+
+  # conditions
+  c1 <- sum(data$type_of_event_accounted == 1)
+  c2 <- sum(data$type_of_event_accounted == 2)
+
+  if (c1 == 0) {
+    ae_prob <- 0
+    ae_prob_var <- 0
+  }
+
+  if (c2 == 0) {
+    ce_prob <- 0
+    ce_prob_var <- 0
+  }
+
+  # define auxiliary objects
+  help <- data.frame(id = data$id)
+  help$from <- 0
+  help$time <- ifelse(time == 0, 0.001, time)
+  tra <- matrix(FALSE, 2, 2)
+  tra[1, 2] <- TRUE
+  state.names <- as.character(0:1)
+
+  if (c1 == 0 & c2 != 0) {
+    help$to <- ifelse(type2 != 2, "cens", type2 - 1)
+    etmmm <- etm::etm(help, state.names, tra, "cens", s = 0)
+    setmm <- summary(etmmm)[[2]]
+    ce_prob <- setmm[sum(setmm$time <= tau), ]$P
+    ce_prob_var <- setmm[sum(setmm$time <= tau), ]$var
+  }
+
+  if (c1 != 0 & c2 == 0) {
+    help$to <- ifelse(type2 != 1, "cens", type2)
+    etmmm <- etm::etm(help, state.names, tra, "cens", s = 0)
+    setmm <- summary(etmmm)[[2]]
+
+    ae_prob <- setmm[sum(setmm$time <= tau), ]$P
+    ae_prob_var <- setmm[sum(setmm$time <= tau), ]$var
+  }
+
+  if (c1 != 0 & c2 != 0) {
+    help$to <- ifelse(!(type2 %in% c(1, 2)), "cens", type2)
+
+    tra <- matrix(FALSE, 3, 3)
+    tra[1, 2:3] <- TRUE
+    state.names <- as.character(0:2)
+    etmmm <- etm::etm(help, state.names, tra, "cens", s = 0)
+    setmm <- summary(etmmm)
+
+    ae_prob <- setmm[[2]][sum(setmm[[2]]$time <= tau), ]$P
+    ae_prob_var <- setmm[[2]][sum(setmm[[2]]$time <= tau), ]$var
+
+    ce_prob <- setmm[[3]][sum(setmm[[3]]$time <= tau), ]$P
+    ce_prob_var <- setmm[[3]][sum(setmm[[3]]$time <= tau), ]$var
+  }
+
+
+  c(
+    "ae_prob" = ae_prob, "ae_prob_var" = ae_prob_var, "ce_prob" = ce_prob,
+    "ce_prob_var" = ce_prob_var
+  )
+}

_pkgdown.yml

@@ -16,3 +16,4 @@ navbar:
         - prop_trans_inc_dens
         - prop_trans_inc_dens_ce
         - one_minus_kaplan_meier
+        - aalen_johansen

tests/testthat/test-aalen_johansen.R

@@ -0,0 +1,13 @@
+test_that("Aalen Johansen works as expected", {
+  set.seed(23)
+  df <- generate_data(
+    n = 25,
+    cens = c(0.2, 3),
+    haz_ae = 0.2,
+    haz_death = 0.3,
+    haz_soft = 0.5
+  )
+  result <- aalen_johansen(data = df, ce = 2, tau = 4)
+  expected <- c(ae_prob = 0.2719, ae_prob_var = 0.0119, ce_prob = 0.7281, ce_prob_var = 0.0119)
+  expect_equal(result, expected, tolerance = 1e-4)
+})