R/simulate.R

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

## simulate.R -------------------------------------------------------------
## Trajectory simulation and panel observation.
##
## Design note. Trajectories are generated exactly, by drawing holding times
## and competing destinations, and the panel record is then formed by
## evaluating the trajectory at the observation times. The state recorded at
## observation time u is the state occupied at u, i.e. the destination of the
## last jump at or before u. Discretising the process directly onto the
## observation grid -- for instance by treating a jump inside an interval as
## having happened at the start or the end of that interval -- shifts the
## implied holding-time distribution by up to one inter-visit interval and
## biases every leaving rate. Panel construction is therefore kept strictly
## separate from trajectory generation, and the correctness of the pair is
## checked against the matrix exponential in tests/testthat/test-simulate.R.

.ms_traj <- function(Q, t_max, start = 1L, sojourn = c("exponential", "weibull"),
                     shape = 1) {
  sojourn <- match.arg(sojourn)
  K <- nrow(Q)
  s <- as.integer(start)
  tt <- 0
  times <- 0
  states <- s
  repeat {
    rate <- -Q[s, s]
    if (!is.finite(rate) || rate <= 0) break          # absorbing state
    wait <- if (sojourn == "exponential") {
      stats::rexp(1L, rate)
    } else {
      ## Weibull sojourn calibrated so that the mean holding time equals the
      ## Markov mean 1/rate at every shape. shape == 1 therefore reproduces
      ## the exponential case exactly, which makes it a built-in correctness
      ## check on the whole misspecification analysis.
      stats::rweibull(1L, shape = shape,
                      scale = (1 / rate) / gamma(1 + 1 / shape))
    }
    tt <- tt + wait
    if (tt > t_max) break
    p <- Q[s, ]
    p[s] <- 0
    s <- sample.int(K, 1L, prob = p / sum(p))
    times <- c(times, tt)
    states <- c(states, s)
  }
  list(times = times, states = states)
}

.ms_observe <- function(traj, obs_times)
  traj$states[findInterval(obs_times, traj$times)]

.ms_visit_times <- function(t_max, schedule, by, visit_rate, p_miss) {
  if (schedule == "regular") {
    v <- seq(0, t_max, by = by)
  } else {
    gaps <- stats::rexp(ceiling(3 * t_max * visit_rate) + 5L, visit_rate)
    v <- cumsum(gaps)
    v <- c(0, v[v < t_max])
  }
  if (p_miss > 0 && length(v) > 1L)
    v <- v[c(TRUE, stats::runif(length(v) - 1L) > p_miss)]
  v
}

#' Simulate panel multi-state data
#'
#' Generates exact continuous-time Markov trajectories under a built-in or
#' user-supplied generator and records the occupied state at a set of
#' observation times, producing a panel (interval-censored) dataset in long
#' format.
#'
#' @param process Character. One of \code{\link{ms_structures}()}. May be
#'   \code{NULL} if \code{Q} is supplied.
#' @param n Integer. Number of subjects.
#' @param t Numeric. Length of the observation window.
#' @param Q Optional user-supplied generator matrix. When both \code{process}
#'   and \code{Q} are given, \code{Q} is checked against the structural zeros
#'   of \code{process}.
#' @param schedule Either \code{"regular"} (visits on a common grid, the
#'   default) or \code{"random"} (subject-specific visit times from a Poisson
#'   process). Panel-data methods are motivated by irregular observation, so
#'   \code{"random"} is the setting in which the estimator should be checked
#'   before it is applied to real cohort data.
#' @param by Numeric. Spacing of the regular grid. Ignored when
#'   \code{schedule = "random"}.
#' @param visit_rate Numeric. Visit intensity for the random schedule.
#' @param p_miss Numeric in \code{[0, 1)}. Probability that any visit after
#'   baseline is missed, applied to either schedule.
#' @param start_state Integer state occupied at time zero, or a vector of
#'   probabilities over states from which the initial state is drawn.
#' @param exact_absorption Logical. If \code{TRUE}, the exact entry time into
#'   an absorbing state is appended as an extra record for subjects who reach
#'   one within the window. Real cohorts usually record death exactly even
#'   when everything else is panel-observed; set this to \code{TRUE} to
#'   reproduce that structure, and pass the absorbing state to the
#'   \code{deathexact} argument of \code{\link{fit_msm}}.
#' @param sojourn \code{"exponential"} for the Markov process, or
#'   \code{"weibull"} for a semi-Markov process with Weibull holding times.
#' @param shape Weibull shape parameter, used only when
#'   \code{sojourn = "weibull"}. The scale is calibrated so that the mean
#'   holding time in each state matches the Markov mean, so
#'   \code{shape = 1} recovers the Markov process exactly.
#'
#' @return A data frame with columns \code{subject}, \code{time} and
#'   \code{state}, one row per subject-visit, ordered by subject and time.
#'   The generating generator is attached as attribute \code{"Q"} and the
#'   simulation settings as attribute \code{"settings"}.
#'
#' @examples
#' dat <- sim_mspdata("illness_death_3state", n = 50, t = 10)
#' head(dat)
#' table(dat$state)
#'
#' ## Irregular visits with dropout
#' irr <- sim_mspdata("illness_death_3state", n = 50, t = 10,
#'                    schedule = "random", visit_rate = 1.2, p_miss = 0.2)
#' range(table(irr$subject))
#' @export
sim_mspdata <- function(process = NULL, n = 100L, t = 10,
                        Q = NULL,
                        schedule = c("regular", "random"),
                        by = 1, visit_rate = 1.2, p_miss = 0,
                        start_state = 1L,
                        exact_absorption = FALSE,
                        sojourn = c("exponential", "weibull"),
                        shape = 1) {

  schedule <- match.arg(schedule)
  sojourn  <- match.arg(sojourn)
  Q <- .ms_resolve_Q(process, Q)
  K <- nrow(Q)

  if (!is.numeric(n) || length(n) != 1L || n < 1)
    stop("'n' must be a single positive integer.", call. = FALSE)
  if (!is.numeric(t) || length(t) != 1L || t <= 0)
    stop("'t' must be a single positive number.", call. = FALSE)
  if (p_miss < 0 || p_miss >= 1)
    stop("'p_miss' must lie in [0, 1).", call. = FALSE)
  if (sojourn == "weibull" && (!is.numeric(shape) || shape <= 0))
    stop("'shape' must be a positive number.", call. = FALSE)

  absorbing <- which(diag(Q) == 0)
  n <- as.integer(n)

  draw_start <- if (length(start_state) == 1L) {
    function() as.integer(start_state)
  } else {
    if (length(start_state) != K)
      stop("'start_state' must be a single state or a vector of length ", K,
           ".", call. = FALSE)
    function() sample.int(K, 1L, prob = start_state)
  }

  out <- vector("list", n)
  for (i in seq_len(n)) {
    v <- .ms_visit_times(t, schedule, by, visit_rate, p_miss)
    tr <- .ms_traj(Q, t, draw_start(), sojourn, shape)
    di <- data.frame(subject = i, time = v, state = .ms_observe(tr, v))
    if (exact_absorption) {
      last <- length(tr$states)
      if (tr$states[last] %in% absorbing && tr$times[last] <= t) {
        di <- di[di$time < tr$times[last], , drop = FALSE]
        di <- rbind(di, data.frame(subject = i, time = tr$times[last],
                                   state = tr$states[last]))
      }
    }
    ## Follow-up stops at absorption: an absorbing state observed at several
    ## consecutive visits would otherwise be read as a sequence of
    ## absorbing-to-absorbing transitions.
    abs_hit <- which(di$state %in% absorbing)
    if (length(abs_hit)) di <- di[seq_len(abs_hit[1]), , drop = FALSE]
    out[[i]] <- di
  }

  dat <- do.call(rbind, out)
  dat <- dat[order(dat$subject, dat$time), ]
  rownames(dat) <- NULL
  attr(dat, "Q") <- Q
  attr(dat, "settings") <- list(process = process, n = n, t = t,
                                schedule = schedule, by = by,
                                visit_rate = visit_rate, p_miss = p_miss,
                                exact_absorption = exact_absorption,
                                sojourn = sojourn, shape = shape)
  dat
}

#' Simulate panel data from a semi-Markov process with Weibull sojourns
#'
#' Convenience wrapper on \code{\link{sim_mspdata}} for assessing what happens
#' when the exponential sojourn assumption implied by the Markov property is
#' false. Holding times in each transient state are Weibull with the supplied
#' shape and a scale calibrated so that the mean holding time is unchanged;
#' destinations are drawn from the same competing-risk probabilities as the
#' Markov process. Shape one reproduces the Markov process, so fitting a
#' Markov model to data from this function at \code{shape = 1} must recover
#' the generator, and any departure there indicates a fault in the pipeline
#' rather than a consequence of misspecification.
#'
#' @inheritParams sim_mspdata
#' @param shape Weibull shape. Values below one give a decreasing hazard with
#'   rapid early exits; values above one give an increasing hazard with
#'   delayed exits.
#'
#' @return As \code{\link{sim_mspdata}}.
#' @examples
#' dat <- sim_semimarkov("illness_death_3state", n = 50, t = 10, shape = 1.5)
#' head(dat)
#' @export
sim_semimarkov <- function(process = NULL, n = 100L, t = 10, shape = 1,
                           Q = NULL, schedule = c("regular", "random"),
                           by = 1, visit_rate = 1.2, p_miss = 0,
                           start_state = 1L) {
  sim_mspdata(process = process, n = n, t = t, Q = Q,
              schedule = match.arg(schedule), by = by,
              visit_rate = visit_rate, p_miss = p_miss,
              start_state = start_state,
              sojourn = "weibull", shape = shape)
}

#' Theoretical state occupation probabilities
#'
#' Row \code{r} of \code{expm(Q t)} for each requested time, computed with
#' \code{msm::MatrixExp}. Provided so that simulated data can be checked
#' against the process they are supposed to come from.
#'
#' @param Q A generator matrix, or \code{NULL} to use \code{process}.
#' @param times Numeric vector of times.
#' @param process Character. One of \code{\link{ms_structures}()}.
#' @param from Integer. Starting state.
#'
#' @return A matrix with one row per element of \code{times} and one column
#'   per state.
#' @examples
#' ms_occupancy(process = "illness_death_3state", times = 0:5)
#' @export
ms_occupancy <- function(Q = NULL, times, process = NULL, from = 1L) {
  Q <- .ms_resolve_Q(process, Q)
  out <- t(vapply(times, function(u) msm::MatrixExp(Q * u)[from, ],
                  numeric(nrow(Q))))
  dimnames(out) <- list(paste0("t=", times), colnames(Q))
  out
}