inst/scripts/run_paper_analyses.R

In modMStates: Simulation and Estimation of Continuous-Time Multi-State Markov Models for Panel Data

## run_paper_analyses.R ----------------------------------------------------
## Every number reported in the manuscript is produced by this script.
## Run order matters: the pipeline check comes first, because the Monte Carlo
## results are only meaningful once the simulator has been shown to generate
## the process it claims.
##
##   1  Pipeline check   simulated occupancy against the matrix exponential
##   2  Validation       fit_msm() against a hand-written msm() call
##   3  Monte Carlo      bias, RMSE, coverage, with Monte Carlo standard errors
##   4  Irregular visits the same study under a random observation schedule
##   5  Misspecification semi-Markov Weibull sojourns
##   6  Application      pbcseq

library(modMStates)
library(msm)
library(survival)

B_MAIN <- as.integer(Sys.getenv("MODM_B", "1000"))
OUT <- Sys.getenv("MODM_OUT", file.path(tempdir(), "modMStates-results"))
dir.create(OUT, showWarnings = FALSE, recursive = TRUE)
set.seed(2026)

occupancy <- function(dat, times, K) {
  n <- length(unique(dat$subject))
  t(vapply(times, function(u) {
    rows <- dat[dat$time <= u, ]
    last <- rows$state[!duplicated(rows$subject, fromLast = TRUE)]
    tabulate(last, nbins = K) / n
  }, numeric(K)))
}

## -- 1. Pipeline check ----------------------------------------------------
message("1. pipeline check")
n_chk <- 20000
times <- 0:10
pipeline <- do.call(rbind, lapply(ms_structures(), function(p) {
  K <- length(ms_states(p))
  dat <- sim_mspdata(p, n = n_chk, t = 10)
  emp <- occupancy(dat, times, K)
  theo <- ms_occupancy(process = p, times = times)
  z <- (emp - theo) / pmax(sqrt(theo * (1 - theo) / n_chk), 1e-8)
  data.frame(process = p,
             max_abs_diff = max(abs(emp - theo)),
             max_abs_z = max(abs(z[is.finite(z)])))
}))
print(pipeline)
write.csv(pipeline, file.path(OUT, "01_pipeline_check.csv"), row.names = FALSE)

## -- 2. Validation against msm -------------------------------------------
message("2. validation against msm")
validation <- do.call(rbind, lapply(ms_structures(), function(p) {
  dat <- sim_mspdata(p, n = 500, t = 10)
  wrapped <- fit_msm(dat, p, t = 5)
  qinit <- crudeinits.msm(state ~ time, subject = subject, data = dat,
                          qmatrix = ms_allowed(p) * 1)
  direct <- msm(state ~ time, subject = subject, data = dat, qmatrix = qinit)
  data.frame(process = p,
             free_intensities = sum(ms_allowed(p)),
             max_abs_diff_Q = max(abs(qmatrix.msm(direct, ci = "none") -
                                        wrapped$qmatrix$estimates)),
             abs_diff_loglik = abs(as.numeric(logLik(direct)) - wrapped$loglik))
}))
print(validation)
write.csv(validation, file.path(OUT, "02_validation_vs_msm.csv"), row.names = FALSE)

## -- 3. Monte Carlo -------------------------------------------------------
message("3. Monte Carlo, B = ", B_MAIN)
mc <- do.call(rbind, lapply(ms_structures(), function(p)
  ms_montecarlo(p, n = c(100, 300, 500), B = B_MAIN, t = 10,
                horizon = 5, seed = 2026)))
write.csv(mc, file.path(OUT, "03_montecarlo.csv"), row.names = FALSE)

## -- 4. Irregular observation --------------------------------------------
message("4. irregular observation schedule")
mc_irr <- do.call(rbind, lapply(ms_structures(), function(p)
  ms_montecarlo(p, n = c(100, 300, 500), B = B_MAIN, t = 10, horizon = 5,
                sim_args = list(schedule = "random", visit_rate = 1.2,
                                p_miss = 0.20),
                seed = 4026)))
write.csv(mc_irr, file.path(OUT, "04_montecarlo_irregular.csv"), row.names = FALSE)

## -- 5. Misspecification: semi-Markov Weibull sojourns -------------------
## shape = 1 is the Markov case and must reproduce the corresponding rows of
## the main study; it is retained as a correctness check, not as a result.
message("5. Weibull sojourn misspecification")
weib <- do.call(rbind, lapply(c(0.7, 1.0, 1.5), function(k) {
  out <- ms_montecarlo("illness_death_3state", n = 300, B = B_MAIN, t = 10,
                       horizon = 5,
                       sim_args = list(sojourn = "weibull", shape = k),
                       seed = 6026 + round(100 * k))
  out$shape <- k
  out
}))
write.csv(weib, file.path(OUT, "05_weibull_misspecification.csv"), row.names = FALSE)

## -- 6. Clinical application ---------------------------------------------
## pbc is one row per patient and cannot support a multi-state analysis.
## pbcseq is its longitudinal companion: repeated clinic visits at irregular
## times with serum bilirubin measured at each. States are defined from that
## measurement; no transition times are imputed.
message("6. pbcseq application")
data(pbc, package = "survival")   # loads both pbc and pbcseq

BILI_CUT <- 3.0                   # mg/dL; conventional marker of advanced disease

visits <- data.frame(subject = pbcseq$id,
                     time = pbcseq$day / 365.25,
                     state = ifelse(pbcseq$bili >= BILI_CUT, 2L, 1L))
visits <- visits[!is.na(visits$state), ]

deaths <- pbc[pbc$status == 2, c("id", "time")]
deaths <- data.frame(subject = deaths$id, time = deaths$time / 365.25,
                     state = 3L)
deaths <- deaths[deaths$subject %in% visits$subject, ]

ms <- rbind(visits, deaths)
ms <- ms[order(ms$subject, ms$time), ]
ms <- ms[!duplicated(ms[, c("subject", "time")]), ]
## Follow-up ends at death.
ms <- do.call(rbind, lapply(split(ms, ms$subject), function(x) {
  d <- which(x$state == 3L)
  if (length(d)) x[seq_len(d[1]), , drop = FALSE] else x
}))
rownames(ms) <- NULL

app_summary <- c(subjects = length(unique(ms$subject)),
                 observations = nrow(ms),
                 median_visits = median(table(ms$subject)),
                 deaths = sum(ms$state == 3L))
print(app_summary)

## Bilirubin can fall back below threshold, so recovery is observable and an
## irreversible structure would be misspecified by construction.
fit_panel <- fit_msm(ms, "reversible_illness_death", t = 5)
fit_exact <- fit_msm(ms, "reversible_illness_death", t = 5, deathexact = 3)
print(fit_exact)

qp <- fit_panel$qmatrix$estimates
qe <- fit_exact$qmatrix$estimates
idx <- which(ms_allowed("reversible_illness_death") == 1, arr.ind = TRUE)
idx <- idx[order(idx[, 1], idx[, 2]), , drop = FALSE]
app <- data.frame(transition = sprintf("%d->%d", idx[, 1], idx[, 2]),
                  panel = qp[idx],
                  exact_death = qe[idx],
                  ci_lower = fit_exact$qmatrix$ci.lower[idx],
                  ci_upper = fit_exact$qmatrix$ci.upper[idx],
                  count = fit_exact$counts[idx])
print(app)
write.csv(app, file.path(OUT, "06_pbcseq_intensities.csv"), row.names = FALSE)
write.csv(as.data.frame(fit_exact$sojourn),
          file.path(OUT, "06_pbcseq_sojourn.csv"))
write.csv(round(fit_exact$pmatrix, 4),
          file.path(OUT, "06_pbcseq_pmatrix5.csv"))
write.csv(fit_exact$counts, file.path(OUT, "06_pbcseq_counts.csv"))

## Markov assumption: observed against expected state prevalence.
prev <- prevalence.msm(fit_exact$fit, times = seq(0, 10, by = 1))
capture.output(print(prev), file = file.path(OUT, "06_pbcseq_prevalence.txt"))

sessionInfo()