Introduction to RFmstate

In RFmstate: Random Forest-Based Multistate Survival Analysis

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

# One configuration is used throughout the vignette. These values keep the
# source build practical; they are not inferential adequacy recommendations.
canonical_seed <- 42L
canonical_n <- 300L
canonical_trees <- 200L
canonical_min_events <- 3L
canonical_sparse_warning <- 20L
canonical_covariates <- c("age", "sex", "BMI", "treatment")

Overview

RFmstate fits clock-reset cause-specific random survival forests for acyclic, non-recurrent multistate processes. For each transient state, competing exits are modeled by separate forests. Patient/profile entry-conditioned state probabilities are assembled from predicted cumulative hazards by semi-Markov entry-mass and sojourn convolution. The package also provides calendar-time Aalen-Johansen point estimates as a covariate-free descriptive baseline. The supported one-row-per-subject contract uses a common initial state, one recorded entry per state, baseline covariates, right censoring, and competing exits; it does not support left truncation, recurrent visits, directed cycles, or time-dependent covariates.

The package provides:

• State space and transition structure definition
• Wide-to-long data conversion for counting-process format
• Cause-specific random forest fitting per origin state
• Entry-conditioned state probabilities via clock-reset semi-Markov convolution
• Aalen-Johansen calendar-time point estimation (covariate-free baseline)
• Per-transition feature importance
• Genuine edge OOB concordance and patient-level cross-validated IPCW scoring
• Comprehensive visualizations

Quick Start

1. Define the Multistate Structure

library(RFmstate)

# Use the built-in clinical trial structure
ms <- clinical_states()
print(ms)

Or define another supported single-root, acyclic, non-recurrent structure:

# A simple 3-state illness-death model
ms_simple <- define_multistate(
  state_names = c("Healthy", "Sick", "Dead"),
  absorbing = "Dead",
  transitions = list(
    Healthy = c("Sick", "Dead"),
    Sick = c("Dead")
  )
)

# A 4-state model with recovery
ms_recovery <- define_multistate(
  state_names = c("Healthy", "Sick", "Recovered", "Dead"),
  absorbing = "Dead",
  transitions = list(
    Healthy = c("Sick", "Dead"),
    Sick = c("Recovered", "Dead"),
    Recovered = c("Dead")
  )
)

The same workflow applies to a validated DAG with one common initial state and at least one absorbing state. Cycles and recurrent visits are rejected.

2. Simulate Data

dat <- sim_clinical_data(
  n = canonical_n, structure = ms, seed = canonical_seed
)
head(dat)

3. Prepare Multistate Data

Convert wide-format data to long format:

msdata <- prepare_data(
  data = dat, id = "ID", structure = ms,
  time_map = list(
    Responded = "time_Responded",
    Unresponded = "time_Unresponded",
    Stabilized = "time_Stabilized",
    Progressed = "time_Progressed",
    Death = "time_Death"
  ),
  censor_col = "time_censored",
  covariates = canonical_covariates
)
print(msdata)
head(msdata)

print(msdata) is a concise validation summary in the state-definition order; head(msdata) displays the first six ordinary data rows.

4. Aalen-Johansen Nonparametric Baseline

Compute the covariate-free calendar-time point-estimate benchmark:

aj <- aalen_johansen(msdata)
print(aj)

plot(aj, type = "state_occupation")

The figure is produced by the immediately preceding plot() call and shows occupation from the recorded common baseline. It has no confidence band.

plot(aj, type = "cumulative_hazard")

This second figure shows Nelson--Aalen cumulative cause-specific hazards. Requested destination states can be selected with states =.

5. Fit Random Forest Model

fit <- rfmstate(
  msdata,
  num.trees = canonical_trees,
  min_events = canonical_min_events,
  sparse_warning = canonical_sparse_warning,
  seed = canonical_seed
)
print(fit)

No covariate vector is repeated here: rfmstate(covariates = NULL) uses the explicit predictor contract stored by prepare_data(). An explicit vector may select a nonempty subset of that contract, but it cannot add structural, outcome-time, censoring, ID, or arbitrary long-format columns. The fitted schema is rebuilt from the rows used for the actual fit. Only the documented ranger whitelist can be forwarded through ...; sampling settings that leave no genuine OOB observations are rejected.

6. Model Summary

summary(fit)

The summary reports the exact fit controls and separate ranger OOB error and OOB concordance for every edge, together with verified OOB coverage and separate target-event, competing-exit, and external-censoring counts. Those edge metrics are not full-state validation.

7. Feature Importance

imp <- importance(fit)
print(imp)
plot(imp, type = "barplot")

Permutation importance is the transition-specific change in ranger OOB predictive loss after permuting a predictor. Negative values can arise from Monte Carlo noise, sparse events, correlated predictors, or irrelevant variables; they are not causal or protective effects. Event counts are stored beside the long-form importance values and should be considered when comparing edges.

plot(imp, type = "heatmap")

The heatmap contains the same edge-specific values as the preceding bar plot.

8. Predict for New Patients

newdata <- data.frame(
  age = c(50, 70),
  sex = c(0, 1),
  BMI = c(24, 32),
  treatment = c(1, 0)
)

prediction_horizon <- min(fit$max_duration_by_origin)
pred <- predict(fit, newdata = newdata,
                times = seq(0, prediction_horizon, length.out = 37))

# Plot for patient 1 (young, treated)
plot(pred, type = "state_occupation", subject = 1)

# Plot for patient 2 (older, untreated)
plot(pred, type = "state_occupation", subject = 2)

Both curves come from the same pred object and canonical fit. They are conditional on fresh entry into the initial state at elapsed duration zero; the public starting-state dimension contains only that requested state. They are not ongoing-sojourn dynamic predictions and have no confidence bands.

9. Diagnostics

diag <- diagnose(fit)
print(diag)

plot(diag, type = "concordance")

This figure visualizes genuine ranger OOB concordance separately for each binary edge endpoint.

Full-state Brier scores require patient-level cross-validation and refitting; they are never assembled from incompatible edge-level OOB predictions. Every fold rebuilds its predictor schema from training subjects only. A validation- only factor level stops the procedure rather than leaking full-data levels, and successful results retain exact subject assignments and refit seeds:

cv_diag <- diagnose(fit, method = "cv", folds = 5,
                    eval_times = seq(0, prediction_horizon * 0.8,
                                     length.out = 9))
plot(cv_diag, type = "brier")

10. Transition Diagram

plot_transition_diagram(ms, msdata)

The diagram uses the original display order and annotates each allowed edge with its observed event count.

Advanced probability assembly

compute_trans_prob() is the advanced public route for combining a complete, named set of clock-reset cumulative cause-specific hazard curves. The same validated solver is used by predict.rfmstate().

simple_ms <- define_multistate(c("A", "B"), "B", list(A = "B"))
elapsed_grid <- seq(0, 2, length.out = 2001)
simple_hazards <- list(
  "A->B" = data.frame(time = elapsed_grid,
                       hazard = 0.4 * elapsed_grid)
)
simple_prob <- compute_trans_prob(
  simple_hazards, simple_ms, times = c(0, 1, 2),
  target_grid_points = 512
)
simple_prob$state_occ

The output rows correspond to the requested elapsed durations and the columns to occupied states. The solver evaluates cumulative hazards as step functions, checks probability mass, and refines a regular grid without clipping or row normalization.

Input and prediction edge cases

• Left truncation/delayed entry and s != 0 are unsupported.
• Cycles, recurrent visits, tied/decreasing transition times, events after censoring or absorption, and nonpositive durations are errors.
• Fitting covariates must be complete. Prediction profiles must contain every stored predictor with compatible classes, finite values, and no unseen factor levels.
• Prediction times must lie within the conservative follow-up support. The explicit extrapolate = "flat" sensitivity option assumes zero additional hazard beyond support and is unsuitable for primary reported analyses.

Methodology

Separate probability constructions

RFmstate forests use duration since fresh entry into the current state. Their predicted cause-specific cumulative hazards are combined by semi-Markov entry-mass and sojourn convolution on a validated regular duration grid. The output is an entry-conditioned state-occupation array, not a general Markov $P(s,t)$ matrix.

The Aalen-Johansen baseline is separate: it uses calendar-time risk sets and a product integral from the recorded common study origin.

Aalen-Johansen Estimator (Nonparametric Baseline)

The Aalen-Johansen (AJ) estimator uses calendar-time risk sets and a product integral from the common baseline. RFmstate exposes point estimates as a descriptive population benchmark; it does not use AJ as the covariate-free form of the clock-reset forest solver. It estimates hazard increments via the Nelson--Aalen formula:

$$d\hat{A}{hj}(u) = \frac{dN{hj}(u)}{Y_h(u)}$$

where $dN_{hj}(u)$ counts the observed $h \to j$ transitions at time $u$ and $Y_h(u)$ is the number at risk in state $h$ just before time $u$. This provides population-level transition probabilities without covariate adjustment and serves as a covariate-free baseline in the package.

Random Forest Multistate Approach

For covariate-adjusted predictions, we decompose the multistate model into per-origin-state competing risks problems:

1. For each transient state $h$, identify all outgoing transitions
2. Fit a cause-specific RSF: For each destination state $j$, fit a random survival forest treating transition $h \to j$ as the target event. Other observed exits end the origin-state risk interval and receive a non-target indicator; they are competing events, not loss to follow-up
3. Extract cumulative hazards: Use ranger's predicted cumulative hazard for each declared edge on its observed duration support
4. Assemble by convolution: combine state-entry masses, origin-state sojourn survival, and cause-specific next-exit distributions in topological order

This approach leverages the flexibility of random forests to capture nonlinear covariate effects and interactions while maintaining the interpretability of the approved acyclic, non-recurrent multistate scope. Analytic, probability-invariant, and grid-refinement checks validate the numerical approximation without clipping or row normalization.

Diagnostics

• OOB Error: Out-of-bag prediction error from the random forest ensemble
• Edge OOB C-index: genuine ranger OOB concordance for each binary cause-specific endpoint; it is not a full-model validation score
• Full-state Brier/IBS: IPCW scores from patient-level held-out complete pipeline predictions

References

• Aalen, O.O. & Johansen, S. (1978). An empirical transition matrix for non-homogeneous Markov chains based on censored observations. Scandinavian Journal of Statistics, 5(3), 141-150.
• Ishwaran, H. et al. (2008). Random survival forests. Annals of Applied Statistics, 2(3), 841-860.
• Putter, H., Fiocco, M. & Geskus, R.B. (2007). Tutorial in biostatistics: Competing risks and multi-state models. Statistics in Medicine, 26, 2389-2430.

Reproducibility

This vignette renders from its source in a clean package checkout. It uses no external comparison_results.rds, private cache, or precomputed numerical result. Every displayed table and figure is generated by the code block that immediately precedes it using the canonical configuration declared in the hidden setup chunk. sessionInfo() records the rendering environment below.

sessionInfo()

RFmstate documentation built on Sept. 10, 2026, 1:09 a.m.