R/methods.R

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

## methods.R ---------------------------------------------------------------

#' @export
print.modMStates_fit <- function(x, digits = 4, ...) {
  cat("Continuous-time Markov multi-state model (modMStates)\n")
  if (!is.null(x$process)) cat("Structure      :", x$process, "\n")
  cat("Subjects       :", x$n_subjects,
      " Observations:", x$n_observations, "\n")
  cat("Free intensities:", x$npar,
      " log-likelihood:", format(x$loglik, digits = digits), "\n")
  cat("Converged      :", x$converged,
      paste0("(code ", x$convergence_code, ")"), "\n")
  if (!is.null(x$deathexact))
    cat("Exact entry states:", paste(x$deathexact, collapse = ", "), "\n")
  cat("\nTransition intensities (95% CI):\n")
  K <- nrow(x$qmatrix$estimates)
  lab <- colnames(x$qmatrix$estimates)
  for (r in seq_len(K)) for (s in seq_len(K)) {
    if (r == s || x$qmatrix$estimates[r, s] == 0) next
    cat(sprintf("  %-12s -> %-12s %8.4f  (%.4f, %.4f)  n = %d\n",
                lab[r], lab[s], x$qmatrix$estimates[r, s],
                x$qmatrix$ci.lower[r, s], x$qmatrix$ci.upper[r, s],
                x$counts[r, s]))
  }
  sparse <- which(x$counts > 0 & x$counts < 10 &
                    x$qmatrix$estimates > 0, arr.ind = TRUE)
  if (nrow(sparse) > 0)
    cat("  Note: fewer than ten observed transitions for ",
        paste(sprintf("%d->%d", sparse[, 1], sparse[, 2]), collapse = ", "),
        "; the corresponding intervals are unreliable.\n", sep = "")
  if (!is.null(x$sojourn)) {
    cat("\nMean sojourn time per visit to state:\n")
    print(round(as.data.frame(x$sojourn), digits))
  }
  if (!is.null(x$pmatrix)) {
    cat("\nTransition probabilities at t =", x$horizon, ":\n")
    print(round(x$pmatrix, digits))
  }
  invisible(x)
}

#' @export
summary.modMStates_fit <- function(object, ...) {
  structure(list(fit = object), class = "summary.modMStates_fit")
}

#' @export
print.summary.modMStates_fit <- function(x, ...) {
  print(x$fit, ...)
  cat("\nObserved transition counts:\n")
  print(x$fit$counts)
  invisible(x)
}

#' Monte Carlo evaluation of the estimator, with Monte Carlo standard errors
#'
#' Repeatedly simulates panel data from a known generator, fits the model, and
#' summarises bias, relative bias, empirical standard error, root mean squared
#' error and interval coverage. Every summary is reported with its Monte Carlo
#' standard error, so that differences smaller than the simulation noise are
#' not read as findings. Convergence is taken from the optimiser's own code
#' rather than from the absence of an error, since a run that fails to
#' converge returns a value rather than signalling a condition.
#'
#' @param process Character. One of \code{\link{ms_structures}()}.
#' @param n Integer vector of sample sizes.
#' @param B Integer. Replicates per cell.
#' @param t Numeric. Observation window.
#' @param Q Optional generator; defaults to the reference generator for
#'   \code{process}.
#' @param truth Optional vector of true intensities in row-major order of the
#'   permitted transitions; defaults to the values in \code{Q}.
#' @param horizon Numeric. Horizon for the transition probability matrix.
#' @param sim_args A list of further arguments passed to
#'   \code{\link{sim_mspdata}}, for example
#'   \code{list(schedule = "random", p_miss = 0.2)} to evaluate the estimator
#'   under an irregular visit schedule, or
#'   \code{list(sojourn = "weibull", shape = 1.5)} for the semi-Markov
#'   misspecification study.
#' @param seed Optional integer seed.
#' @param verbose Logical. Report progress by cell.
#'
#' @return A data frame with one row per parameter and sample size, holding
#'   the truth, the mean estimate, bias, relative bias, empirical standard
#'   error, RMSE and coverage, each with its Monte Carlo standard error, plus
#'   the number of converged replicates.
#'
#' @examples
#' \donttest{
#' ms_montecarlo("two_state", n = 100, B = 25, seed = 1)
#' }
#' @export
ms_montecarlo <- function(process, n = c(100, 300, 500), B = 1000, t = 10,
                          Q = NULL, truth = NULL, horizon = 5,
                          sim_args = list(), seed = NULL, verbose = TRUE) {

  if (!is.null(seed)) set.seed(seed)
  Qref <- .ms_resolve_Q(process, Q)
  idx <- which(ms_allowed(process) == 1, arr.ind = TRUE)
  idx <- idx[order(idx[, 1], idx[, 2]), , drop = FALSE]
  P <- nrow(idx)
  qtrue <- if (is.null(truth)) Qref[idx] else truth
  pname <- sprintf("q%d%d", idx[, 1], idx[, 2])

  res <- list()
  for (nn in n) {
    est <- lo <- hi <- matrix(NA_real_, B, P)
    ok <- logical(B)
    for (b in seq_len(B)) {
      dat <- do.call(sim_mspdata,
                     c(list(process = process, n = nn, t = t, Q = Q), sim_args))
      f <- tryCatch(fit_msm(dat, process = process, t = horizon),
                    error = function(e) NULL, warning = function(w) NULL)
      if (is.null(f) || !isTRUE(f$converged)) next
      ok[b] <- TRUE
      est[b, ] <- f$qmatrix$estimates[idx]
      lo[b, ]  <- f$qmatrix$ci.lower[idx]
      hi[b, ]  <- f$qmatrix$ci.upper[idx]
    }
    Bc <- sum(ok)
    if (Bc < 2L) {
      warning("Fewer than two converged replicates at n = ", nn, ".",
              call. = FALSE)
      next
    }
    e <- est[ok, , drop = FALSE]
    cov_ind <- sweep(lo[ok, , drop = FALSE], 2, qtrue, "<=") &
               sweep(hi[ok, , drop = FALSE], 2, qtrue, ">=")
    sqerr <- sweep(e, 2, qtrue, "-")^2

    mean_q <- colMeans(e)
    ese <- apply(e, 2, stats::sd)
    rmse <- sqrt(colMeans(sqerr))
    cp <- colMeans(cov_ind)

    res[[length(res) + 1L]] <- data.frame(
      process = process, n = nn, parameter = pname, true = qtrue,
      mean = mean_q,
      bias = mean_q - qtrue,
      mcse_bias = ese / sqrt(Bc),
      rbias_pct = 100 * (mean_q - qtrue) / qtrue,
      mcse_rbias_pct = 100 * (ese / sqrt(Bc)) / qtrue,
      ese = ese,
      rmse = rmse,
      mcse_rmse = apply(sqerr, 2, stats::sd) / (2 * rmse * sqrt(Bc)),
      coverage_pct = 100 * cp,
      mcse_coverage_pct = 100 * sqrt(cp * (1 - cp) / Bc),
      nonconvergence_pct = 100 * (1 - Bc / B),
      B_converged = Bc,
      row.names = NULL, stringsAsFactors = FALSE)
    if (verbose)
      message(sprintf("%s  n = %d: %d/%d converged", process, nn, Bc, B))
  }
  do.call(rbind, res)
}