--- title: "Random Hazard Forests with ggRandomForests" author: "John Ehrlinger" date: today format: html: fig-format: png fig-dpi: 96 toc: true toc-depth: 3 html-math-method: mathjax code-fold: false bibliography: ggRandomForests.bib vignette: > %\VignetteIndexEntry{Random Hazard Forests with ggRandomForests} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} # vignette dir under R CMD check, package root under an interactive knit .fo <- c("_fig_optim.R", file.path("vignettes", "_fig_optim.R")) if (any(file.exists(.fo))) source(.fo[file.exists(.fo)][1]) knitr::opts_chunk$set( fig.align = "center", fig.width = 7, fig.height = 4.5, message = FALSE, warning = FALSE, comment = "#>" ) options(mc.cores = 1, rf.cores = 1) ``` ```{r packages, include=FALSE} .has_rhf_extractor <- requireNamespace( "ggRandomForests", quietly = TRUE ) && "gg_rhf" %in% getNamespaceExports("ggRandomForests") if (.has_rhf_extractor) { library(ggRandomForests) } else if (requireNamespace("pkgload", quietly = TRUE)) { pkgload::load_all( export_all = FALSE, helpers = FALSE, attach_testthat = FALSE ) } else { stop( "Install ggRandomForests (or pkgload for dev builds) to render this vignette." ) } ``` ```{r load-bundle, include=FALSE, cache=FALSE} .required_bundle_names <- c( "data", "fit", "auct_cumulative", "auct_incident", "importance", "tune_risk", "tune_iauc", "seed", "settings", "versions" ) .required_setting_names <- c( "formula", "fit", "auct_cumulative", "auct_incident", "importance_cache", "importance_time_index", "tune_risk", "tune_iauc" ) .valid_rhf_bundle <- function(x) { if (!is.list(x)) return(FALSE) if (!identical(names(x), .required_bundle_names)) return(FALSE) checks <- c( is.data.frame(x$data), inherits(x$fit, "rhf"), inherits(x$auct_cumulative, "auct.rhf"), inherits(x$auct_incident, "auct.rhf"), inherits(x$importance, "importance.rhf"), inherits(x$tune_risk, "tune.treesize.rhf"), inherits(x$tune_iauc, "tune.treesize.rhf"), is.integer(x$seed), length(x$seed) == 1L, is.list(x$settings), identical(names(x$settings), .required_setting_names), is.character(x$versions), length(x$versions) > 0L ) all(checks) } .bundle_paths <- c( "rhf_precomputed.rds", file.path("vignettes", "rhf_precomputed.rds") ) bundle <- NULL for (.bundle_path in .bundle_paths[file.exists(.bundle_paths)]) { .candidate_bundle <- tryCatch( readRDS(.bundle_path), error = function(e) NULL ) if (.valid_rhf_bundle(.candidate_bundle)) { bundle <- .candidate_bundle break } } if (is.null(bundle)) { stop( paste( "The RHF vignette needs a valid precomputed bundle.", "Run vignettes/precompute_rhf.R from the package root." ) ) } ``` Most survival data begin with one row per subject: a follow-up time, an event indicator, and covariates measured at baseline. That row cannot show a covariate changing during follow-up. A Random Hazard Forest (RHF) instead accepts repeated time intervals for each subject, so a tree can use the covariate values that belong to each part of the subject's path [@Ishwaran:RHF:2026; @Ishwaran:RHF:software:2026]. This vignette uses simulated counting-process data to show how **ggRandomForests** reads an RHF fit. We first make the time-dependent data structure explicit, then extract case-specific hazard and cumulative hazard curves. The fitting, performance, importance, and tuning work has already been run by `vignettes/precompute_rhf.R`. The article evaluates only the inexpensive **ggRandomForests** extraction and plotting steps. ## What changes with counting-process data? One subject can contribute several rows. Each row covers an interval $(\mathit{start}, \mathit{stop}]$, and `event` records whether the event occurs at that interval's right endpoint. The first simulated subject makes the layout concrete: ```{r subject-intervals} rhf_data <- bundle$data one_subject <- rhf_data[rhf_data$id == rhf_data$id[1L], c("id", "start", "stop", "event", "xtd")] knitr::kable(one_subject, digits = 3) ``` The intervals do more than divide follow-up into pieces. They define a covariate path. A time-dependent covariate must be *predictable*: the value used at time $t$ is known just before the event decision at $t$. In the first RHF simulation, the actual identity is `xtd = (x.4 + x.5) * stop`. Because this is a continuous function of time and the subject's fixed `x.4` and `x.5`, its value at `stop` is also its left-hand limit there. It is available for the interval ending at that time. In data collected in practice, a lab value measured after an interval ends cannot be copied backward into that interval. At any candidate split time, the subject must be routed using the active record, not a later measurement. This is the no-lookahead rule: future covariate values never select an earlier branch. The active-record stitching and tree routing are behavior of `randomForestRHF::rhf()`, not **ggRandomForests**. **ggRandomForests** reads the fitted estimates after the upstream model has applied those rules [@Ishwaran:RHF:2026; @Ishwaran:RHF:software:2026]. Here is a compact view of the full simulated data set: ```{r data-summary} records_per_subject <- as.integer(table(rhf_data$id)) events_by_subject <- tapply(rhf_data$event, rhf_data$id, max) data_summary <- data.frame( measure = c( "Subjects", "Counting-process records", "Records per subject", "Subjects with an event", "xtd range" ), value = c( length(unique(rhf_data$id)), nrow(rhf_data), sprintf( "%d to %d (median %g)", min(records_per_subject), max(records_per_subject), stats::median(records_per_subject) ), sum(events_by_subject), sprintf("%.3f to %.3f", min(rhf_data$xtd), max(rhf_data$xtd)) ) ) knitr::kable(data_summary, col.names = c("Measure", "Value")) ``` ## Where did the fitted objects come from? The preparation below is shown but not evaluated while the vignette renders. These are the fixed calls used to create the saved bundle. Keeping the upstream objects in `rhf_precomputed.rds` makes CRAN builds predictable, while `vignettes/precompute_rhf.R` lets you reproduce the bundle from the same public simulation whenever you want to rerun the longer calculations. ```{r preparation, eval=FALSE} set.seed(20260825L) sim <- randomForestRHF::hazard.simulation(1) data <- sim$dta formula <- stats::as.formula("Surv(id, start, stop, event) ~ .") with(data, stopifnot(identical(xtd, (x.4 + x.5) * stop))) fit <- randomForestRHF::rhf( formula, data, ntree = 50L, seed = -1L ) auct_cumulative <- randomForestRHF::auct.rhf( fit, marker = "cumhaz", method = "cumulative", verbose = FALSE ) auct_incident <- randomForestRHF::auct.rhf( fit, marker = "hazard", method = "incident", riskset = "subject", verbose = FALSE ) cache <- randomForestRHF::varpro.cache.rhf( fit, max.rules.tree = 30L, max.tree = 20L, verbose = FALSE ) time_index <- unique(as.integer(round(seq.int( 1L, cache$K, length.out = 5L )))) importance <- randomForestRHF::importance.rhf( fit, cache = cache, time.index = time_index, verbose = FALSE ) tune_risk <- randomForestRHF::tune.treesize.rhf( formula, data, ntree = 20L, perf = "risk", lower = 2L, upper = 6L, max.evals = 5L, seed = 20260825L, verbose = FALSE, forest = FALSE ) tune_iauc <- randomForestRHF::tune.iAUC.rhf( formula, data, ntree = 20L, lower = 2L, upper = 6L, max.evals = 5L, seed = 20260825L, verbose = FALSE, forest = FALSE ) ``` The loaded fit can be checked without growing another forest: ```{r fit-summary} fit_summary <- data.frame( measure = c( "Object class", "Trees", "Subjects", "Working time points", "Saved ensemble source" ), value = c( paste(class(bundle$fit), collapse = ", "), bundle$fit$ntree, length(bundle$fit$ensemble.id), length(bundle$fit$time.interest), "out-of-bag hazard and cumulative hazard" ) ) knitr::kable(fit_summary, col.names = c("Fit component", "Value")) ``` ## How do hazard and cumulative hazard differ? `gg_rhf()` extracts both estimates into one tidy object, with one row for each subject and point on the RHF working time grid. This extraction runs live when the article renders: ```{r extract-rhf} rhf_curves <- ggRandomForests::gg_rhf(bundle$fit) head(as.data.frame(rhf_curves)) ``` The `source` column is `"oob"` here, so each displayed curve uses the trees for which that subject was out of bag. We select six subjects only to keep the lines readable; the extracted object contains all 500. ### Hazard is local ```{r hazard-plot} case_ids <- c(1, 2, 3, 11, 14, 22) plot(rhf_curves, idx = case_ids) ``` Hazard is a local event rate. A high point says that events are occurring at a high instantaneous rate among subjects still at risk around that time. It is not the probability that a subject has the event at that time, so neither its height nor its scale should be read as a percentage. RHF estimates this local quantity directly [@Ishwaran:RHF:2026]. ### Why do the hazard curves stop early? Each line ends where that subject's follow-up ends. This is upstream behavior: from version 2.0.0, `randomForestRHF::rhf()` returns a pointwise hazard only where a grid point falls inside one of the subject's supplied `(start, stop]` records, and `NA` everywhere else, both in the gaps between records and after the final stop time. Earlier versions carried an estimate across the whole grid. `gg_rhf()` passes those `NA` values through unchanged, and `plot()` drops them before drawing, so a curve covers the range the forest was asked about and no more. Cumulative hazard is masked too, but on a different rule, and version 2.0.3 changed it. It accumulates the exact overlap between the working grid and the supplied records, so it stays flat through a gap between two records rather than going missing there. After a subject's final stop, though, there is no more overlap to accumulate, and from 2.0.3 the value is `NA` instead of the last level carried forward. So the cumulative-hazard panel below ends each curve where the hazard panel does, and the difference between the two masks only shows up on a fit with time-dependent covariates, where the gaps are real. If you summarize either column yourself, remember to pass `na.rm = TRUE`. ### Cumulative hazard adds the local rates The same `gg_rhf` object carries `chf`. Set the public `hazard.only` argument to `FALSE` to draw it: ```{r cumulative-hazard-plot} plot(rhf_curves, idx = case_ids, hazard.only = FALSE) ``` The cumulative hazard function (CHF) adds hazard over follow-up. It summarizes accumulated event pressure, but it is unbounded and is not an event probability. Its scale answers "how much hazard has accumulated?", while the hazard plot answers "where is the event rate high right now?" There is one wrinkle in these time-dependent curves. `randomForestRHF` stitches record-specific ensemble estimates as each subject's active record changes. A new record's estimated CHF can start below the preceding record's value, so a stitched case curve may contain a downward step. That step does not mean a negative hazard or that past event pressure disappeared. It marks the upstream model's transition to a different active-record estimate [@Ishwaran:RHF:software:2026]. ## How does discrimination change over time? Survival AUC needs a definition of a case and a control at each time. The two definitions below answer different questions, and each uses the marker that matches its question. The saved AUC calculations are supplied to `gg_auct()`; only the extraction and plotting run during this render. Each call below still names its own `method`, which the saved fit already fixes, so that dropping `auct_fit` and letting `gg_auct()` compute leaves the chunk asking for the same estimand instead of falling back to the `"cumulative"` default. ### Cumulative/dynamic AUC uses cumulative hazard Suppose the question is whether a subject who has experienced the event by a horizon ranks above a subject who remains event-free at that horizon. This is the cumulative/dynamic target. The corresponding marker is cumulative hazard, named `"chf"` in `gg_auct()` and `"cumhaz"` in the retained upstream object. ```{r auct-cumulative} auct_cumulative <- gg_auct( bundle$fit, marker = "chf", method = "cumulative", auct_fit = bundle$auct_cumulative ) plot(auct_cumulative) ``` The saved curve has a Uno iAUC of `r sprintf("%.3f", attr(auct_cumulative, "iauc")$uno)`. Its finite AUC values range from `r sprintf("%.3f", min(auct_cumulative$auc, na.rm = TRUE))` to `r sprintf("%.3f", max(auct_cumulative$auc, na.rm = TRUE))` over this fitted time grid. AUC is a ranking probability on a 0 to 1 scale, with 0.5 shown as the chance reference. The retained calculation did not use a bootstrap, so the plot has no confidence ribbon. This curve needs `randomForestRHF` 2.0.3 or newer, which is why the package asks for that version. Before 2.0.3, `auct.rhf()` could put the cumulative/dynamic curve below the chance line on data the forest fits well. That was upstream rather than in `gg_auct()`, which passes the values through unchanged either way. ### Incident/dynamic AUC uses hazard A different question asks whether a subject who experiences the event near time $t$ ranks above a subject in the relevant risk set at $t$. This is the incident/dynamic target. It uses the local hazard marker, named `"haz"` in `gg_auct()` and `"hazard"` upstream. ```{r auct-incident} auct_incident <- gg_auct( bundle$fit, marker = "haz", method = "incident", riskset = "subject", auct_fit = bundle$auct_incident ) plot(auct_incident) ``` Here the Uno iAUC is `r sprintf("%.3f", attr(auct_incident, "iauc")$uno)`. The curve is more variable than the cumulative/dynamic curve, with finite AUC values from `r sprintf("%.3f", min(auct_incident$auc, na.rm = TRUE))` to `r sprintf("%.3f", max(auct_incident$auc, na.rm = TRUE))`. These values do not show that one marker is a better version of the other. Cumulative/dynamic AUC ranks accumulated risk through a horizon, while incident/dynamic AUC ranks local failures within a risk set. They estimate different targets. ## Which variables matter, and when? The retained `importance.rhf` object is the starting point here. Its expensive upstream work was performed once: `varpro.cache.rhf()` collected the rule and near-miss information, then `importance.rhf()` evaluated five selected time windows. Supplying that result avoids repeating either calculation. ```{r rhf-importance, fig.width=8, fig.height=5.25} rhf_priority <- gg_rhf_importance( bundle$fit, importance_fit = bundle$importance ) plot(rhf_priority) + ggplot2::theme(plot.margin = ggplot2::margin(l = 80)) ``` RHF priority is a time-local rule-release contrast. Within a window, it compares the log integrated-hazard working response for a rule with its near-miss set and asks how much that fitted response changes when rules involving a variable are released. Larger nonnegative values in this supplied object mean a larger contrast. They do not say that a covariate raises or lowers hazard, and they are not z-scores, p-values, or automatic variable selection thresholds. The dot size and color show the same priority value, with a 99th-percentile cap applied to the display only. The extracted values do not change. In this simulation, `xtd` has the largest early-window value, about `r sprintf("%.2f", max(rhf_priority$priority, na.rm = TRUE))`, while `x.1` has the largest value in several later windows. The number at risk falls from `r max(rhf_priority$n_risk)` to `r min(rhf_priority$n_risk)` across the five displayed windows, so the late-window comparisons use less information. One late `x.4` result is unavailable in the upstream object and appears without a point. For work you plan to revisit, retain both the upstream cache and the `importance.rhf` result; supply the latter to `gg_rhf_importance()` for each display. ## Which tree size did the upstream searches select? The two retained `tune.treesize.rhf` objects contain the evaluated paths and their selected sizes. `gg_tune_rhf()` preserves the upstream evaluation order, so the connecting line follows the order in which candidates were evaluated, not a newly sorted sequence. The emphasized point marks the size selected by the upstream search. ### OOB risk tuning ```{r tune-risk} risk_tuning <- gg_tune_rhf(bundle$tune_risk) risk_tuning plot(risk_tuning) ``` The risk search evaluated tree sizes `r paste(risk_tuning$treesize, collapse = ", ")` in that order and selected size `r risk_tuning$treesize[risk_tuning$selected]`, where the saved mean OOB risk was smallest. The retained risk path has no standard-error field. `gg_tune_rhf()` sets its `se` values to missing and does not fabricate uncertainty, so a risk plot has no ribbon. ### OOB iAUC tuning ```{r tune-iauc} iauc_tuning <- gg_tune_rhf(bundle$tune_iauc) iauc_tuning plot(iauc_tuning) ``` The iAUC search evaluated tree sizes `r paste(iauc_tuning$treesize, collapse = ", ")` and selected size `r iauc_tuning$treesize[iauc_tuning$selected]`, where OOB iAUC was largest. An iAUC ribbon is drawn only when finite upstream standard errors are present. This saved path has `r sum(is.finite(iauc_tuning$se))` finite standard errors, so its plot has no ribbon. Retaining each tuning result preserves the actual candidate path, search order, and selected point without rerunning the fits. ## What does each extractor return? All four families keep extraction separate from plotting. The intermediate object is a data frame with an added class, so you can inspect it before calling `plot()`. ```{r extractor-support} support <- data.frame( extractor = c( "gg_rhf()", "gg_auct()", "gg_rhf_importance()", "gg_tune_rhf()" ), upstream_input = c( "rhf fit", "rhf fit plus an optional auct.rhf object", "rhf fit plus an optional importance.rhf object", "tune.treesize.rhf object" ), returned_class = c( "gg_rhf, data.frame", "gg_auct, data.frame", "gg_rhf_importance, data.frame", "gg_tune_rhf, data.frame" ), scale = c( paste( "Hazard is a local event rate; cumulative hazard is accumulated", "event pressure. Neither is a probability." ), paste( "AUC(t) and iAUC are ranking measures from 0 to 1; interpretation", "depends on the case definition and marker." ), paste( "Time-local rule-release priority; larger values mean a larger", "working-response contrast, not statistical significance." ), paste( "OOB risk is minimized or OOB iAUC is maximized; the selected", "tree size is marked." ) ), check.names = FALSE ) knitr::kable( support, col.names = c( "Extractor", "Upstream input class or object", "Returned class", "Output scale and interpretation" ) ) ``` ## How can this analysis be reproduced? The bundle records the public simulation, seed, calculation settings, and package versions used for this render. These are the saved versions: ```{r bundle-versions} version_table <- data.frame( package = names(bundle$versions), version = unname(bundle$versions), row.names = NULL ) knitr::kable( version_table, col.names = c("Software", "Saved version") ) ``` Run `Rscript vignettes/precompute_rhf.R` from the package root to regenerate all ten bundle components. That script repeats the sole longitudinal simulation and uses `bundle$seed` and the settings shown earlier. The exact numerical result can still change when the recorded package versions change. For routine use, retain the fitted forest and each expensive upstream result: the two `auct.rhf` calculations, the importance cache and result, and both tuning objects. Supplying those objects is the default workflow for a report you expect to render more than once. It separates a long model calculation from a short, repeatable extraction and plot. ## Further reading The RHF method, including its counting-process construction and hazard-based forest, is described by Ishwaran, Hsich, Kogalur, and Lee [@Ishwaran:RHF:2026]. The `randomForestRHF` software reference documents the R implementation and the upstream objects used throughout this article [@Ishwaran:RHF:software:2026].