tests/testthat/test-fit.R

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

test_that("fit_msm agrees numerically with a direct msm call", {
  skip_on_cran()
  skip_if_not_installed("msm")
  set.seed(21)
  for (p in ms_structures()) {
    dat <- sim_mspdata(p, n = 400, t = 10)
    wrapped <- fit_msm(dat, p, t = 5)
    qinit <- msm::crudeinits.msm(state ~ time, subject = subject, data = dat,
                                 qmatrix = ms_allowed(p) * 1)
    direct <- msm::msm(state ~ time, subject = subject, data = dat,
                       qmatrix = qinit)
    expect_lt(max(abs(msm::qmatrix.msm(direct, ci = "none") -
                        wrapped$qmatrix$estimates)), 1e-6)
    expect_equal(as.numeric(logLik(direct)), wrapped$loglik, tolerance = 1e-8)
  }
})

test_that("the estimator recovers the generating intensities at large n", {
  skip_on_cran()
  set.seed(22)
  dat <- sim_mspdata("illness_death_3state", n = 3000, t = 10)
  f <- fit_msm(dat, "illness_death_3state", t = 5)
  expect_true(f$converged)
  est <- c(f$qmatrix$estimates[1, 2], f$qmatrix$estimates[1, 3],
           f$qmatrix$estimates[2, 3])
  expect_equal(est, c(0.20, 0.10, 0.40), tolerance = 0.06)
})

test_that("output objects have the documented structure", {
  skip_on_cran()
  set.seed(23)
  dat <- sim_mspdata("illness_death_3state", n = 300, t = 10)
  f <- fit_msm(dat, "illness_death_3state", t = 5)
  expect_s3_class(f, "modMStates_fit")
  expect_named(f$qmatrix, c("estimates", "se", "ci.lower", "ci.upper"))
  expect_equal(dim(f$pmatrix), c(3L, 3L))
  expect_equal(unname(rowSums(f$pmatrix)), rep(1, 3), tolerance = 1e-8)
  expect_equal(dim(f$counts), c(3L, 3L))
  expect_true(is.finite(f$loglik))
  expect_true(all(f$qmatrix$ci.lower[1, 2:3] <= f$qmatrix$estimates[1, 2:3]))
  expect_true(all(f$qmatrix$ci.upper[1, 2:3] >= f$qmatrix$estimates[1, 2:3]))
})

test_that("the horizon argument changes pmatrix but not the estimates", {
  skip_on_cran()
  set.seed(24)
  dat <- sim_mspdata("illness_death_3state", n = 300, t = 10)
  f5 <- fit_msm(dat, "illness_death_3state", t = 5)
  f2 <- fit_msm(dat, "illness_death_3state", t = 2)
  expect_equal(f5$qmatrix$estimates, f2$qmatrix$estimates)
  expect_false(isTRUE(all.equal(f5$pmatrix, f2$pmatrix)))
})

test_that("input problems are caught before the optimiser is called", {
  set.seed(25)
  dat <- sim_mspdata("reversible_illness_death", n = 120, t = 10)
  ## Recovery is present in these data but forbidden by the irreversible model.
  expect_error(fit_msm(dat, "illness_death_3state"), "cannot produce")
  expect_error(fit_msm(dat, "illness_death_3state", state = "nope"),
               "not found")
  bad <- dat; bad$state[1] <- 9L
  expect_error(fit_msm(bad, "reversible_illness_death"), "consecutive integers")
  dup <- rbind(dat[1, ], dat)
  expect_error(fit_msm(dup, "reversible_illness_death"), "Duplicated")
  expect_error(fit_msm(as.list(dat), "reversible_illness_death"), "data frame")
})

test_that("deathexact changes the exit intensity", {
  skip_on_cran()
  set.seed(26)
  dat <- sim_mspdata("illness_death_3state", n = 600, t = 10,
                     exact_absorption = TRUE)
  panel_only <- fit_msm(dat, "illness_death_3state", t = 5)
  exact <- fit_msm(dat, "illness_death_3state", t = 5, deathexact = 3)
  expect_true(exact$converged)
  expect_false(isTRUE(all.equal(panel_only$qmatrix$estimates,
                                exact$qmatrix$estimates, tolerance = 1e-4)))
})

test_that("ms_montecarlo returns MCSEs and counts convergence", {
  skip_on_cran()
  set.seed(27)
  mc <- ms_montecarlo("two_state", n = 150, B = 15, seed = 5, verbose = FALSE)
  expect_true(all(c("mcse_bias", "mcse_rmse", "mcse_coverage_pct",
                    "nonconvergence_pct", "B_converged") %in% names(mc)))
  expect_equal(nrow(mc), 1L)
  expect_true(mc$B_converged <= 15L)
  expect_true(mc$mcse_bias > 0)
})

test_that("deprecated aliases still work but warn", {
  skip_on_cran()
  set.seed(28)
  expect_warning(d <- sim.mspdata("two_state", n = 20, t = 5), "deprecated")
  expect_warning(fit.msm(d, "two_state", t = 2), "deprecated")
})

test_that("a likelihood overflow is recovered by rescaling", {
  skip_on_cran()
  set.seed(101)
  ## At this size msm's objective overflows on its natural scale for the
  ## reversible structure; the fit must still succeed and recover the truth.
  dat <- sim_mspdata("reversible_illness_death", n = 3000, t = 10)
  f <- suppressWarnings(fit_msm(dat, "reversible_illness_death", t = 5))
  expect_true(f$converged)
  est <- c(f$qmatrix$estimates[1, 2], f$qmatrix$estimates[2, 1],
           f$qmatrix$estimates[1, 3], f$qmatrix$estimates[2, 3])
  expect_equal(est, c(0.20, 0.15, 0.05, 0.10), tolerance = 0.05)
})