testing/scratch.R

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

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


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

library(tidyverse)
# data <- data %>%
#   arrange(id, state) %>%
#   group_by(id) %>%
#   mutate(time = cumsum(time))

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)


# -----------------------------------------------------------------------------
# TRUE
# -----------------------------------------------------------------------------
t_seq <- seq(0, 5, .001)
qt <- mcgen(t_seq, Q)

p <- illness_death_prob(t_seq[-1], Q)
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)


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

p2 %>%
  ggplot(aes(x = t, y = p)) +
  geom_line() +
  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")


# write_rds(etm_p, "etm_05b.rds")
# write_rds(p2, "p_02.rds")