testing/aux_mixture_plot.R

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

rm(list = ls())
set.seed(323235)
library(tidyverse)

# -----------------------------------------------------------------------------
# ETM
# -----------------------------------------------------------------------------
Q1 <- array(NA, dim = c(3, 3, 2))
Q1[1,2,] <- c(1, 1)
Q1[1,3,] <- c(1, 1)
Q1[2,3,] <- c(1, 1)

Q2 <- array(NA, dim = c(3, 3, 2))
Q2[1,2,] <- c(2, 1)
Q2[1,3,] <- c(2, 1)
Q2[2,3,] <- c(2, 1)

Q = (.5*Q1 + .5*Q2)


data1 <- simulate_process(N = 5000, gen = Q1)
data2 <- simulate_process(N = 5000, gen = Q2)


data3 <- simulate_process(N = 10000, gen = Q)
dt_etm3 <- data3 %>%
  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()


dt_etm1 <- data1 %>%
  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()

dt_etm2 <- data2 %>%
  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()


dt_etm <- dt_etm1 %>%
  bind_rows(dt_etm2 %>% mutate(id = id + max(dt_etm1$id)))


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

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

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


t_seq <- seq(0, 5, .001)

p1 <- illness_death_prob(t_seq[-1], Q1)
p2 <- illness_death_prob(t_seq[-1], Q2)
p <- 0.5*p1 + 0.5*p2

p <- data.frame(p = as.vector(p),
                t = rep(t_seq[-1], each = n_states^2),
                trans = rep(trans, times = length(t_seq)-1),
                stringsAsFactors = F) %>%
  filter(trans %in% poss_trans)

etm_p %>%
  ggplot(aes(x = t, y = p)) +
  geom_step(color = "black") +
  geom_step(aes(x = t, y = p), color = "black", linetype = "dotted", data = p) +
  geom_step(aes(x = t, y = p), color = "red", data = etm_p3) +
  geom_step(aes(x = t, y = p), color = "gray60", data = etm_p1) +
  geom_step(aes(x = t, y = p), color = "gray60", data = etm_p2) +
  facet_wrap(~ trans, scales = "free")