first commit Aalen-Johansen

DESCRIPTION

@@ -48,4 +48,4 @@ biocViews:
 Encoding: UTF-8
 Language: en-US
 LazyData: true
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(aalen_johansen)
 export(generate_data)
 export(inc_prop)
 export(one_minus_kaplan_meier)

R/aalen_johansen.R

@@ -0,0 +1,126 @@
+#' Aalen Johansen Estimator
+#'
+#' This function calculates the Aalen Johansen estimator of
+#' adverse events  observed in `[0, tau]`.
+#' Please also refer to  \insertCite{stegherr_estimating_2021;textual}{savvyr}.
+#'
+#' @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.
+#'
+#' @typed ce: number
+#'  code for competing event.
+#'
+#' @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_estimating_2021}{savvyr}
+#'
+#' @examples
+#' set.seed(123)
+#' dat <- generate_data(n = 5, cens = c(2, 5), haz_ae = 2, haz_death = 3, haz_soft = 5)
+#' aalen_johansen(dat, ce = 2, tau = 4)
+aalen_johansen <- function(data,
+                           ce,
+                           tau) {
+  assert_data_frame(data, any.missing = FALSE, min.rows = 1, min.cols = 2)
+  assert_subset(ce, c(2, 3))
+  assert_true(tau > 0)
+  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)
+
+  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]]
+    trans_mat_tau <- setmm[sum(setmm$time <= tau), ]
+    ce_prob <- trans_mat_tau$P
+    ce_prob_var <- trans_mat_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]]
+    trans_mat_tau <- setmm[sum(setmm$time <= tau), ]
+    ae_prob <- trans_mat_tau$P
+    ae_prob_var <- trans_mat_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)
+
+    trans_mat_tau_ae <- setmm[[2]][sum(setmm[[2]]$time <= tau), ]
+    ae_prob <- trans_mat_tau_ae$P
+    ae_prob_var <- trans_mat_tau_ae$var
+
+    trans_mat_tau_ce <- setmm[[3]][sum(setmm[[3]]$time <= tau), ]
+    ce_prob <- trans_mat_tau_ce$P
+    ce_prob_var <- trans_mat_tau_ce$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

inst/WORDLIST

@@ -26,3 +26,5 @@ Stegherr
 Ulm
 Universität
 VarYing
+Aalen
+Johansen

tests/testthat/test-aalen_johansen.R

@@ -0,0 +1,61 @@
+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)
+})
+
+test_that("Aalen Johansen works without events", {
+  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
+  )
+  df <- df[df$type_of_event != 1, ]
+  result <- aalen_johansen(data = df, ce = 2, tau = 4)
+  expected <- c(ae_prob = 0, ae_prob_var = 0, ce_prob = 1, ce_prob_var = 0)
+  expect_equal(result, expected, tolerance = 1e-4)
+})
+
+test_that("Aalen Johansen works without competing events", {
+  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
+  )
+  df <- df[df$type_of_event != 2, ]
+  df <- df[df$type_of_event != 3, ]
+  result <- aalen_johansen(data = df, ce = 2, tau = 4)
+  expected <- c(ae_prob = 0.5897, ae_prob_var = 0.0404, ce_prob = 0, ce_prob_var = 0)
+  expect_equal(result, expected, tolerance = 1e-4)
+})
+
+test_that("Aalen Johansen works without any events", {
+  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
+  )
+  df <- df[df$type_of_event != 1, ]
+  df <- df[df$type_of_event != 2, ]
+  df <- df[df$type_of_event != 3, ]
+  result <- aalen_johansen(data = df, ce = 2, tau = 4)
+  expected <- c(ae_prob = 0, ae_prob_var = 0, ce_prob = 0, ce_prob_var = 0)
+  expect_equal(result, expected, tolerance = 1e-4)
+})