tests/testthat/test-simulate.R

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

## The tests below are the ones that would have caught a panel-construction
## fault: they compare simulated occupancy against the matrix exponential and
## against a closed form, rather than only checking that the code runs.

## Occupancy is computed by carrying the last recorded state forward, because
## follow-up stops at absorption and absorbed subjects therefore have no rows
## at later visit times.
.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)))
}

test_that("panel observation reproduces the matrix exponential", {
  skip_on_cran()
  set.seed(11)
  n <- 6000
  dat <- sim_mspdata("illness_death_3state", n = n, t = 10)
  Q <- ms_generator("illness_death_3state")
  times <- 0:10
  emp <- .occupancy(dat, times, 3)
  theo <- ms_occupancy(process = "illness_death_3state", times = times)
  ## Three binomial standard errors, with a floor for the near-degenerate cells.
  tol <- 3 * sqrt(0.25 / n) + 0.005
  expect_lt(max(abs(emp - theo)), tol)
})

test_that("two-state failure matches the closed form and is not time-shifted", {
  skip_on_cran()
  set.seed(12)
  n <- 8000
  dat <- sim_mspdata("two_state", n = n, t = 10)
  times <- 0:10
  emp <- .occupancy(dat, times, 2)[, 2]
  theo <- 1 - exp(-0.20 * times)
  shift_down <- 1 - exp(-0.20 * pmax(times - 1, 0))
  shift_up <- 1 - exp(-0.20 * (times + 1))
  expect_lt(max(abs(emp - theo)), 3 * sqrt(0.25 / n) + 0.005)
  ## The correct curve must fit strictly better than either one-interval shift.
  expect_lt(sum((emp - theo)^2), sum((emp - shift_down)^2))
  expect_lt(sum((emp - theo)^2), sum((emp - shift_up)^2))
})

test_that("Weibull shape one reproduces the Markov process exactly", {
  skip_on_cran()
  set.seed(13)
  n <- 6000
  dat <- sim_semimarkov("illness_death_3state", n = n, t = 10, shape = 1)
  times <- 0:10
  emp <- .occupancy(dat, times, 3)
  theo <- ms_occupancy(process = "illness_death_3state", times = times)
  expect_lt(max(abs(emp - theo)), 3 * sqrt(0.25 / n) + 0.005)
})

test_that("irregular schedules produce subject-specific visit times", {
  set.seed(14)
  dat <- sim_mspdata("illness_death_3state", n = 40, t = 10,
                     schedule = "random", visit_rate = 1.2, p_miss = 0.2)
  n_visits <- table(dat$subject)
  expect_gt(stats::sd(as.numeric(n_visits)), 0)
  expect_true(all(tapply(dat$time, dat$subject, min) == 0))
  expect_true(all(dat$time <= 10))
})

test_that("simulated data have the documented shape", {
  set.seed(15)
  dat <- sim_mspdata("competing_risks", n = 25, t = 6)
  expect_named(dat, c("subject", "time", "state"))
  expect_equal(length(unique(dat$subject)), 25L)
  expect_true(all(dat$state %in% 1:3))
  expect_false(is.null(attr(dat, "Q")))
})

test_that("exact absorption times are appended when requested", {
  set.seed(16)
  dat <- sim_mspdata("two_state", n = 400, t = 10, exact_absorption = TRUE)
  non_integer <- dat$time[dat$time %% 1 != 0]
  expect_gt(length(non_integer), 0)
  expect_true(all(dat$state[dat$time %% 1 != 0] == 2))
})