R/structures.R

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

## structures.R -----------------------------------------------------------
## Canonical process structures, generator construction, and validation.
## Every structure is defined once, here, and both the simulator and the
## fitting function read from this single definition. Nothing downstream
## constructs a generator matrix independently.

.ms_defs <- list(

  two_state = list(
    label  = "Two-state failure",
    states = c("Event-free", "Event"),
    rates  = list(c(1, 2, 0.20))
  ),

  recurrent_event = list(
    label  = "Recurrent event",
    states = c("Event-free", "Event"),
    rates  = list(c(1, 2, 0.40), c(2, 1, 0.30))
  ),

  illness_death_3state = list(
    label  = "Illness-death (three-state)",
    states = c("Healthy", "Ill", "Dead"),
    rates  = list(c(1, 2, 0.20), c(1, 3, 0.10), c(2, 3, 0.40))
  ),

  illness_death_4state = list(
    label  = "Illness-death (four-state)",
    states = c("Healthy", "Mild", "Severe", "Dead"),
    rates  = list(c(1, 2, 0.15), c(1, 4, 0.10), c(2, 3, 0.20),
                  c(2, 4, 0.15), c(3, 4, 0.50))
  ),

  reversible_illness_death = list(
    label  = "Reversible illness-death",
    states = c("Healthy", "Ill", "Dead"),
    rates  = list(c(1, 2, 0.20), c(2, 1, 0.15), c(1, 3, 0.05), c(2, 3, 0.10))
  ),

  competing_risks = list(
    label  = "Competing risks",
    states = c("At risk", "Cause 1", "Cause 2"),
    rates  = list(c(1, 2, 0.10), c(1, 3, 0.05))
  ),

  complex_hybrid = list(
    label  = "Complex hybrid",
    states = c("Healthy", "Ill", "Cause 1", "Cause 2"),
    rates  = list(c(1, 2, 0.15), c(1, 3, 0.10), c(1, 4, 0.05),
                  c(2, 1, 0.10), c(2, 3, 0.15), c(2, 4, 0.10))
  )
)

#' Names of the built-in process structures
#'
#' @return A character vector of the seven structure names accepted by the
#'   \code{process} argument of \code{\link{sim_mspdata}} and
#'   \code{\link{fit_msm}}.
#' @examples
#' ms_structures()
#' @export
ms_structures <- function() names(.ms_defs)

.ms_lookup <- function(process) {
  if (!is.character(process) || length(process) != 1L)
    stop("'process' must be a single character string; see ms_structures().",
         call. = FALSE)
  if (!process %in% names(.ms_defs))
    stop("Unknown process '", process, "'. Available structures: ",
         paste(names(.ms_defs), collapse = ", "), ".", call. = FALSE)
  .ms_defs[[process]]
}

#' Generator matrix for a built-in process structure
#'
#' Returns the structurally valid infinitesimal generator for one of the
#' seven canonical structures, with the reference intensities used
#' throughout the package documentation and simulation study.
#'
#' @param process Character. One of \code{\link{ms_structures}()}.
#' @param rates Optional named or unnamed numeric vector of replacement
#'   off-diagonal intensities, in the order in which the permitted
#'   transitions appear in \code{ms_allowed(process)} (row-major). Use this
#'   to keep the structure but change the parameter values.
#'
#' @return A \code{K} by \code{K} generator matrix with zero row sums and
#'   dimnames taken from the clinical state labels.
#' @examples
#' ms_generator("illness_death_3state")
#' ms_generator("two_state", rates = 0.35)
#' @export
ms_generator <- function(process, rates = NULL) {
  def <- .ms_lookup(process)
  K <- length(def$states)
  Q <- matrix(0, K, K, dimnames = list(def$states, def$states))
  for (r in def$rates) Q[r[1], r[2]] <- r[3]
  if (!is.null(rates)) {
    idx <- which(ms_allowed(process) == 1, arr.ind = TRUE)
    idx <- idx[order(idx[, 1], idx[, 2]), , drop = FALSE]
    if (length(rates) != nrow(idx))
      stop("'rates' has length ", length(rates), " but structure '", process,
           "' has ", nrow(idx), " free intensities.", call. = FALSE)
    Q[] <- 0
    Q[idx] <- rates
  }
  diag(Q) <- 0
  diag(Q) <- -rowSums(Q)
  Q
}

#' Permitted-transition indicator matrix for a built-in structure
#'
#' @param process Character. One of \code{\link{ms_structures}()}.
#' @return A \code{K} by \code{K} matrix of zeros and ones; a one in position
#'   \code{(r, s)} means the transition from state \code{r} to state \code{s}
#'   is permitted by the structure.
#' @examples
#' ms_allowed("competing_risks")
#' @export
ms_allowed <- function(process) {
  Q <- ms_generator(process)
  A <- matrix(0L, nrow(Q), ncol(Q), dimnames = dimnames(Q))
  A[Q > 0] <- 1L
  A
}

#' State labels for a built-in structure
#'
#' @param process Character. One of \code{\link{ms_structures}()}.
#' @return A character vector of clinical state labels.
#' @examples
#' ms_states("complex_hybrid")
#' @export
ms_states <- function(process) .ms_lookup(process)$states

#' Validate a generator matrix
#'
#' Checks that a user-supplied generator is usable: square, numeric, finite,
#' with non-negative off-diagonal elements and at least one transient state.
#' Row sums are repaired on the diagonal rather than being required to be
#' exactly zero on input, so a matrix of intensities with an arbitrary
#' diagonal is accepted. Unreachable transient states are reported.
#'
#' @param Q A square numeric matrix of transition intensities.
#' @param strict Logical. If \code{TRUE}, an unreachable transient state or a
#'   state with no outgoing and no incoming transitions raises an error rather
#'   than a warning.
#'
#' @return The validated generator, invisibly, with the diagonal set to minus
#'   the row sum of the off-diagonal elements.
#' @examples
#' Q <- rbind(c(0, 0.2, 0.1), c(0, 0, 0.4), c(0, 0, 0))
#' ms_validate_Q(Q)
#' @export
ms_validate_Q <- function(Q, strict = FALSE) {
  if (!is.matrix(Q) || !is.numeric(Q))
    stop("'Q' must be a numeric matrix.", call. = FALSE)
  if (nrow(Q) != ncol(Q))
    stop("'Q' must be square; got ", nrow(Q), " by ", ncol(Q), ".",
         call. = FALSE)
  if (nrow(Q) < 2L)
    stop("'Q' must have at least two states.", call. = FALSE)
  if (any(!is.finite(Q)))
    stop("'Q' contains non-finite values.", call. = FALSE)

  off <- Q
  diag(off) <- 0
  if (any(off < 0))
    stop("Off-diagonal elements of 'Q' must be non-negative; offending ",
         "positions: ",
         paste(apply(which(off < 0, arr.ind = TRUE), 1, paste,
                     collapse = "->"), collapse = ", "), ".", call. = FALSE)

  diag(Q) <- 0
  diag(Q) <- -rowSums(Q)

  if (all(rowSums(off) == 0))
    stop("'Q' has no permitted transitions.", call. = FALSE)

  ## Reachability from state 1 over the directed transition graph.
  K <- nrow(Q)
  reach <- rep(FALSE, K)
  reach[1] <- TRUE
  repeat {
    new <- reach | apply(off[reach, , drop = FALSE] > 0, 2, any)
    if (identical(new, reach)) break
    reach <- new
  }
  if (!all(reach)) {
    msg <- paste0("States unreachable from state 1: ",
                  paste(which(!reach), collapse = ", "),
                  ". No transitions into them can be estimated.")
    if (strict) stop(msg, call. = FALSE) else warning(msg, call. = FALSE)
  }
  invisible(Q)
}

.ms_resolve_Q <- function(process = NULL, Q = NULL) {
  if (is.null(process) && is.null(Q))
    stop("Supply either 'process' or 'Q'.", call. = FALSE)
  if (!is.null(Q)) {
    Q <- ms_validate_Q(Q)
    if (!is.null(process)) {
      A <- ms_allowed(process)
      if (!identical(dim(A), dim(Q)))
        stop("Supplied 'Q' is ", nrow(Q), " by ", ncol(Q), " but structure '",
             process, "' has ", nrow(A), " states.", call. = FALSE)
      off <- Q; diag(off) <- 0
      bad <- which(off > 0 & A == 0, arr.ind = TRUE)
      if (nrow(bad) > 0)
        stop("Supplied 'Q' has positive intensities on transitions the '",
             process, "' structure forbids: ",
             paste(apply(bad, 1, paste, collapse = "->"), collapse = ", "),
             ".", call. = FALSE)
    }
    return(Q)
  }
  ms_generator(process)
}