BIOS 667 - Lecture 19: Transition (Markov) Models

# Hypothetical transition probabilities
pi_0 <- 0.15  # P(sick | was healthy)
pi_1 <- 0.70  # P(sick | was sick)

P <- matrix(c(1 - pi_0, pi_0,
              1 - pi_1, pi_1),
            nrow = 2, byrow = TRUE,
            dimnames = list(c("Healthy", "Sick"),
                            c("Healthy", "Sick")))
P

suppressPackageStartupMessages({
  library(ggplot2)
  library(dplyr)
  library(tidyr)
})

trans_df <- as.data.frame(P) |>
  tibble::rownames_to_column("From") |>
  pivot_longer(-From, names_to = "To", values_to = "Probability")

ggplot(trans_df, aes(To, From, fill = Probability)) +
  geom_tile(color = "white", linewidth = 1) +
  geom_text(aes(label = round(Probability, 2)), size = 6) +
  scale_fill_gradient(low = "white", high = "steelblue") +
  labs(title = "Transition Probability Matrix",
       x = "To State", y = "From State") +
  theme_minimal(base_size = 14) +
  theme(legend.position = "none")

# Solve for stationary distribution
# pi = pi * P means (I - P')pi = 0
# Add constraint: sum(pi) = 1
A <- rbind(t(diag(2) - P), c(1, 1))
b <- c(0, 0, 1)
pi_star <- qr.solve(A, b)
names(pi_star) <- c("Healthy", "Sick")
round(pi_star, 3)

set.seed(667)

N <- 200      # subjects
n_times <- 6  # observations per subject

# Subject-level covariates
trt <- rep(rbinom(N, 1, 0.5), each = n_times)
id <- rep(1:N, each = n_times)
time <- rep(0:(n_times-1), N)

# Initialize
y <- numeric(N * n_times)

# Parameters
beta_0 <- -1.5   # baseline log-odds
beta_1 <- 1.8    # lag effect (strong positive dependence)
gamma_trt <- -0.6  # treatment effect
gamma_time <- 0.1  # time effect

# Generate first observation (use marginal model)
first_idx <- seq(1, N * n_times, by = n_times)
p_first <- plogis(beta_0 + gamma_trt * trt[first_idx])
y[first_idx] <- rbinom(N, 1, p_first)

# Generate subsequent observations using transition model
for (i in 1:N) {
  for (j in 2:n_times) {
    idx <- (i - 1) * n_times + j
    idx_prev <- idx - 1
    lp <- beta_0 + beta_1 * y[idx_prev] + gamma_trt * trt[idx] + gamma_time * time[idx]
    y[idx] <- rbinom(1, 1, plogis(lp))
  }
}

dat <- data.frame(
  id = factor(id),
  time = time,
  trt = factor(trt, labels = c("Control", "Treatment")),
  y = y
)

head(dat, 12)

# Add lagged response within each subject
dat <- dat |>
  group_by(id) |>
  mutate(y_lag = lag(y, 1)) |>
  ungroup()

# Remove first observation (no lag available)
dat_model <- dat |> filter(!is.na(y_lag))

head(dat_model, 12)

fit_trans <- glm(y ~ y_lag + trt + time,
                 data = dat_model,
                 family = binomial(link = "logit"))

summary(fit_trans)

suppressPackageStartupMessages(library(broom))

tidy_trans <- tidy(fit_trans, conf.int = TRUE, exponentiate = TRUE)
tidy_trans |>
  mutate(across(where(is.numeric), ~round(., 3))) |>
  select(term, estimate, conf.low, conf.high, p.value)

newdat <- expand.grid(
  y_lag = c(0, 1),
  trt = c("Control", "Treatment"),
  time = 3  # mid-study
)

newdat$prob <- predict(fit_trans, newdata = newdat, type = "response")

ggplot(newdat, aes(x = factor(y_lag), y = prob, fill = trt)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.6) +
  labs(title = "Transition Probabilities by Treatment",
       x = "Previous Response (Y_lag)",
       y = "P(Y = 1 | Y_lag, Treatment)",
       fill = "Group") +
  scale_y_continuous(limits = c(0, 1)) +
  theme_minimal(base_size = 14)

trans_table <- dat_model |>
  group_by(trt, y_lag, y) |>
  summarise(n = n(), .groups = "drop") |>
  group_by(trt, y_lag) |>
  mutate(prop = n / sum(n)) |>
  filter(y == 1) |>
  select(trt, y_lag, prop)

trans_table

suppressPackageStartupMessages(library(sandwich))
suppressPackageStartupMessages(library(lmtest))

# Cluster-robust standard errors
robust_se <- coeftest(fit_trans, vcov = vcovCL(fit_trans, cluster = dat_model$id))
robust_se

suppressPackageStartupMessages(library(geepack))

fit_gee <- geeglm(y ~ y_lag + trt + time,
                  data = dat_model,
                  id = id,
                  family = binomial,
                  corstr = "independence")

summary(fit_gee)

se_comparison <- data.frame(
  term = names(coef(fit_trans)),
  GLM_SE = summary(fit_trans)$coefficients[, "Std. Error"],
  Robust_SE = robust_se[, "Std. Error"],
  GEE_SE = summary(fit_gee)$coefficients[, "Std.err"]
)

se_comparison[, -1] <- round(se_comparison[, -1], 4)
se_comparison

# Create second lag
dat_order2 <- dat |>
  group_by(id) |>
  mutate(
    y_lag1 = lag(y, 1),
    y_lag2 = lag(y, 2)
  ) |>
  ungroup() |>
  filter(!is.na(y_lag2))

fit_order2 <- glm(y ~ y_lag1 + y_lag2 + trt + time,
                  data = dat_order2,
                  family = binomial)

summary(fit_order2)$coefficients

# Refit first-order on same data
fit_order1 <- glm(y ~ y_lag1 + trt + time,
                  data = dat_order2,
                  family = binomial)

anova(fit_order1, fit_order2, test = "LRT")

fit_tv <- glm(y ~ y_lag * time + trt,
              data = dat_model,
              family = binomial)

summary(fit_tv)$coefficients

times <- 0:5
coefs <- coef(fit_tv)

# Effect of y_lag at each time point
effect_at_t <- coefs["y_lag"] + coefs["y_lag:time"] * times

eff_df <- data.frame(time = times, log_or = effect_at_t)

ggplot(eff_df, aes(time, log_or)) +
  geom_line(linewidth = 1.2, color = "steelblue") +
  geom_point(size = 3, color = "steelblue") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Time-Varying Dependence on Previous Response",
       x = "Time", y = "Log Odds Ratio for Y_lag") +
  theme_minimal(base_size = 14)

suppressPackageStartupMessages(library(lme4))

# Transition model (already fit)
# fit_trans

# Marginal model (GEE)
fit_marginal <- geeglm(y ~ trt + time,
                       data = dat_model,
                       id = id,
                       family = binomial,
                       corstr = "ar1")

# Mixed model (GLMM with random intercept)
fit_mixed <- glmer(y ~ trt + time + (1 | id),
                   data = dat_model,
                   family = binomial)

# Compare treatment effects
data.frame(
  Model = c("Transition", "Marginal (GEE)", "Mixed (GLMM)"),
  Estimate = c(coef(fit_trans)["trtTreatment"],
               coef(fit_marginal)["trtTreatment"],
               fixef(fit_mixed)["trtTreatment"]),
  SE = c(summary(fit_trans)$coefficients["trtTreatment", "Std. Error"],
         summary(fit_marginal)$coefficients["trtTreatment", "Std.err"],
         summary(fit_mixed)$coefficients["trtTreatment", "Std. Error"])
) |> mutate(across(where(is.numeric), ~round(., 3)))

# Sample 20 subjects
set.seed(667)
sample_ids <- sample(unique(dat$id), 20)

dat_sample <- dat |>
  filter(id %in% sample_ids)

ggplot(dat_sample, aes(time, y, group = id, color = trt)) +
  geom_line(alpha = 0.6) +
  geom_point(size = 1.5) +
  facet_wrap(~trt) +
  labs(title = "Individual Binary Trajectories",
       x = "Time", y = "Response (0/1)") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")