RFmstate Demo: Custom Clinical Data with Date Variables
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)
#> record_id gender trt weight trt_date cured_date heart_failure_date
#> 1 1 Male Drug A 74.9 2021-11-08 <NA> 2022-04-17
#> 2 2 Male Drug B 86.4 2021-04-26 2021-10-23 <NA>
#> 3 3 Male Drug B 75.6 2021-11-18 2022-03-06 <NA>
#> 4 4 Male Drug B 86.0 2021-10-31 2022-04-14 <NA>
#> 5 5 Female Drug A 72.8 2021-12-29 2022-10-06 2022-01-04
#> 6 6 Female Drug A 74.1 2020-04-03 2020-09-03 2020-04-06
#> 7 7 Female Drug B 82.2 2020-07-19 <NA> <NA>
#> 8 8 Female Drug A 89.9 2020-10-19 2021-02-06 <NA>
#> 9 9 Male Drug A 56.3 2021-12-21 <NA> 2022-06-04
#> 10 10 Female Drug B 74.5 2020-08-23 2020-12-01 <NA>
#> death_date last_followup_date
#> 1 <NA> 2022-12-21
#> 2 <NA> <NA>
#> 3 <NA> <NA>
#> 4 <NA> <NA>
#> 5 <NA> <NA>
#> 6 <NA> <NA>
#> 7 2020-08-24 <NA>
#> 8 <NA> <NA>
#> 9 <NA> 2023-05-24
#> 10 <NA> <NA>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)
#> record_id gender trt weight time_HeartFailure time_Cured time_Death
#> 1 1 1 1 74.9 160 NA NA
#> 2 2 1 0 86.4 NA 180 NA
#> 3 3 1 0 75.6 NA 108 NA
#> 4 4 1 0 86.0 NA 165 NA
#> 5 5 0 1 72.8 6 281 NA
#> 6 6 0 1 74.1 3 153 NA
#> 7 7 0 0 82.2 NA NA 36
#> 8 8 0 1 89.9 NA 110 NA
#> 9 9 1 1 56.3 165 NA NA
#> 10 10 0 0 74.5 NA 100 NA
#> time_censored
#> 1 408
#> 2 NA
#> 3 NA
#> 4 NA
#> 5 NA
#> 6 NA
#> 7 NA
#> 8 NA
#> 9 519
#> 10 NAQuick summary of event counts:
cat("Total patients:", nrow(dat), "\n")
#> Total patients: 300
cat("Heart failure observed:", sum(!is.na(dat$time_HeartFailure)), "\n")
#> Heart failure observed: 146
cat("Cured:", sum(!is.na(dat$time_Cured)), "\n")
#> Cured: 198
cat("Death:", sum(!is.na(dat$time_Death)), "\n")
#> Death: 95
cat("Censored (no absorbing state):", sum(!is.na(dat$time_censored)), "\n")
#> Censored (no absorbing state): 73. 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)
#> Multistate Structure
#> States: Treatment -> HeartFailure -> Cured -> Death
#> Absorbing: Cured, Death
#> Common initial state: Treatment
#> Computational order: Treatment -> HeartFailure -> Cured -> Death
#> Transitions: 5
#> 1: Treatment -> HeartFailure
#> 2: Treatment -> Cured
#> 3: Treatment -> Death
#> 4: HeartFailure -> Cured
#> 5: HeartFailure -> Death4. 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)
#> Multistate Data (msdata)
#> Patients: 300
#> Intervals: 446
#> Transitions observed: 439
#> Externally censored intervals: 7
#> Initial state: Treatment
#> States: Treatment, HeartFailure, Cured, Death
#> Approved baseline predictors: gender, trt, weight
#>
#> Transition counts:
#> to
#> from Treatment HeartFailure Cured Death
#> Treatment 0 146 103 49
#> HeartFailure 0 0 95 46
#> Cured 0 0 0 0
#> Death 0 0 0 0
#>
#> Per-edge outcomes (target / competing / external censoring):
#> from to n_events n_competing_exits n_external_censored
#> Treatment HeartFailure 146 152 2
#> Treatment Cured 103 195 2
#> Treatment Death 49 249 2
#> HeartFailure Cured 95 46 5
#> HeartFailure Death 46 95 55. Transition Diagram
plot_transition_diagram(ms, msdata)
State transition diagram with event counts.
6. Aalen-Johansen Nonparametric Estimates
aj <- aalen_johansen(msdata)
print(aj)
#> Aalen-Johansen Estimate
#> Time range: [1, 632]
#> Event times: 251
#> States: Treatment, HeartFailure, Cured, Death
#> Common initial state: Treatment
#> Uncertainty: point estimates only
#>
#> Event counts per transition:
#> from to n_events
#> Treatment HeartFailure 146
#> Treatment Cured 103
#> Treatment Death 49
#> HeartFailure Cured 95
#> HeartFailure Death 46
#>
#> Final state occupation probabilities:
#> Treatment: 0
#> HeartFailure: 0
#> Cured: 0.6771
#> Death: 0.3229plot(aj, type = "state_occupation")
State occupation probabilities (Aalen-Johansen).
plot(aj, type = "cumulative_hazard")
Nelson-Aalen cumulative hazards by transition.
plot(aj, type = "stacked_transition_prob")
Transition probabilities from Treatment state (AJ).
plot(aj, type = "hazard_increment")
Nelson-Aalen hazard increments over time (AJ).
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)
#> Clock-Reset Semi-Markov Random-Forest Model
#> Common initial state: Treatment
#> Time scale: duration since fresh state entry
#> Covariates: gender, trt, weight
#> Trees per edge: 200
#> min_events safeguard: 3
#>
#> Fitted edge models:
#> Treatment->HeartFailure: 300 sojourns, 146 target, 152 competing, 2 externally censored; OOB C = 0.6117, OOB coverage = 1.000
#> Treatment->Cured: 300 sojourns, 103 target, 195 competing, 2 externally censored; OOB C = 0.4815, OOB coverage = 1.000
#> Treatment->Death: 300 sojourns, 49 target, 249 competing, 2 externally censored; OOB C = 0.3867, OOB coverage = 1.000
#> HeartFailure->Cured: 146 sojourns, 95 target, 46 competing, 5 externally censored; OOB C = 0.4608, OOB coverage = 1.000
#> HeartFailure->Death: 146 sojourns, 46 target, 95 competing, 5 externally censored; OOB C = 0.4809, OOB coverage = 1.0008. Model Summary
s <- summary(fit)9. Feature Importance
imp <- importance(fit)
print(imp)
#> Feature Importance per Transition
#> ============================================================
#>
#> Treatment->HeartFailure Treatment->Cured Treatment->Death
#> gender -0.0015 -0.0085 -0.0185
#> trt 0.0263 0.0340 -0.0180
#> weight 0.0455 -0.0025 -0.0147
#> HeartFailure->Cured HeartFailure->Death
#> gender -0.0001 -0.0079
#> trt 0.0057 -0.0120
#> weight -0.0088 0.0175
#>
#> Top variables per transition:
#> Treatment->HeartFailure: weight (0.0455)
#> Treatment->Cured: trt (0.034)
#> Treatment->Death: weight (-0.0147)
#> HeartFailure->Cured: trt (0.0057)
#> HeartFailure->Death: weight (0.0175)
plot(imp, type = "barplot")
Feature importance per transition.
plot(imp, type = "heatmap")
Feature importance 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)
#> gender trt weight
#> Light male, Drug A 1 1 65
#> Heavy female, Drug B 0 0 90
#> Average male, Drug A 1 1 75
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)
#> Entry-Conditioned RFmstate Predictions
#> Profiles: 3
#> Starting state: Treatment
#> Conditioning: fresh entry at duration zero
#> Time scale: clock-reset
#> Prediction type: new profiles
#> Elapsed-time range: [0, 541]
#> Extrapolation: noneplot(pred, type = "state_occupation", subject = 1)
Predicted state occupation: light male on Drug A.
plot(pred, type = "state_occupation", subject = 2)
Predicted state occupation: heavy female on Drug B.
plot(pred, type = "state_occupation", subject = 3)
Predicted state occupation: average male on Drug A.
plot(pred, type = "transition_prob", subject = 1)
Predicted transition probabilities: light male on Drug A.
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)
#> RFmstate Diagnostics
#> Validation label: edge-level ranger OOB only
#>
#> Genuine ranger edge OOB concordance:
#> transition n_target_events n_competing_exits n_external_censored
#> Treatment->HeartFailure 146 152 2
#> Treatment->Cured 103 195 2
#> Treatment->Death 49 249 2
#> HeartFailure->Cured 95 46 5
#> HeartFailure->Death 46 95 5
#> prediction_error oob_concordance oob_fraction replace sample_fraction
#> 0.3882907 0.6117093 1 TRUE 1
#> 0.5185383 0.4814617 1 TRUE 1
#> 0.6133056 0.3866944 1 TRUE 1
#> 0.5391995 0.4608005 1 TRUE 1
#> 0.5191125 0.4808875 1 TRUE 1plot(diag, type = "concordance")
Concordance index by transition.
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)
}