Introduction to RFmstate

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:

Quick Start

1. Define the Multistate Structure

library(RFmstate)

# Use the built-in clinical trial structure
ms <- clinical_states()
print(ms)
#> Multistate Structure
#>   States: Baseline -> Responded -> Unresponded -> Stabilized -> Progressed -> Death 
#>   Absorbing: Death 
#>   Common initial state: Baseline 
#>   Computational order: Baseline -> Responded -> Unresponded -> Stabilized -> Progressed -> Death 
#>   Transitions: 12 
#>     1: Baseline -> Responded
#>     2: Baseline -> Unresponded
#>     3: Baseline -> Death
#>     4: Responded -> Stabilized
#>     5: Responded -> Progressed
#>     6: Responded -> Death
#>     7: Unresponded -> Stabilized
#>     8: Unresponded -> Progressed
#>     9: Unresponded -> Death
#>     10: Stabilized -> Progressed
#>     11: Stabilized -> Death
#>     12: Progressed -> Death

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)
#>   ID  age sex  BMI treatment time_Responded time_Unresponded time_Stabilized
#> 1  1 76.5   0 19.5         0      23.785816               NA        32.56660
#> 2  2 53.2   0 24.7         1       1.955828               NA        63.72278
#> 3  3 64.4   1 26.9         0             NA               NA              NA
#> 4  4 67.6   0 24.0         1      26.344666               NA              NA
#> 5  5 64.9   1 26.5         0             NA         19.76615              NA
#> 6  6 58.7   0 24.4         1      23.633832               NA              NA
#>   time_Progressed time_Death time_censored
#> 1              NA         NA      123.2743
#> 2       194.16850  251.42542            NA
#> 3              NA   53.25271            NA
#> 4        66.17295  130.56225            NA
#> 5        59.07904         NA      203.1828
#> 6              NA   42.33299            NA

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)
#> Multistate Data (msdata)
#>   Patients: 300 
#>   Intervals: 865 
#>   Transitions observed: 748 
#>   Externally censored intervals: 117 
#>   Initial state: Baseline 
#>   States: Baseline, Responded, Unresponded, Stabilized, Progressed, Death 
#>   Approved baseline predictors: age, sex, BMI, treatment 
#> 
#> Transition counts:
#>              to
#> from          Baseline Responded Unresponded Stabilized Progressed Death
#>   Baseline           0       190          85          0          0    25
#>   Responded          0         0           0        118         46    18
#>   Unresponded        0         0           0         37         30    17
#>   Stabilized         0         0           0          0         59    33
#>   Progressed         0         0           0          0          0    90
#>   Death              0         0           0          0          0     0
#> 
#> Per-edge outcomes (target / competing / external censoring):
#>         from          to n_events n_competing_exits n_external_censored
#>     Baseline   Responded      190               110                   0
#>     Baseline Unresponded       85               215                   0
#>     Baseline       Death       25               275                   0
#>    Responded  Stabilized      118                64                   8
#>    Responded  Progressed       46               136                   8
#>    Responded       Death       18               164                   8
#>  Unresponded  Stabilized       37                47                   1
#>  Unresponded  Progressed       30                54                   1
#>  Unresponded       Death       17                67                   1
#>   Stabilized  Progressed       59                33                  63
#>   Stabilized       Death       33                59                  63
#>   Progressed       Death       90                 0                  45
head(msdata)
#>   id       from         to    Tstart      Tstop status trans_id  age sex  BMI
#> 1  1   Baseline  Responded  0.000000  23.785816      1        1 76.5   0 19.5
#> 2  1  Responded Stabilized 23.785816  32.566596      1        4 76.5   0 19.5
#> 3  1 Stabilized       <NA> 32.566596 123.274260      0       NA 76.5   0 19.5
#> 4  2   Baseline  Responded  0.000000   1.955828      1        1 53.2   0 24.7
#> 5  2  Responded Stabilized  1.955828  63.722784      1        4 53.2   0 24.7
#> 6  2 Stabilized Progressed 63.722784 194.168502      1       10 53.2   0 24.7
#>   treatment   duration
#> 1         0  23.785816
#> 2         0   8.780781
#> 3         0  90.707664
#> 4         1   1.955828
#> 5         1  61.766956
#> 6         1 130.445717

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)
#> Aalen-Johansen Estimate
#>   Time range: [0.04004214, 331.811]
#>   Event times: 748 
#>   States: Baseline, Responded, Unresponded, Stabilized, Progressed, Death 
#>   Common initial state: Baseline 
#>   Uncertainty: point estimates only
#> 
#> Event counts per transition:
#>         from          to n_events
#>     Baseline   Responded      190
#>     Baseline Unresponded       85
#>     Baseline       Death       25
#>    Responded  Stabilized      118
#>    Responded  Progressed       46
#>    Responded       Death       18
#>  Unresponded  Stabilized       37
#>  Unresponded  Progressed       30
#>  Unresponded       Death       17
#>   Stabilized  Progressed       59
#>   Stabilized       Death       33
#>   Progressed       Death       90
#> 
#> Final state occupation probabilities:
#>   Baseline: 0
#>   Responded: 0
#>   Unresponded: 0.0033
#>   Stabilized: 0.1234
#>   Progressed: 0.0297
#>   Death: 0.8435
plot(aj, type = "state_occupation")
State occupation probabilities from Aalen-Johansen estimator
State occupation probabilities from Aalen-Johansen estimator

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")
Nelson-Aalen cumulative hazards
Nelson-Aalen cumulative hazards

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
)
#> Warning: Sparse transition(s); interpret edge diagnostics cautiously:
#> Responded->Death (18 events), Unresponded->Death (17 events). The threshold is
#> descriptive, not a universal adequacy rule.
print(fit)
#> Clock-Reset Semi-Markov Random-Forest Model
#>   Common initial state: Baseline 
#>   Time scale: duration since fresh state entry
#>   Covariates: age, sex, BMI, treatment 
#>   Trees per edge: 200 
#>   min_events safeguard: 3 
#> 
#> Fitted edge models:
#>   Baseline->Responded: 300 sojourns, 190 target, 110 competing, 0 externally censored; OOB C = 0.6050, OOB coverage = 1.000
#>   Baseline->Unresponded: 300 sojourns, 85 target, 215 competing, 0 externally censored; OOB C = 0.5058, OOB coverage = 1.000
#>   Baseline->Death: 300 sojourns, 25 target, 275 competing, 0 externally censored; OOB C = 0.6004, OOB coverage = 1.000
#>   Responded->Stabilized: 190 sojourns, 118 target, 64 competing, 8 externally censored; OOB C = 0.5267, OOB coverage = 1.000
#>   Responded->Progressed: 190 sojourns, 46 target, 136 competing, 8 externally censored; OOB C = 0.6175, OOB coverage = 1.000
#>   Responded->Death: 190 sojourns, 18 target, 164 competing, 8 externally censored; OOB C = 0.6228, OOB coverage = 1.000
#>   Unresponded->Stabilized: 85 sojourns, 37 target, 47 competing, 1 externally censored; OOB C = 0.5228, OOB coverage = 1.000
#>   Unresponded->Progressed: 85 sojourns, 30 target, 54 competing, 1 externally censored; OOB C = 0.5910, OOB coverage = 1.000
#>   Unresponded->Death: 85 sojourns, 17 target, 67 competing, 1 externally censored; OOB C = 0.6903, OOB coverage = 1.000
#>   Stabilized->Progressed: 155 sojourns, 59 target, 33 competing, 63 externally censored; OOB C = 0.6060, OOB coverage = 1.000
#>   Stabilized->Death: 155 sojourns, 33 target, 59 competing, 63 externally censored; OOB C = 0.4534, OOB coverage = 1.000
#>   Progressed->Death: 135 sojourns, 90 target, 0 competing, 45 externally censored; OOB C = 0.5948, OOB coverage = 1.000

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)
#> Random Forest Multistate Model Summary
#> ================================================== 
#> 
#> Call: rfmstate(msdata = msdata, num.trees = canonical_trees, min_events = canonical_min_events, 
#>     sparse_warning = canonical_sparse_warning, seed = canonical_seed)
#> 
#> Data:
#>   Patients: 300 
#>   Total transitions: 748 
#>   Total intervals: 865 
#> 
#> Covariates: age, sex, BMI, treatment 
#>   Time scale: clock-reset duration; semi-Markov assembly
#> 
#> Forest parameters:
#>   Trees: 200 
#>   mtry: 2 
#>   Min node size: 15 
#>   min_events safeguard: 3 
#>   Forwarded ranger arguments: none 
#>   Effective OOB sampling: replace = TRUE ; sample.fraction = 1 ; oob.error = TRUE; keep.inbag = TRUE
#> 
#> Transition-specific models:
#> ---------------------------------------------------------------------------------------------------- 
#> Transition            Total  Target Compete  Extern  OOB Error      OOB C OOB Frac
#> ---------------------------------------------------------------------------------------------------- 
#> Baseline -> Responded    300     190     110       0     0.3950     0.6050    1.000
#> Baseline -> Unresponded    300      85     215       0     0.4942     0.5058    1.000
#> Baseline -> Death       300      25     275       0     0.3996     0.6004    1.000
#> Responded -> Stabilized    190     118      64       8     0.4733     0.5267    1.000
#> Responded -> Progressed    190      46     136       8     0.3825     0.6175    1.000
#> Responded -> Death      190      18     164       8     0.3772     0.6228    1.000
#> Unresponded -> Stabilized     85      37      47       1     0.4772     0.5228    1.000
#> Unresponded -> Progressed     85      30      54       1     0.4090     0.5910    1.000
#> Unresponded -> Death     85      17      67       1     0.3097     0.6903    1.000
#> Stabilized -> Progressed    155      59      33      63     0.3940     0.6060    1.000
#> Stabilized -> Death     155      33      59      63     0.5466     0.4534    1.000
#> Progressed -> Death     135      90       0      45     0.4052     0.5948    1.000
#> ----------------------------------------------------------------------------------------------------

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)
#> Feature Importance per Transition
#> ============================================================ 
#> 
#>           Baseline->Responded Baseline->Unresponded Baseline->Death
#> age                    0.0014                0.0045         -0.0086
#> sex                   -0.0004               -0.0004         -0.0048
#> BMI                    0.0171               -0.0173          0.0284
#> treatment              0.0466                0.0087         -0.0127
#>           Responded->Stabilized Responded->Progressed Responded->Death
#> age                     -0.0089                0.0329           0.0790
#> sex                     -0.0034                0.0095          -0.0122
#> BMI                      0.0118                0.0463          -0.0052
#> treatment                0.0061                0.0051           0.0038
#>           Unresponded->Stabilized Unresponded->Progressed Unresponded->Death
#> age                        0.0158                  0.0206             0.1281
#> sex                        0.0011                 -0.0049            -0.0118
#> BMI                       -0.0125                  0.0020             0.0052
#> treatment                  0.0036                  0.0326            -0.0100
#>           Stabilized->Progressed Stabilized->Death Progressed->Death
#> age                       0.0435            0.0028            0.0341
#> sex                      -0.0029           -0.0091           -0.0033
#> BMI                       0.0173            0.0331           -0.0024
#> treatment                 0.0070           -0.0111            0.0137
#> 
#> Top variables per transition:
#>   Baseline->Responded: treatment (0.0466)
#>   Baseline->Unresponded: treatment (0.0087)
#>   Baseline->Death: BMI (0.0284)
#>   Responded->Stabilized: BMI (0.0118)
#>   Responded->Progressed: BMI (0.0463)
#>   Responded->Death: age (0.079)
#>   Unresponded->Stabilized: age (0.0158)
#>   Unresponded->Progressed: treatment (0.0326)
#>   Unresponded->Death: age (0.1281)
#>   Stabilized->Progressed: age (0.0435)
#>   Stabilized->Death: BMI (0.0331)
#>   Progressed->Death: age (0.0341)
plot(imp, type = "barplot")
Feature importance per transition
Feature importance per transition

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")
Feature importance heatmap
Feature importance 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)
Predicted state occupation for two patient profiles
Predicted state occupation for two patient profiles

# Plot for patient 2 (older, untreated)
plot(pred, type = "state_occupation", subject = 2)
Predicted state occupation for two patient profiles
Predicted state occupation for two patient profiles

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)
#> RFmstate Diagnostics
#>   Validation label: edge-level ranger OOB only 
#> 
#> Genuine ranger edge OOB concordance:
#>               transition n_target_events n_competing_exits n_external_censored
#>      Baseline->Responded             190               110                   0
#>    Baseline->Unresponded              85               215                   0
#>          Baseline->Death              25               275                   0
#>    Responded->Stabilized             118                64                   8
#>    Responded->Progressed              46               136                   8
#>         Responded->Death              18               164                   8
#>  Unresponded->Stabilized              37                47                   1
#>  Unresponded->Progressed              30                54                   1
#>       Unresponded->Death              17                67                   1
#>   Stabilized->Progressed              59                33                  63
#>        Stabilized->Death              33                59                  63
#>        Progressed->Death              90                 0                  45
#>  prediction_error oob_concordance oob_fraction replace sample_fraction
#>         0.3949897       0.6050103            1    TRUE               1
#>         0.4942148       0.5057852            1    TRUE               1
#>         0.3995732       0.6004268            1    TRUE               1
#>         0.4732574       0.5267426            1    TRUE               1
#>         0.3824935       0.6175065            1    TRUE               1
#>         0.3772016       0.6227984            1    TRUE               1
#>         0.4772036       0.5227964            1    TRUE               1
#>         0.4090136       0.5909864            1    TRUE               1
#>         0.3097463       0.6902537            1    TRUE               1
#>         0.3939962       0.6060038            1    TRUE               1
#>         0.5465740       0.4534260            1    TRUE               1
#>         0.4051637       0.5948363            1    TRUE               1
plot(diag, type = "concordance")
Concordance index per transition
Concordance index per transition

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)
Transition diagram with event counts
Transition diagram with event counts

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
#>             occupied_state
#> elapsed_time        A             B
#>            0 1.000000 -1.856154e-16
#>            1 0.670320  3.296800e-01
#>            2 0.449329  5.506710e-01

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

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

References

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()
#> R version 4.6.0 (2026-04-24)
#> Platform: aarch64-apple-darwin23
#> Running under: macOS Tahoe 26.6.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
#> 
#> time zone: America/Chicago
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] RFmstate_0.1.9
#> 
#> loaded via a namespace (and not attached):
#>  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   Matrix_1.7-5   
#>  [5] ranger_0.18.0   xfun_0.60       lattice_0.22-9  splines_4.6.0  
#>  [9] cachem_1.1.0    knitr_1.51      htmltools_0.5.9 rmarkdown_2.31 
#> [13] lifecycle_1.0.5 cli_3.6.6       grid_4.6.0      sass_0.4.10    
#> [17] jquerylib_0.1.4 compiler_4.6.0  tools_4.6.0     evaluate_1.0.5 
#> [21] bslib_0.12.0    survival_3.8-6  Rcpp_1.1.2      yaml_2.3.12    
#> [25] otel_0.2.0      rlang_1.3.0     jsonlite_2.0.0