R/fit.R

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

## fit.R -------------------------------------------------------------------
## Single-call fitting interface.
##
## Naming note. Earlier versions of this package exported the fitting
## function as fit.msm(). That name is read by R's S3 dispatch as the "msm"
## method of a generic called fit, and msm::msm() returns objects of class
## "msm", so the two collide whenever both packages are attached. The
## function is now fit_msm(); fit.msm() remains as a deprecated alias.

.ms_reachable <- function(A) {
  ## Transitive closure of the directed transition graph, with the diagonal
  ## set to TRUE so that staying in the same state is always admissible.
  K <- nrow(A)
  R <- A > 0
  for (k in seq_len(K)) R <- R | (R[, k, drop = FALSE] %*% R[k, , drop = FALSE] > 0)
  diag(R) <- TRUE
  R
}

.ms_counts_matrix <- function(state, subject, K, labels) {
  tab <- msm::statetable.msm(state = state, subject = subject)
  out <- matrix(0L, K, K, dimnames = list(labels, labels))
  rn <- as.integer(rownames(tab))
  cn <- as.integer(colnames(tab))
  for (i in seq_along(rn)) for (j in seq_along(cn))
    out[rn[i], cn[j]] <- as.integer(tab[i, j])
  out
}

#' Fit a continuous-time Markov multi-state model to panel data
#'
#' Fits the continuous-time Markov model implied by one of the seven built-in
#' process structures, or by a user-supplied generator, by maximum likelihood,
#' and returns all primary inferential summaries in a single object.
#' Likelihood evaluation is delegated to \code{\link[msm]{msm}}; this function
#' supplies the structurally valid generator, crude starting values, input
#' validation and the derived summaries.
#'
#' @param data A data frame in long format, one row per subject-visit.
#' @param process Character. One of \code{\link{ms_structures}()}. May be
#'   \code{NULL} if \code{Q} is supplied.
#' @param state,subject,time Character. Column names in \code{data}.
#' @param t Numeric. Horizon at which the transition probability matrix is
#'   evaluated. This affects the reported \code{pmatrix} only; it has no
#'   effect on estimation and is unrelated to the observation window.
#' @param Q Optional user-supplied generator giving the permitted transitions
#'   and, if \code{inits = "user"}, the starting values.
#' @param deathexact Integer vector of states whose entry times are recorded
#'   exactly rather than being interval-censored, passed through to
#'   \code{\link[msm]{msm}}. Cohorts that panel-observe clinical state but
#'   record death exactly need this; omitting it treats the death time as
#'   known only to lie in the last interval and biases the exit intensities.
#' @param inits Either \code{"crude"} (the default; starting values from
#'   \code{\link[msm]{crudeinits.msm}}) or \code{"user"} (use \code{Q}).
#' @param ci Interval type passed to \code{\link[msm]{qmatrix.msm}}. The
#'   default \code{"normal"} gives Wald intervals on the log intensity scale,
#'   back-transformed, so the reported intervals are asymmetric on the
#'   intensity scale.
#' @param ... Further arguments passed to \code{\link[msm]{msm}}.
#'
#' @return An object of class \code{"modMStates_fit"}: a list with elements
#'   \code{qmatrix} (a list of \code{estimates}, \code{ci.lower},
#'   \code{ci.upper}), \code{sojourn} (a data frame of mean sojourn times with
#'   intervals), \code{pmatrix}, \code{counts}, \code{converged},
#'   \code{convergence_code}, \code{loglik}, \code{npar}, \code{n_subjects},
#'   \code{n_observations}, \code{process}, \code{horizon} and \code{fit}, the
#'   underlying \code{msm} object.
#'
#' @details
#' Mean sojourn times are per visit to a state, not total time occupied. In
#' structures where a state can be re-entered -- the recurrent event,
#' reversible illness-death and complex hybrid structures -- the total time a
#' subject spends in that state over follow-up is larger than the reported
#' sojourn time by a factor equal to the expected number of visits.
#'
#' @examples
#' dat <- sim_mspdata("illness_death_3state", n = 200, t = 10)
#' fit <- fit_msm(dat, "illness_death_3state", t = 5)
#' fit
#' fit$counts
#' @export
fit_msm <- function(data, process = NULL, state = "state",
                    subject = "subject", time = "time", t = NULL,
                    Q = NULL, deathexact = NULL, inits = c("crude", "user"),
                    ci = c("normal", "none"), ...) {

  inits <- match.arg(inits)
  ci    <- match.arg(ci)

  if (!is.data.frame(data))
    stop("'data' must be a data frame.", call. = FALSE)
  miss <- setdiff(c(state, subject, time), names(data))
  if (length(miss))
    stop("Column(s) not found in 'data': ", paste(miss, collapse = ", "),
         ".", call. = FALSE)

  st <- data[[state]]
  sb <- data[[subject]]
  tm <- data[[time]]

  if (!is.numeric(tm))
    stop("Column '", time, "' must be numeric.", call. = FALSE)
  if (any(!is.finite(tm)))
    stop("Column '", time, "' contains missing or non-finite values.",
         call. = FALSE)
  if (any(is.na(st)))
    stop("Column '", state, "' contains missing values.", call. = FALSE)

  st <- as.integer(st)
  Qstruct <- .ms_resolve_Q(process, Q)
  K <- nrow(Qstruct)

  if (min(st) < 1L || max(st) > K)
    stop("States must be coded as consecutive integers 1..", K,
         "; observed range ", min(st), "..", max(st), ".", call. = FALSE)

  ord <- order(sb, tm)
  if (!identical(ord, seq_along(tm))) {
    data <- data[ord, , drop = FALSE]
    st <- st[ord]; sb <- sb[ord]; tm <- tm[ord]
  }
  dup <- duplicated(data.frame(sb, tm))
  if (any(dup))
    stop("Duplicated observation times within subject(s): ",
         paste(utils::head(unique(sb[dup]), 5), collapse = ", "),
         if (length(unique(sb[dup])) > 5) ", ...", ".", call. = FALSE)

  ## Observed pairs of consecutive states that the assumed structure cannot
  ## produce. Under panel observation a subject may pass through intermediate
  ## states between visits, so an observed pair (r, s) is legitimate whenever
  ## s is reachable from r in the transition graph; only unreachable pairs
  ## indicate the wrong structure. Left unchecked these surface as an opaque
  ## optimiser failure inside msm.
  counts <- .ms_counts_matrix(st, sb, K, colnames(Qstruct))
  reachable <- .ms_reachable(Qstruct > 0)
  bad <- which(counts > 0 & !reachable, arr.ind = TRUE)
  if (nrow(bad) > 0)
    stop("The data contain state changes the assumed structure cannot ",
         "produce: ",
         paste(sprintf("%d->%d (n=%d)", bad[, 1], bad[, 2],
                       counts[bad]), collapse = ", "),
         ". Choose a structure that permits them, or supply 'Q'.",
         call. = FALSE)

  qinit <- if (inits == "user") {
    Qstruct
  } else {
    msm::crudeinits.msm(st ~ tm, subject = sb, qmatrix = (Qstruct > 0) * 1)
  }

  args <- list(formula = st ~ tm, subject = sb, qmatrix = qinit)
  if (!is.null(deathexact)) args$deathexact <- deathexact
  dots <- list(...)

  ## The optimiser minimises -2 log L on its natural scale, which grows with
  ## the number of observed intervals. Above roughly a thousand subjects the
  ## objective can overflow during the search and the fit fails with
  ## "numerical overflow in calculating likelihood", even though the model is
  ## perfectly well identified. Rescaling the objective fixes this and leaves
  ## the maximum unchanged. The rescaling is applied only after an unscaled
  ## attempt has failed, so a fit that succeeds without it is bit-for-bit the
  ## same as the equivalent hand-written msm() call.
  rescaled <- FALSE
  fit <- tryCatch(do.call(msm::msm, c(args, dots)),
                  error = function(e) e)
  if (inherits(fit, "error")) {
    if (!grepl("overflow|non-finite|NA/NaN", conditionMessage(fit)) ||
        !is.null(dots$control$fnscale))
      stop(fit)
    scale <- max(100, length(st) - length(unique(sb)))
    ctrl <- dots$control
    ctrl$fnscale <- scale
    ctrl$maxit <- if (is.null(ctrl$maxit)) 10000L else ctrl$maxit
    dots$control <- ctrl
    fit <- tryCatch(do.call(msm::msm, c(args, dots)),
                    error = function(e)
                      stop("Optimisation failed both unscaled and with ",
                           "fnscale = ", scale, ". Last message: ",
                           conditionMessage(e), call. = FALSE))
    rescaled <- TRUE
    warning("The likelihood overflowed on its natural scale; refitted with ",
            "control$fnscale = ", scale, ". The maximum is unchanged.",
            call. = FALSE)
  }

  lab <- colnames(Qstruct)
  as_mat <- function(z) matrix(as.numeric(z), K, K, dimnames = list(lab, lab))

  qm <- if (ci == "none") {
    msm::qmatrix.msm(fit, ci = "none")
  } else {
    ## The simulation-based normal interval is on the log intensity scale and
    ## is therefore asymmetric on the intensity scale. It fails for models
    ## with a single free parameter, where msm's covariance is a scalar; the
    ## delta-method interval is used there instead.
    tryCatch(msm::qmatrix.msm(fit, ci = "normal"),
             error = function(e) msm::qmatrix.msm(fit, ci = "delta"))
  }

  if (is.list(qm)) {
    est <- as_mat(qm$estimates)
    L <- as_mat(qm$L)
    U <- as_mat(qm$U)
    SE <- if (is.null(qm$SE)) as_mat(rep(NA_real_, K * K)) else as_mat(qm$SE)
  } else {
    est <- as_mat(qm)
    L <- U <- SE <- as_mat(rep(NA_real_, K * K))
  }
  ci_type <- if (ci == "none") "none" else
    if (identical(attr(qm, "ci"), "delta")) "delta" else "normal (log scale)"

  soj <- try(msm::sojourn.msm(fit), silent = TRUE)
  if (inherits(soj, "try-error")) soj <- NULL

  horizon <- if (is.null(t)) NA_real_ else t
  pmat <- if (is.null(t)) NULL else {
    pm <- msm::pmatrix.msm(fit, t = t)
    matrix(as.numeric(pm), K, K, dimnames = dimnames(est))
  }

  code <- tryCatch(as.integer(fit$opt$convergence), error = function(e) NA_integer_)

  out <- list(
    qmatrix = list(estimates = est, se = SE, ci.lower = L, ci.upper = U),
    ci_type = ci_type,
    sojourn = soj,
    pmatrix = pmat,
    counts = counts,
    converged = isTRUE(code == 0L),
    convergence_code = code,
    rescaled = rescaled,
    loglik = as.numeric(stats::logLik(fit)),
    npar = sum(Qstruct > 0),
    n_subjects = length(unique(sb)),
    n_observations = length(st),
    process = process,
    horizon = horizon,
    deathexact = deathexact,
    fit = fit
  )
  class(out) <- "modMStates_fit"
  out
}

#' @rdname fit_msm
#' @param ... Passed to \code{fit_msm}.
#' @export
fit.msm <- function(...) {
  .Deprecated("fit_msm",
              msg = paste("fit.msm() is deprecated: the name collides with S3",
                          "dispatch on objects of class 'msm'. Use fit_msm()."))
  fit_msm(...)
}

#' @rdname sim_mspdata
#' @param ... Passed to \code{sim_mspdata}.
#' @export
sim.mspdata <- function(...) {
  .Deprecated("sim_mspdata")
  sim_mspdata(...)
}