In vsritter/sim.etm: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(
  echo = FALSE,
  cache = TRUE,
  collapse = TRUE,
  comment = ">"
)

library(sim.etm)
library(tidyverse)

General setup

Illness-death model without recovery
States: (1) healthy, (2) ill, (3) dead
Two strata occurring with proportions 30/70
Weibull transition intensities, i.e. transition from $i$ to $j$ for an individual from strata $h$ is $$q_{ij}^{(h)}(t) = \lambda_{h,ij} \gamma_{h,ij} t^{\gamma_{h,ij} - 1},\qquad t > 0$$

Event history simulation steps

Let $(X,S) = {(x_i, s_i)}{i = 0,1,2}$, with $(x_i, s_i) \subset [0,\tau] \times {0,1,2,3}$, be numerical vectors storing an individual's transition times and status until death or censoring. Here, $\tau < \infty$ is a fixed stopping time and status values correspond to _censured (0), healthy (1), ill (2), and dead (3). Also, let $A_{ij}(t) = \int_0^t q_{ij}(u)du$ be the cumulative transition intensity from state $i$ to $j$ at time $t$ (strata indexes omitted). In order to simulate an individual's path, we do

Initialize $t_0 = 0$ and $s_0 = 1$.
Simulate $u_2,u_3 \sim U(0,1)$, independent, and compute sojourn times $$t_{1j} = \left[{\dfrac{-\log(u_j)}{\lambda_{1j}}}\right]^{1/\gamma_{1j}}, \qquad j=2,3.$$
If $t_{12} \wedge t_{13} > \tau$, set $x_1 = \tau$, $s_1 = 0$, and $x_2$ to null.
If $t_{13} < t_{12}$, set $x_1 = t_{13}$, $s_1 = 3$, and $x_2$ to null. Else, set $x_1 = t_{12}$ and $s_1 = 2$; simulate $u \sim U(0,1)$ and compute the sojourn time $$t_{23} = \left[{\dfrac{A_{23}(t_{12})-\log(u)}{\lambda_{23}}}\right]^{1/\gamma_{23}}.$$ If $t_{23} < \tau$, set $x_2 = t_{23}$ and $s_2 = 3$; otherwise, $x_2 = \tau$ and $s_2 = 0$
Simulate censoring time $c \sim U(a,b)$.
If $x_2$ is not null and $x_1 < c < x_2$, set $x_2 = c$ and $s_2 = 0$.
If $c < x_1$, set $x_1 = c$ and $s_1 = 0$.
Return $(X,S)$ where $x_i$ is not null.

Population composition I (fixed population size)

Strata 1 (30% of the population):
H to I: $\lambda_{1,12} = 1$; $\gamma_{1,12} = 1$
H to D: $\lambda_{1,13} = 1$; $\gamma_{1,13} = 0.8$
I to D: $\lambda_{1,23} = 1$; $\gamma_{1,23} = 1.5$
Strata 2 (70% of the population):
H to I: $\lambda_{2,12} = 0.5$; $\gamma_{2,12} = 1$
H to D: $\lambda_{2,13} = 0.5$; $\gamma_{2,13} = 0.8$
I to D: $\lambda_{2,23} = 0.5$; $\gamma_{2,23} = 1.5$

Population composition II (subgroups of random size)

Strata 1 (30% of the population):
H to I: $\lambda_{1,12} = 1$; $\gamma_{1,12} = 1$
H to D: $\lambda_{1,13} = 1$; $\gamma_{1,13} = 1 - 0.2x_i$
I to D: $\lambda_{1,23} = 1$; $\gamma_{1,23} = 1 + 0.5x_i$
Strata 2 (70% of the population):
H to I: $\lambda_{2,12} = 0.5$; $\gamma_{2,12} = 1$
H to D: $\lambda_{2,13} = 0.5$; $\gamma_{2,13} = 1 - 0.2x_i$
I to D: $\lambda_{2,23} = 0.5$; $\gamma_{2,23} = 1 + 0.5x_i$

with $x_i \sim \mathrm{Bernoulli}(0.5)$

Sampling

Stratified sampling with allocation proportions 80/20.

# -----------------------------------------------------------------------------
# ETM
# -----------------------------------------------------------------------------
Q <- array(NA, dim = c(3, 3, 2))
Q[1,2,] <- c(1, 1)
Q[1,3,] <- c(2, .5)
Q[2,3,] <- c(1, 1.5)

data <- simulate_process(N = 2000, gen = Q)

dt_etm <- data %>%
  arrange(id, time) %>%
  group_by(id) %>%
  mutate(from = paste0("T", state),
         to = paste0("T", c(state[-1], NA)),
         entry = time,
         exit = c(time[-1], NA)) %>%
  filter(!is.na(exit)) %>%
  select(id, from, to, entry, exit) %>%
  ungroup()

tra <- apply(Q, 1:2, function(x) sum(!is.na(x)) > 0)
n_states <- nrow(tra)
Q_index <- which(tra)
P_index <- sort(c(seq(1, n_states^2, n_states + 1), Q_index))
state_names <- paste0("T", 1:n_states)
trans <- paste0("T", rep(1:n_states, times = n_states), rep(1:n_states, each = n_states))
poss_trans <- trans[P_index]

etm_p <- etm::etm(dt_etm, state_names, tra, "cens", s=0, covariance = F)
etm_p <- data.frame(p = as.vector(etm_p$est),
                    t = rep(etm_p$time, each = n_states^2),
                    trans = rep(trans, times = length(t)),
                    stringsAsFactors = F) %>%
  filter(trans %in% poss_trans)

etm_p %>%
  ggplot(aes(x = t, y = p)) +
  geom_step() +
  # geom_line(aes(x = t, y = p), color = "red", data = p) +
  # geom_step(aes(x = t, y = p), color = "blue", data = etm_p) +
  facet_wrap(~ trans, scales = "free")