RFmstate Demo: Custom Clinical Data with Date Variables

In RFmstate: Random Forest-Based Multistate Survival Analysis

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

canonical_seed <- 42L
canonical_n <- 300L
canonical_trees <- 200L
canonical_min_events <- 3L
canonical_sparse_warning <- 20L
canonical_covariates <- c("gender", "trt", "weight")

All sections use one configuration: seed 42, 300 subjects, 200 trees per edge, and the predictor set declared during data preparation. These values keep the source render practical and are not event-count or tuning recommendations.

1. Simulate Raw Clinical Data

We simulate a dataset with date-based event times, as one might receive from a clinical database. Patients start treatment and may experience heart failure, be cured, or die. Cured and Death are absorbing states; Heart Failure is a transient intermediate state.

set.seed(canonical_seed)
n <- canonical_n

record_id <- seq_len(n)
gender <- sample(c("Male", "Female"), n, replace = TRUE)
trt <- sample(c("Drug A", "Drug B"), n, replace = TRUE)
weight <- round(rnorm(n, mean = 75, sd = 15), 1)

# Treatment start dates spread over 2 years
trt_date <- as.Date("2020-01-01") + sample(0:730, n, replace = TRUE)

# Simulate trajectories manually
heart_failure_date <- rep(as.Date(NA), n)
cured_date <- rep(as.Date(NA), n)
death_date <- rep(as.Date(NA), n)
last_followup_date <- rep(as.Date(NA), n)

for (i in seq_len(n)) {
  # From Treatment state, competing events:
  #   -> Heart Failure (rate depends on weight, treatment)
  #   -> Cured (rate depends on treatment)
  #   -> Death (rate depends on weight)
  trt_effect <- ifelse(trt[i] == "Drug A", 0.7, 1.0)
  wt_effect <- exp((weight[i] - 75) / 50)

  t_hf <- rweibull(1, shape = 1.3, scale = 300 * trt_effect * wt_effect)
  t_cured <- rweibull(1, shape = 1.5, scale = 250 / trt_effect)
  t_death <- rweibull(1, shape = 1.0, scale = 800 * (1 / wt_effect))

  first_wait_days <- pmax(1L, ceiling(c(t_hf, t_cured, t_death)))
  first_event <- which.min(first_wait_days)
  first_day <- first_wait_days[first_event]

  # Censoring at 3 years. Integer-day waits are rounded upward and constrained
  # to at least one day so the date representation cannot create zero-duration
  # sojourns.
  cens_day <- max(1L, ceiling(runif(1, 400, 1095)))

  if (first_day >= cens_day) {
    # Censored from Treatment state
    last_followup_date[i] <- trt_date[i] + cens_day
    next
  }

  if (first_event == 1) {
    # Heart Failure reached
    heart_failure_date[i] <- trt_date[i] + first_day

    # From Heart Failure: -> Cured or -> Death
    t_cured2 <- rweibull(1, shape = 1.4, scale = 200 / trt_effect)
    t_death2 <- rweibull(1, shape = 1.2, scale = 400 * (1 / wt_effect))

    second_wait_days <- pmax(1L, ceiling(c(t_cured2, t_death2)))
    second_event <- which.min(second_wait_days)
    second_day <- first_day + second_wait_days[second_event]

    if (second_day >= cens_day) {
      last_followup_date[i] <- trt_date[i] + cens_day
    } else if (second_event == 1) {
      cured_date[i] <- trt_date[i] + second_day
    } else {
      death_date[i] <- trt_date[i] + second_day
    }
  } else if (first_event == 2) {
    cured_date[i] <- trt_date[i] + first_day
  } else {
    death_date[i] <- trt_date[i] + first_day
  }
}

raw_data <- data.frame(
  record_id = record_id,
  gender = gender,
  trt = trt,
  weight = weight,
  trt_date = trt_date,
  cured_date = cured_date,
  heart_failure_date = heart_failure_date,
  death_date = death_date,
  last_followup_date = last_followup_date,
  stringsAsFactors = FALSE
)

head(raw_data, 10)

2. Compute Time-to-Event from Treatment Date

Convert date columns to days since treatment start.

dat <- data.frame(
  record_id = raw_data$record_id,
  gender = as.integer(raw_data$gender == "Male"),
  trt = as.integer(raw_data$trt == "Drug A"),
  weight = raw_data$weight,
  time_HeartFailure = as.numeric(
    difftime(raw_data$heart_failure_date, raw_data$trt_date, units = "days")
  ),
  time_Cured = as.numeric(
    difftime(raw_data$cured_date, raw_data$trt_date, units = "days")
  ),
  time_Death = as.numeric(
    difftime(raw_data$death_date, raw_data$trt_date, units = "days")
  ),
  time_censored = as.numeric(
    difftime(raw_data$last_followup_date, raw_data$trt_date, units = "days")
  ),
  stringsAsFactors = FALSE
)

head(dat, 10)

Quick summary of event counts:

cat("Total patients:", nrow(dat), "\n")
cat("Heart failure observed:", sum(!is.na(dat$time_HeartFailure)), "\n")
cat("Cured:", sum(!is.na(dat$time_Cured)), "\n")
cat("Death:", sum(!is.na(dat$time_Death)), "\n")
cat("Censored (no absorbing state):", sum(!is.na(dat$time_censored)), "\n")

3. Define Multistate Structure

library(RFmstate)

ms <- define_multistate(
  state_names = c("Treatment", "HeartFailure", "Cured", "Death"),
  absorbing = c("Cured", "Death"),
  transitions = list(
    Treatment = c("HeartFailure", "Cured", "Death"),
    HeartFailure = c("Cured", "Death")
  )
)
print(ms)

4. Prepare Multistate Data

msdata <- prepare_data(
  data = dat,
  id = "record_id",
  structure = ms,
  time_map = list(
    HeartFailure = "time_HeartFailure",
    Cured = "time_Cured",
    Death = "time_Death"
  ),
  censor_col = "time_censored",
  covariates = canonical_covariates
)
print(msdata)

5. Transition Diagram

plot_transition_diagram(ms, msdata)

6. Aalen-Johansen Nonparametric Estimates

aj <- aalen_johansen(msdata)
print(aj)

plot(aj, type = "state_occupation")

plot(aj, type = "cumulative_hazard")

plot(aj, type = "stacked_transition_prob")

plot(aj, type = "hazard_increment")

7. 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)

8. Model Summary

s <- summary(fit)

9. Feature Importance

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

plot(imp, type = "heatmap")

10. Predict for New Patients

new_patients <- data.frame(
  gender = c(1, 0, 1),
  trt = c(1, 0, 1),
  weight = c(65, 90, 75)
)
rownames(new_patients) <- c("Light male, Drug A",
                             "Heavy female, Drug B",
                             "Average male, Drug A")
print(new_patients)

prediction_horizon <- floor(min(fit$max_duration_by_origin))
prediction_times <- sort(unique(c(
  0, seq(30, prediction_horizon, by = 30), prediction_horizon
)))
# This date-based example has integer-day event times, so an integer-day
# initial grid aligns with the observed hazard jumps before refinement.
pred <- predict(fit, newdata = new_patients, times = prediction_times,
                grid_step = 1)
print(pred)

plot(pred, type = "state_occupation", subject = 1)

plot(pred, type = "state_occupation", subject = 2)

plot(pred, type = "state_occupation", subject = 3)

plot(pred, type = "transition_prob", subject = 1)

All four prediction figures come from the same fitted object and prediction grid. They are conditional on fresh entry into Treatment at elapsed day zero; the transition-probability view exposes only that selected starting-state row, not a general all-start-state Markov matrix.

11. Diagnostics

diag <- diagnose(fit)
print(diag)

plot(diag, type = "concordance")

Patient-level cross-validation is required for full-state IPCW Brier scores:

Each refit learns factor levels and numeric ranges from its training subjects only. A held-out-only factor level fails explicitly, while successful results retain exact subject assignments, assignment/refit seeds, per-edge event counts, support, and censoring-stability metadata. Confirmatory analyses should supply eval_times explicitly as below.

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

12. Comparing Treatment Arms

We can compare predicted outcomes between Drug A and Drug B for an average patient.

drug_a <- data.frame(gender = 1, trt = 1, weight = 75)
drug_b <- data.frame(gender = 1, trt = 0, weight = 75)

pred_a <- predict(fit, newdata = drug_a, times = prediction_times,
                  grid_step = 1)
pred_b <- predict(fit, newdata = drug_b, times = prediction_times,
                  grid_step = 1)

times <- pred_a$time
states <- ms$state_names

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
cols <- c("#1b9e77", "#d95f02")
for (j in seq_along(states)) {
  occ_a <- pred_a$state_occ[1, j, ]
  occ_b <- pred_b$state_occ[1, j, ]
  plot(times, occ_a, type = "l", col = cols[1], lwd = 2,
       ylim = c(0, max(c(occ_a, occ_b)) * 1.1),
       xlab = "Days", ylab = "Probability",
       main = states[j])
  lines(times, occ_b, col = cols[2], lwd = 2)
  legend("topright", legend = c("Drug A", "Drug B"),
         col = cols, lwd = 2, bty = "n", cex = 0.8)
}

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