Package {randomForestRHF}

Time-varying AUC(t) and iAUC for Random Hazard Forests (TDC)

Description

Compute the time-varying AUC(t) and its time-aggregated summaries for Random Hazard Forests (RHF) with time-dependent covariates (TDC).

Usage

auct.rhf(
  object,
  marker = c("cumhaz", "hazard", "chf", "haz"),
  method = c("cumulative", "incident"),
  tau = NULL,
  riskset = c("subject", "record"),
  min.controls = 25,
  nfrac.controls = 0.10,
  min.cases = 1,
  g.floor = 0.10,
  g.floor.q = NULL,
  power = 2,
  ydata = NULL,
  winsor.q = NULL,
  eps = 1e-12,
  bootstrap.rep = 0L,
  bootstrap.refit = FALSE,
  bootstrap.conf = 0.95,
  bootstrap.seed = NULL,
  verbose = TRUE
)

## S3 method for class 'auct.rhf'
print(x, digits = 4, max.rows = 8, ...)
## S3 method for class 'auct.rhf'
plot(x, bass = 10, xlab = "Time", ylab = NULL,
         main = NULL, ylim = NULL, pch = 16, alpha = .05, ...)

Arguments

Details

What is estimated. At each evaluation time t, the function computes an AUC for the chosen marker (either cumulative hazard or hazard), using the specified case/control definition:

• Incident/dynamic (method = "incident"): cases fail at t; controls are at risk at t.

• Cumulative/dynamic (method = "cumulative"): cases have failed by t; controls are event-free at t.

The per-time AUC(t) values are combined in two ways:

• iAUC.uno: an inverse-probability-of-censoring weighted (IPCW) average of AUC(t) over time, using a global KM estimate G(t) and a weight proportional to 1 / G(t)^{power} (Uno-style time weighting). Stabilization parameters include g.floor, g.floor.q, and winsor.q.

• iAUC.std: a simple time-standardized mean of AUC(t) over the evaluation grid (trapezoidal rule divided by span). This quantity is sensitive to the time window and the shape of AUC(t) across time.

Markers. For method = "cumulative", using marker = "cumhaz" matches standard time-dependent AUCs. Using marker = "hazard" with "cumulative" often yields smaller AUC(t), as it targets instantaneous risk rather than cumulative risk.

Risk-set choice. When method = "incident" there are two ways to define controls at time t:

• riskset = "subject": at risk when entry < t <= Tstop (one interval per subject). Fast and appropriate when subject rows tile continuous follow-up (no gaps).

• riskset = "record": at risk when a counting-process row satisfies start < t <= stop. Exact under gaps but may yield fewer controls at sparse times.

If rows tile the entire follow-up, the two definitions coincide (up to ties). With gaps, "subject" can over-include controls (lower variance but a conceptual mismatch), while "record" is exact but can be sparser. For method = "cumulative" the risk-set choice is ignored.

Bootstrap. When bootstrap.rep > 0, two modes are available:

• Plug-in (bootstrap.refit = FALSE): subject-level pairs bootstrap holding the fitted marker matrix fixed. Weights the within-time AUC(t) using resampled multiplicities and recomputes Uno's time weights via a weighted KM.

• Refit (bootstrap.refit = TRUE): for each replicate, resample complete subject histories with replacement, retaining repeated draws as separate bootstrap subjects, refit the RHF using the original tuning parameters and saved hazard/CHF configuration, and re-evaluate AUC(t). The returned per-time SEs are matched to the original evaluation grid. For predict objects, the refit is performed using the data carried in the object; predictions are then recomputed for those data. This is computationally heavier but captures variability from model fitting. If every refit fails, the function reports the first captured failure rather than silently returning an all-NA bootstrap result.

The plot method draws a base-R shaded band for AUC(t) when bootstrap SEs are present.

Value

An object of class "auct.rhf" with elements:

• AUC.by.time: data frame with columns time, AUC, n.cases, n.ctrl, G, W.

• iAUC.uno: IPCW (Uno-style) time-averaged AUC.

• iAUC.std: time-standardized mean of AUC(t).

• marker, method, riskset, power, g.floor, g.floor.q, winsor.q, times: metadata.

• diag.riskset: Quick diagnostic comparing control counts under riskset = "subject" vs "record" on a small subset of times (list with times, n.ctrl.subject, n.ctrl.record, n.diff, prop.times.different).

• boot: present when bootstrap.rep > 0 with AUC.se (per-time SE matched to the reporting grid), iAUC.uno.se, iAUC.std.se, conf.level, AUC.lower, AUC.upper, rep, and mode ("plug-in" or "refit").

Generic methods print and plot are provided for compact summaries and visualization.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

Heagerty, P. J., Lumley, T., and Pepe, M. S. (2000). Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics, 56(2), 337–344.

Heagerty, P. J., and Zheng, Y. (2005). Survival model predictive accuracy and ROC curves. Biometrics, 61(1), 92–105.

Uno, H., Tian, L., Cai, T., Kohane, I. S., Wei, L.-J. (2013). A unified inference procedure for a class of measures to assess improvement in risk prediction systems with survival data. Statistics in Medicine, 32(14), 2430–2442.

See Also

rhf, predict.rhf

Examples

## ------------------------------------------------------------
## Peak VO2 example
## ------------------------------------------------------------

data(peakVO2, package = "randomForestSRC")
d <- convert.counting(Surv(ttodead, died) ~ ., peakVO2)
f <- "Surv(id, start, stop, event) ~ ."

o <- rhf(f, d, ntree = 25, nodesize = 5)

## AUC(t) with cumulative/dynamic definition and cumhaz marker
a.chf <- auct.rhf(o, marker = "cumhaz", method = "cumulative")

## AUC(t) with incident/dynamic definition and hazard marker
a.haz <- auct.rhf(o, marker = "hazard", method = "incident")

print(a.chf)
print(a.haz)

oldpar <- par(mfrow = c(1, 2))
plot(a.chf, main = "AUC(t): cumulative + cumhaz")
plot(a.haz, main = "AUC(t): incident + hazard")
par(oldpar)

## ------------------------------------------------------------
##  TDC illustration with training/testing
## ------------------------------------------------------------

trn <- hazard.simulation(1)$dta
tst <- hazard.simulation(1)$dta
f <- "Surv(id, start, stop, event) ~ ."
o <- rhf(f, trn, ntree = 25)
p <- predict(o, tst)
a.trn <- auct.rhf(o)
a.tst <- auct.rhf(p)

oldpar <- par(mfrow = c(1, 2))
plot(a.trn, main = "AUC(t): chf marker, train", ylim = c(0.5,1))
plot(a.tst, main = "AUC(t): chf marker, test",  ylim = c(0.5,1))
par(oldpar)

## ------------------------------------------------------------
## Bootstrap SEs and shaded band
## ------------------------------------------------------------
d <- hazard.simulation(1)$dta
f <- "Surv(id, start, stop, event) ~ ."
o <- rhf(f, d)


oldpar <- par(mfrow = c(1, 2))


## plug-in bootstrap
a.bs1 <- auct.rhf(o, marker = "cumhaz", method = "cumulative",
                  bootstrap.rep = 20, bootstrap.seed = 123)
plot(a.bs1, main = "AUC(t) with bootstrap band (plug-in)")

## refit bootstrap (can be slow)
a.bs2 <- auct.rhf(o, marker = "cumhaz", method = "cumulative",
                  bootstrap.rep = 10, bootstrap.refit = TRUE, bootstrap.seed = 7)
plot(a.bs2, main = "AUC(t) with bootstrap band (refit)")


par(oldpar)

Time-localized VarPro importance for random hazard forests

Description

Computes time-localized variable importance for a random hazard forest (RHF) by using VarPro (variable priority) importance and restricting rule and near-miss memberships to pseudo-individuals whose start-stop intervals overlap selected windows from the master time grid time.interest. The working response used for VarPro importance is now taken directly from the RHF object as the logarithm of the upstream integrated hazard exposure, giving a fast localized view of how variable importance evolves over time.

Usage

importance.rhf(o,
 cache = NULL,
 time.index = NULL,
 trim = 0.1,
 sort = TRUE,
 max.rules.tree,
 max.tree,
 eps = 1e-6,
 y.external = NULL,
 verbose = FALSE,
 ...)

varpro.cache.rhf(o,
 max.rules.tree = 150L,
 max.tree = 150L,
 y.external = NULL,
 eps = 1e-6,
 verbose = FALSE)

## S3 method for class 'importance.rhf'
print(x,
 top = 10L,
 rank.by = c("q90", "median", "mean", "max"),
 digits = 4L,
 scientific.threshold = 1e4,
 ...)

## S3 method for class 'importance.rhf'
as.data.frame(x,
  row.names = NULL,
  optional = FALSE,
  format = c("long", "variable_by_time", "time_by_variable"),
  ...)

dotmatrix.importance.rhf(x,
 vars = NULL,
 top.n.union = 15L,
 variable.labels = NULL,
 time.labels = NULL,
 sort.by = c("q90", "sum", "max", "mean", "median", "alphabetical", "cluster", "none"),
 sort.abs = TRUE,
 transform = c("none", "log10"),
 color.by = c("value", "sign", "single", "none"),
 point.color = "steelblue4",
 value.colors = c("grey85", "steelblue4"),
 sign.colors = c("firebrick3", "grey90", "steelblue4"),
 pch = 16,
 cex.range = c(0.6, 3.2),
 size.cap = 0.99,
 color.cap = 0.99,
 alpha = 0.9,
 show.grid = FALSE,
 grid.col = "grey92",
 legend = TRUE,
 legend.args = list(),
 display.note = FALSE,
 xlab = "",
 ylab = "",
 main = "RHF time-localized VarPro importance",
 axis.cex = 0.7,
 var.cex = 0.7,
 time.label.srt = 45,
 save.plot = FALSE,
 out.file = "rhf_time_varpro_dotmatrix.pdf",
 width = 11,
 height = NULL,
 mar = NULL,
 legend.width = 0.7,
 point.args = list(),
 ...)

barplot.importance.rhf(x,
 vars = NULL,
 top.n.union = 15L,
 variable.labels = NULL,
 time.labels = NULL,
 sort.by = c("q90", "sum", "max", "mean", "median", "alphabetical", "cluster", "none"),
 sort.abs = TRUE,
 transform = c("none", "log10"),
 color.by = c("value", "sign", "single", "none"),
 bar.color = "steelblue4",
 value.colors = c("grey85", "steelblue4"),
 sign.colors = c("firebrick3", "grey90", "steelblue4"),
 bar.width = 0.65,
 bar.max.height = NULL,
 size.cap = 0.99,
 color.cap = 0.99,
 alpha = 0.9,
 show.grid = FALSE,
 grid.col = "grey92",
 zero.line = TRUE,
 zero.col = "grey65",
 legend = TRUE,
 legend.args = list(),
 display.note = FALSE,
 xlab = "",
 ylab = "",
 main = "RHF time-localized VarPro importance",
 axis.cex = 0.7,
 var.cex = 0.7,
 time.label.srt = 45,
 save.plot = FALSE,
 out.file = "rhf_time_varpro_barplot.pdf",
 width = 11,
 height = NULL,
 mar = NULL,
 legend.width = 0.7,
 border = NA,
 bar.args = list(),
 ...)

## S3 method for class 'importance.rhf'
plot(x,
 type = c("barplot", "dotmatrix", "lines"),
 vars = NULL,
 variable.labels = NULL,
 top = 10L,
 rank.by = c("q90", "median", "mean", "max"),
 curve = c("step", "line", "lowess"),
 smooth.f = 2/3,
 display.cap = 0.99,
 display.note = FALSE,
 xlab = NULL,
 ylab = NULL,
 lty = 1,
 lwd = 2,
 line.colors = NULL,
 line.args = list(),
 legend = TRUE,
 legend.args = list(),
 ...)

Arguments

Details

Implements time-localized RHF variable priority. The time grid is given by time.interest. For each selected time window, the method uses the counting process start–stop records whose intervals overlap that window. The forest rules and near-miss sets are taken from the fitted RHF object and are reused across time windows.

For RHF objects, variable priority (VarPro) is applied in a regression-style form. The working response is the logarithm of the integrated hazard exposure assigned by RHF over the relevant interval-window overlap. Local priority for a variable is computed by comparing this working response among records satisfying a forest rule with the corresponding response among records in the rule's near-miss set, obtained by releasing that variable from the rule. Temporal localization comes from restricting the calculation to records active in each time window, rather than rebuilding the rules or refitting the forest at each time point.

The helper varpro.cache() stores the minimum information needed for repeated localized importance calculations: a regression-style rule template, window metadata, the working response source, and precomputed window-local rule statistics. During cache construction, raw OOB and complementary memberships are converted into compact per-window rule summaries, so the later window sweep does not need to rescan membership vectors.

The returned importance matrix has variables in rows and selected time windows in columns. Column names correspond to the right endpoints of the selected windows. The long-format table contains the same values together with window metadata such as start, stop, midpoint, number at risk, and number of active rules.

Printing and plotting share a robust strategy. Summaries default to a robust over-time ranking based on the 90th percentile, and the plotting helpers apply optional quantile capping for display only. This prevents rare extreme spikes from flattening curves while preserving the original importance matrix for downstream analyses.

For consistent base-graphics behavior, ... customizes the initial plot frame for every plot type. Layer-specific options belong in point.args, bar.args, or line.args. The single legend.args list controls the custom matrix legend or is passed to graphics::legend() for a line plot. Data-driven coordinates, sizes, heights, and colors remain under control of the plotting helper. Line-plot limits may still be supplied as xlim or ylim through ....

Value

varpro.cache() returns an object containing cached rule memberships, the working response used for importance, start-stop information for pseudo-individuals, time-window metadata, and the rule extraction settings.

importance.rhf() returns a list including:

• importance.matrix: matrix of localized importance values with variables in rows and selected time windows in columns.

• importance.long: long-format data frame containing variable, time, window metadata, and localized importance.

• window.info: data frame describing the analyzed windows, including start, stop, midpoint, n.risk, and n.rules.

• y.source: source of the working response. This is "int.haz.oob", "int.haz.test", or "y.external".

• trim: tuning value used in importance aggregation.

print() returns its input invisibly after displaying a short summary that includes robust over-time summaries for the leading variables.

as.data.frame() returns one of the supported data-frame views.

dotmatrix.importance.rhf() produces a base-R dot-matrix plot and returns plotting metadata invisibly, including the resolved legend.args settings.

barplot.importance.rhf() produces a base-R bar-matrix plot and returns plotting metadata invisibly, including the resolved legend.args settings.

plot() returns invisibly the result of the underlying plotting helper. For type = "lines", the returned display-capped matrix has attributes containing the selected variable codes, resolved labels and time values.

References

Ishwaran H. (2025). Multivariate Statistics: Classical Foundations and Modern Machine Learning. Chapman and Hall.

Ishwaran H., Kogalur U.B., Hsich E.M. and Lee D.K. (2026). Random hazard forests.

Lu, M. and Ishwaran, H. (2024). Model-independent variable selection via the rule-based variable priority. arXiv e-prints, pp.arXiv-2409.

Lu, M. and Ishwaran, H. (2025). Individual variable priority: a model-independent local gradient method for variable importance. Artificial Intelligence Review, 58:407.

Zhou L., Lu M. and Ishwaran H. (2026). Variable priority for unsupervised variable selection. Pattern Recognition, 172:112727.

See Also

rhf

Examples

################################################################
##
## simulation model
##
################################################################

## draw simulation (can be modified)
n <- 400
p <- 10
simid <- 2
d <- hazard.simulation(type = simid, n = n, p = p, nrecords = 4)$dta

## fit a RHF model with weighted mtry (use for high-dimension)
f <- "Surv(id, start, stop, event) ~ ."
o <- rhf(f, d, ntree = 50, nsplit = 5, xvar.wt = xvar.wt.rhf(f, d))
print(o)

## time-localized RHF importance across the full time grid
imp.t <- importance.rhf(o)
print(imp.t)

## extract the variable-by-time matrix
print(head(imp.t$importance.matrix))

oldpar <- par(mfrow=c(1,1))

## bar-matrix importance plot (default)
plot(imp.t)

## companion dot-matrix importance plot
plot(imp.t, type = "dotmatrix")

## suppress matrix-legend titles and enlarge its numeric labels
plot(imp.t, legend.args = list(title = FALSE, cex = 1.15))

## retain titles but control their size separately
plot(imp.t, legend.args = list(cex = 1.05, title.cex = 0.9))

## step-style importance line plot for the top variables
## (ranked by the 90th percentile over time and display-capped at q99)
plot(imp.t, type = "lines", top = 10)

## smoothed importance plot for all variables with display capping
plot(imp.t, type = "lines", curve = "lowess", smooth.f = 0.5, display.cap = 0.95)

## robust ordering and display capping
plot(imp.t, sort.by = "q90", size.cap = 0.99, color.cap = 0.99)

## bar-matrix plot using signed colors, if signed importance values are present
plot(imp.t, color.by = "sign")

par(oldpar)


## reuse a cache for repeated calls on subsets of the time grid
cache <- varpro.cache(o)
imp.t.sub <- importance.rhf(
  o,
  cache = cache,
  time.index = seq(1, length(o$time.interest), by = 5),
  verbose = TRUE
)

## long-format export
print(head(as.data.frame(imp.t.sub)))

Advanced Utility Functions for RHF Workflows

Description

These helper functions support advanced random-hazard-forest workflows, including data conversion, simulation, and predictor weighting.

Usage

convert.counting(f, dta, scale = FALSE)

convert.standard.counting(formula, data,
  scale = FALSE,
  rescale.from.attr = FALSE,
  keep.id = FALSE,
  keep.row_index = FALSE,
  sorted = FALSE,
  id.default = "id",
  eps = 1e-8,
  landmark.time = NULL,
  landmark.use.tminus = TRUE,
  return.type = c("survival", "x"),
  keep.landmark.cols = FALSE)

hazard.simulation(type = 1,
  n = 500, p = 10, nrecords = 7,
  scale = FALSE, ngrid = 1e5, ...)

xvar.wt.rhf(f, 
    d,
    scale = 4,
    rank.by = c("q90", "median", "mean", "max", "sum"),
    rhf.args = list(),
    cache.args = list())

Arguments

Details

convert.counting() converts standard right-censored survival data to the counting-process representation expected by rhf().

convert.standard.counting() converts counting-process style RHF data to a conventional survival-analysis data frame or to a predictor-only frame.

hazard.simulation() generates synthetic counting-process data together with the corresponding theoretical hazard function.

xvar.wt.rhf() computes variable weights that can be supplied to rhf(..., xvar.wt = ...).

Value

The returned value depends on the helper:

• convert.counting(): a data frame with columns id, start, stop, event, followed by the predictor columns from the input data.

• convert.standard.counting(): a data frame in the format requested by return.type. With return.type = "survival", the leading columns are time and event; additional covariate and optional helper columns may also be present. The function may return an empty data frame when no valid rows remain after filtering.

• hazard.simulation(): a list with components dta (simulated counting-process data frame), haz (a function that evaluates the true hazard on a supplied time grid), and scale (the time-scaling factor used in the simulation).

• xvar.wt.rhf(): a named numeric vector of predictor weights for the xvar.wt argument of rhf().

See Also

rhf, predict.rhf

Plot smoothed hazard and cumulative hazard plots from RHF analysis

Description

Plot case specific hazard and cumulative hazard (CHF) from a fitted random hazard forest (RHF) object. Hazards are smoothed with stats::supsmu. Optional scaling places hazards on the same scale as the CHF for easier comparison.

Usage

## S3 method for class 'rhf'
plot(x, idx = NULL, scale.hazard = FALSE, 
         ngrid = 30, bass = 0, q = 0.99, grid = FALSE,
         col = NULL, lty = NULL, legend.loc = "topright",
         jitter.factor = 1, lwd = 4,
         hazard.only = TRUE, legend.show = TRUE, ...)

Arguments

Details

The function displays out of bag hazard and CHF estimates on the grid x$time.interest. Hazards are smoothed using supsmu after trimming large values for stability. Scaling the hazard allows direct visual comparison with CHF. Plots are drawn to the current device; no model re-estimation is performed.

Value

Invisibly returns NULL. The function is used for its plotting side effects.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

Ishwaran H. and Kogalur U.B. (2007). Random survival forests for R. R News, 7(2): 25–31.

Ishwaran H., Kogalur U.B., Blackstone E.H. and Lauer M.S. (2008). Random survival forests. Annals of Applied Statistics, 2: 841–860.

Lee D.K., Chen N., and Ishwaran H. (2021). Boosted nonparametric hazards with time dependent covariates. Annals of Statistics, 49: 2101–2128.

Ishwaran H., Lee D.K. and Hsich E.M. (2025). Random hazard forests.

See Also

hazard.simulation, predict.rhf, rhf, smoothed.hazard.rhf

Examples

## ------------------------------------------------------------
##  time static peakVO2
## ------------------------------------------------------------
data(peakVO2, package = "randomForestSRC")
d <- convert.counting(Surv(ttodead, died)~., peakVO2)
f <- "Surv(id, start, stop, event) ~ ."

set.seed(1)
o <- rhf(f, d)

oldpar <- par(mfrow=c(1,1))

## plot selected cases
ids <- o$ensemble.id[1:3]
plot(o, idx = ids)

## custom axis label
plot(o, idx = ids, ylab = "my hazard")

## hazard only
plot(o, idx = ids)

## scaled hazard with CHF
plot(o, idx = ids, scale.hazard = TRUE, hazard.only = FALSE)

## auxiliary grid with smoother control
plot(o, idx = ids, grid = TRUE, ngrid = 60, bass = 2, hazard.only = FALSE)

## multiple cases, no legend
plot(o, idx = o$ensemble.id[1:10], lwd = 0, legend.show = FALSE)

## lowess median smoothed hazard
s <- smoothed.hazard(o)  ## default method="median.loess"
plot(s, idx = o$ensemble.id[1:10])
plot(s, idx = o$ensemble.id[1:10], lwd = 0)


par(oldpar)

## ------------------------------------------------------------
##  complex simulated time dependent covariate hazards
## ------------------------------------------------------------
f <- "Surv(id, start, stop, event) ~ ."
sim2 <- hazard.simulation(2)$dta
o2 <- rhf(f, sim2)

oldpar <- par(mfrow=c(1,1))
plot(o2, 1)
par(oldpar)

## ------------------------------------------------------------
##  same complex simulation, but using smoothed hazard
## ------------------------------------------------------------
f <- "Surv(id, start, stop, event) ~ ."
sim2 <- hazard.simulation(2)$dta
o2 <- rhf(f, sim2)

s2.loess <- smoothed.hazard(o2, method="loess")
s2.med.loess <- smoothed.hazard(o2, method="median.loess")

oldpar <- par(mfrow=c(1,1))
plot(o2, 1)
plot(s2.loess, 1)
plot(s2.med.loess, 1)
plot(s2.med.loess, 1:10, lwd = 0)
par(oldpar)

Prediction on Test Data for Random Hazard Forests

Description

Obtain predicted values on test data using a trained random hazard forests.

Usage

## S3 method for class 'rhf'
predict(object, newdata,  get.tree = NULL,
     block.size = 10, membership = TRUE, adaptive = TRUE,
     seed = NULL, do.trace = FALSE,...)

Arguments

Details

Returns the predicted values for a random hazard forests.

Prediction inherits the hazard calculation configuration saved with the trained forest. In ordinary test-data prediction there is no OOB objective for reselecting the trim value, so the grow-time selected index is used. In restore mode, the stored forest can reconstruct OOB quantities and may reselect the trim value over the active candidate grid. Setting adaptive = FALSE changes that grid to the single value 0.05, unless coe.trim is supplied through ....

Value

An object of class c("rhf", "predict", family). The returned list contains the fitted forest together with prediction summaries on the evaluation grid time.interest. Important components include:

• hazard.test, chf.test, risk.test, and int.haz.test: test-set hazard, cumulative hazard, risk, and integrated-hazard summaries when newdata is supplied.

• hazard.oob, chf.oob, risk.oob, and int.haz.oob: out-of-bag summaries for the training data.

• hazard.inbag, chf.inbag, risk.inbag, and int.haz.inbag: in-bag summaries when available.

• id, yvar, and xvar: identifiers and processed outcome/predictor data used by the returned prediction object.

• pseudo.membership and inbag: terminal-node membership and inbag information when membership = TRUE.

• forest: the fitted forest object used to generate the predictions.

If newdata is omitted, the function restores predictions for the original training data using the stored forest and returns the same class of object.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

Ishwaran H. and Kogalur U.B. (2007). Random survival forests for R, Rnews, 7(2):25-31.

Ishwaran H., Kogalur U.B., Blackstone E.H. and Lauer M.S. (2008). Random survival forests, Ann. App. Statist., 2:841-860.

Lee, D.K. and Chen N. and Ishwaran H (2021). Boosted nonparametric hazards with time-dependent covariates. Annals of Statistics, 49: 2101-2128.

Ishwaran H. (2025). Multivariate Statistics: Classical Foundations and Modern Machine Learning. Chapman and Hall.

Ishwaran H., Kogalur U.B., Hsich E.M. and Lee D.K. (2026). Random hazard forests.

See Also

rhf

Examples

## ------------------------------------------------------------
## canonical train/test example (synthetic data)
## ------------------------------------------------------------

simID <- 1
trn <- hazard.simulation(simID)$dta
tst <- hazard.simulation(simID)$dta
f <- "Surv(id, start, stop, event) ~ ."

## training
o <- rhf(f, trn, ntree = 3)
print(o)

## testing
p <- predict(o, tst)
print(p)

## fixed trim protocol
p.fixed <- predict(o, tst, adaptive = FALSE)



## ------------------------------------------------------------
##  pbc: train/test example
## ------------------------------------------------------------

library("randomForestSRC")
data(pbc, package = "randomForestSRC")
pbc.raw <- na.omit(pbc)
trn <- sample(1:nrow(pbc.raw), size=nrow(pbc.raw) * .75, replace = FALSE)

d.trn <- convert.counting(Surv(days, status) ~ ., pbc.raw[trn,])
d.tst <- convert.counting(Surv(days, status) ~ ., pbc.raw[-trn,])
f <- "Surv(id, start, stop, event) ~ ."

## train/predict
o <- rhf(f, d.trn)
print(predict(o, d.tst))

## restore the forest
print(predict(o))

Print Output from a Random Hazard Forest Analysis

Description

Print a compact summary of a Random Hazard Forest analysis. This is the default print method for objects produced by grow or predict.

Usage

## S3 method for class 'rhf'
print(x, digits = 3, label.width = 34, ...)

Arguments

Details

The printed summary includes basic sample information, forest size characteristics, event counts when available, and risk metrics (in-bag, OOB, or test risk depending on how the object was created).

Value

Returns the object x invisibly.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

Random Hazard Forests

Description

Random Hazard Forest (RHF) is a tree-ensemble survival method that estimates case-specific hazards and cumulative hazard functions on a working time grid and supports time-dependent covariates using counting-process (start/stop) format.

Usage

rhf(formula,
      data,
      ntree = 500,
      nsplit = 10,
      treesize = NULL,
      nodesize = NULL,
      block.size = 10,
      bootstrap = c("by.root", "none", "by.user"),
      samptype = c("swor", "swr"),
      samp = NULL,
      case.wt = NULL,
      membership = TRUE,
      sampsize = if (samptype == "swor") function(x){x * .632} else function(x){x},
      xvar.wt = NULL,
      ntime = 50,
      min.events.per.gap = 10,
      adaptive = TRUE,
      seed = NULL,
      do.trace = FALSE,
      ...)

Arguments

Details

Overview

rhf() is the main fitting function for Random Hazard Forests (RHF). It grows an ensemble of hazard trees from survival data represented in counting-process, start–stop form, Surv(id, start, stop, event). Multiple rows for the same id describe an evolving covariate path. Predictors may include time-static variables, repeated across rows, and time-dependent covariates that change from row to row.

RHF estimates hazards pathwise on the time grid time.interest. For a supplied predictable covariate path, each tree routes the current covariate state through its partition over time and returns the corresponding terminal-node hazard. Averaging across trees yields a case-specific hazard trajectory. A fitted RHF object stores in-sample ensemble estimates of the hazard and cumulative hazard on time.interest. Predictions for new data, or from restored fitted objects, are obtained with predict.rhf().

The default hazard aggregation uses an adaptive trimming protocol: RHF evaluates a grid of candidate trimming fractions and retains the candidate selected by OOB risk. Setting adaptive = FALSE uses the fixed trimming fraction coe.trim = 0.05 instead. This can be useful when another tuning procedure should hold the hazard aggregation protocol fixed; however the value 0.05 is arbitrary and is not guaranteed to give good peformance in all settings.

Predictable covariate paths

The start–stop rows must represent a predictable covariate process. That is, the covariate values used to evaluate the hazard at a given time must be determined by information available before that time. Measurements or treatments recorded later must not be used to construct covariates for an earlier prediction time. rhf() assumes that this no-lookahead convention has already been encoded in the input data.

Required data format

The input data must contain the following four columns. Column names must agree with those specified in the formula; see the examples.

id

A unique integer identifier for each individual. The identifier is repeated across all rows belonging to that individual.

start

The start time of the row. All times must be scaled to [0, 1] and must use the same scale as time.interest.

stop

The stop time of the row. Each row must satisfy stop > start.

event

A binary event indicator for the row (0 = no event, 1 = event). The current implementation supports single-event survival data, so at most one event may be recorded for each individual.

Tree construction

Each tree is grown by a best-first search guided by reduction in empirical hazard risk. Candidate splits are evaluated using a nonparametric hazard likelihood for predictable covariate processes. For computational efficiency, tree induction uses a working model that assigns a time-constant hazard to each current node. Under this restriction, candidate splits can be scored using node-level at-risk exposure and event totals.

The node-constant model is used only to choose the partition. Once the partition is fixed, RHF removes this restriction and estimates a time-varying step-function hazard in each terminal node on time.interest. Both static and time-dependent covariates are eligible for splitting. The option xvar.wt controls how often individual variables are considered, and can be used, for example, to give additional weight to time-dependent covariates.

In-sample fitted trajectories

For the training data, RHF assembles case-specific fitted curves on time.interest using an active-record rule. Each start–stop row supplies the routing covariate state for the portion of the grid on which that row is active. Rows belonging to the same individual are then stitched together to form the fitted pathwise hazard trajectory. At a grid point t, the returned case-specific hazard is defined only when t lies in at least one of the individual's supplied (\mathit{start},\mathit{stop}] intervals. It is returned as NA in gaps and after the final supplied stop time; covariate values are not carried forward outside the observed path.

The corresponding cumulative hazard is calculated using the exact overlap between the working time-grid intervals and the supplied start–stop path. It therefore remains flat through gaps and after follow-up. An observed portion of a final grid interval still contributes to cumulative hazard even when the pointwise hazard at that grid interval's right endpoint lies outside the supplied path. These output conventions do not alter the likelihood contribution of the training data.

Data conversion

The helper function convert.counting() converts conventional survival data to counting-process form. See the examples for additional guidance on data preparation.

Value

An object of class "rhf" containing the fitted forest and related results. Components include (among others) the working time grid time.interest and case-specific ensemble estimates of the hazard and cumulative hazard. When bootstrapping is used, these are typically returned as hazard.oob/chf.oob (out-of-bag) and hazard.inbag/chf.inbag (in-bag). Use predict.rhf to obtain estimates for new data.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

Ishwaran H. and Kogalur U.B. (2007). Random survival forests for R, Rnews, 7(2):25-31.

Ishwaran H., Kogalur U.B., Blackstone E.H. and Lauer M.S. (2008). Random survival forests, Ann. App. Statist., 2:841-860.

Lee, D.K. and Chen N. and Ishwaran H (2021). Boosted nonparametric hazards with time-dependent covariates. Annals of Statistics, 49: 2101-2128.

Ishwaran H. (2025). Multivariate Statistics: Classical Foundations and Modern Machine Learning. Chapman and Hall.

Ishwaran H., Kogalur U.B., Hsich E.M. and Lee D.K. (2026). Random hazard forests.

See Also

auct.rhf, importance.rhf, plot.rhf, predict.rhf, smoothed.hazard.rhf, tune.treesize.rhf

Examples

  
## ------------------------------------------------------------
##  time-static pbc: parameters set for fast CRAN run
## ------------------------------------------------------------

## load the data
data(pbc, package = "randomForestSRC")
## convert the data to counting process
d <- convert.counting(Surv(days, status) ~ ., na.omit(pbc))
## set the formula
f <- "Surv(id, start, stop, event) ~ ."

## rhf call
print((o <- rhf(f, d, ntree = 3)))

## smooth hazard estimator for specific cases
plot(o, idx=c(1,5,10))
plot(o, idx=c(1,5,10), hazard.only=TRUE)
plot(o, idx=c(1,5,10), hazard.only=TRUE, lwd=0)
plot(o, idx=1:10, lwd=0, hazard.only=TRUE, legend.show=FALSE)



## ------------------------------------------------------------
##  time-static pbc: compare RHF to RSF
## ------------------------------------------------------------

data(pbc, package = "randomForestSRC")
d <- convert.counting(Surv(days, status) ~ ., na.omit(pbc))

## first we run rhf
o.rhf <- rhf("Surv(id, start, stop, event) ~ .", d)
h <- o.rhf$hazard.oob
time <- o.rhf$time.interest
H <- t(o.rhf$chf.oob)
S.rhf <- exp(-H)

## next we run rsf, using the same time.interest grid
o.rsf <- randomForestSRC::rfsrc(Surv(days, status) ~ ., pbc, ntime = time)
S.rsf <- t(o.rsf$survival.oob)

## graphical parameters
bass <- 1
oldpar <- par(mfrow=c(2,2))

## plot survival results 
matplot(time,S.rhf,type="l",col=gray(.6))
matplot(time,S.rsf,type="l",col=gray(.6))
boxplot(S.rhf-S.rsf,ylab="S.rhf-S.rsf",xaxt="n",outline=FALSE,range=1e-10)
abline(h=0,lwd=3,col=2)

## plot chf and hazard results
finite <- is.finite(h)
matplot(time,t(h),ylab="hazard",pch=16,type="n",
        ylim=c(0,quantile(h[finite],.99)))
nO <- lapply(which(rowSums(finite) >= 2L), function(i) {
  lines(supsmu(time[finite[i,]],h[i,finite[i,]],bass=bass),
        type="b",pch=16,col=gray(.6),cex=.5)
})

par(oldpar)

## ------------------------------------------------------------
##  TDC illustration (using built in hazard simulation function)
## ------------------------------------------------------------

d1 <- hazard.simulation(1)$dta
d2 <- hazard.simulation(2)$dta
d3 <- hazard.simulation(3)$dta
f <- "Surv(id, start, stop, event) ~ ."

o1 <- rhf(f, d1)
o2 <- rhf(f, d2)
o3 <- rhf(f, d3)

plot(o1,idx=4)
plot(o2,idx=1:3)
plot(o3,idx=2)


## ------------------------------------------------------------
##  peakVO2: demonstrates time-varying auc
## ------------------------------------------------------------

data(peakVO2, package = "randomForestSRC")
d <- convert.counting(Surv(ttodead, died)~., peakVO2)
f <- "Surv(id, start, stop, event) ~ ."

## build the forest
o1 <- rhf(f, d, treesize=10)
o2 <- rhf(f, d, treesize=30)

## acquire auc-t
a1.chf <- auct.rhf(o1) ## same as auct.rhf(o1, marker = "chf")
a1.haz <- auct.rhf(o1, marker = "haz")
a2.chf <- auct.rhf(o2) ## same as auct.rhf(o2, marker = "chf")
a2.haz <- auct.rhf(o2, marker = "haz")

## print auc-t
print(a1.chf)
print(a2.chf)
print(a1.haz)
print(a2.haz)

## plot auc-t
oldpar <- par(mfrow=c(2,2))
plot(a1.chf, main="treesize 1")
plot(a1.haz, main="treesize 1")
plot(a2.chf, main="treesize 15")
plot(a2.haz, main="treesize 15")
par(oldpar)

## ------------------------------------------------------------
##  more detailed example of time-varying performance 
## ------------------------------------------------------------

d <- hazard.simulation(1)$dta
f <- "Surv(id, start, stop, event) ~ ."
treesize <- c(10, 30, 100)

oldpar <- par(mfrow=c(3,2))
perf <- lapply(treesize, function(tz) {
  o <- rhf(f, d, treesize=tz)
  print(o)
  rO <- list(chf = auct.rhf(o, marker="chf"),
             haz = auct.rhf(o, marker="haz"))
  plot(rO$chf, main=paste("chf marker, treesize=", tz))
  plot(rO$haz, main=paste("hazard marker, treesize=", tz))
  rO
})
par(oldpar)

## ------------------------------------------------------------
##  example illustrating tuning the tree size
##  using function 'rhf.tune.treesize'
## ------------------------------------------------------------

d <- hazard.simulation(1)$dta
f <- "Surv(id, start, stop, event) ~ ."

tune <- tune.rhf(f, d)  ## same as rhf.tune.treesize(f,d)
oldpar <- par(mfrow=c(1,1))
plot(tune)
par(oldpar)

## ------------------------------------------------------------
##  example illustrating guided feature selection
##  using built in 'xvar.wt.rhf' helper
## ------------------------------------------------------------

d <- hazard.simulation(1)$dta
f <- "Surv(id, start, stop, event) ~ ."

o <- rhf(f,d)
o.gfs <- rhf(f,d, xvar.wt = xvar.wt.rhf(f, d))

print(o)
print(o.gfs)

Show the NEWS file

Description

Show the NEWS file of the randomForestRHF package.

Usage

rhf.news(...)

Arguments

Value

None.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

Smoothed hazard and cumulative hazard estimates from a random hazard forest

Description

Compute smoothed case-specific hazard curves (and corresponding cumulative hazards) from a Random Hazard Forest (RHF) object. The returned object is a lightweight list designed to be directly usable by plot.rhf.

Usage

smoothed.hazard.rhf(
  o,
  method = c("median.loess", "loess"),
  oob = TRUE,
  span = 0.35,
  degree = 1,
  family = NULL,
  eps = 1e-12,
  suppress.warnings = TRUE,
  trace = FALSE
)

Arguments

Details

Hazard curves are computed/smoothed on the grid o$time.interest. Smoothing is performed on the log-hazard scale: log(hazard + eps), and then back-transformed to the hazard scale.

method="loess"

Applies loess smoothing to an already-constructed hazard matrix stored in o (e.g. o$hazard.oob and/or o$hazard.inbag), independently for each subject.

method="median.loess"

For each subject and time bin, computes the median across trees of the terminal-node (leaf) log-hazard log(h + eps) (using a native C routine), then applies loess smoothing to the resulting median log-hazard curve and back-transforms. This method requires the fitted object to retain tree-specific hazard and membership structures (e.g., o$forest$t.hazard and related arrays).

Value

A list containing:

• id, ensemble.id, yvar, time.interest copied from o.

• Depending on target:

  • hazard.oob, chf.oob (out-of-bag), matrices of dimension nSubj x length(time.interest).

  • hazard.inbag, chf.inbag (in-bag), matrices of dimension nSubj x length(time.interest).

  • hazard.test, chf.test (test), matrices of dimension nSubjTest x length(time.interest).

See Also

plot.rhf, predict.rhf, rhf

Examples

## ------------------------------------------------------------
## canonical example using synthetic data
## includes both train/test scenarios
## ------------------------------------------------------------

simID <- 1
trn <- hazard.simulation(simID)$dta
tst <- hazard.simulation(simID)$dta
f <- "Surv(id, start, stop, event) ~ ."

## training/testing
o <- rhf(f, trn)
p <- predict(o, tst)

## default: median(log leaf hazard) + loess
so <- smoothed.hazard(o)
sp <- smoothed.hazard(p)

## the returned object can be passed directly to plot.rhf()
oldpar <- par(mfrow=c(1,1))
plot.rhf(so)
plot.rhf(sp)
par(oldpar)

## ------------------------------------------------------------
##  peak vo2
## ------------------------------------------------------------

data(peakVO2, package = "randomForestSRC")
d <- convert.counting(Surv(ttodead, died)~., peakVO2)
f <- "Surv(id, start, stop, event) ~ ."

## training
o <- rhf(f, d)

oldpar <- par(mfrow=c(1,1))


## median loess
s0 <- smoothed.hazard(o) 
ids <- o$ensemble.id[1:3]
plot.rhf(s0, idx = ids)

## loess smoothing 
s1 <- smoothed.hazard(o, method = "loess")
plot.rhf(s1, idx = ids)

par(oldpar)

Tune tree size for Random Hazard Forests

Description

Tune the treesize for a Random Hazard Forest using either out-of-bag (OOB) empirical risk or an OOB integrated AUC computed by auct.rhf.

Usage

tune.treesize.rhf(
  formula,
  data,
  ntree = 500,
  nsplit = 10,
  nodesize = NULL,
  ntime = 50,
  min.events.per.gap = 10,
  adaptive = TRUE,
  perf = c("risk", "iAUC"),
  auct.args = NULL,
  lower = 5L,
  upper = NULL,
  C = 3,
  method = c("golden", "bisect"),
  max.evals = 20L,
  bracket.tol = 2L,
  seed = NULL,
  verbose = TRUE,
  forest = TRUE,
  ...
)


tune.iAUC.rhf(
  formula,
  data,
  auct.args = NULL,
  adaptive = TRUE,
  ...
)

## S3 method for class 'tune.treesize.rhf'
plot(x, ylab = NULL, main = NULL,
    se.band = TRUE, se.mult = 1, ylim = NULL, ...)

Arguments

Details

The alias tune.rhf is provided for convenience.

Repeatedly fits a hazard forest with different values of treesize and selects the value that optimizes the chosen criterion perf. When perf = "risk", the criterion is the mean OOB risk. When perf = "iAUC", auct.rhf is applied to each fit and the criterion is 1 - iAUC.uno. tune.iAUC is a wrapper around tune.rhf with perf = "iAUC", using default auct.rhf settings when auct.args = NULL.

Because the default RHF fit adaptively selects the trimming fraction by OOB risk, tuning by OOB risk can otherwise combine the choice of treesize with adaptive trimming. Use adaptive = FALSE for a controlled tree-size search that fixes the trimming fraction at coe.trim = 0.05.

Value

An object of class "tune.treesize.rhf" with components

• best.size: optimal treesize;

• best.err: minimal criterion value (mean OOB risk if perf = "risk", or 1 - iAUC.uno if perf = "iAUC");

• path: data frame with evaluated treesize values and corresponding criterion values (and an iAUC column when perf = "iAUC");

• adaptive: logical value indicating whether adaptive trimming was used during tuning;

• forest: RHF fit at best.size when forest = TRUE.

See Also

rhf, auct.rhf

Examples

## ------------------------------------------------------------
## Example: tuning treesize by OOB empirical risk
## ------------------------------------------------------------

set.seed(7)
d <- hazard.simulation(1)$dta
f <- "Surv(id, start, stop, event) ~ ."

## tune.rhf is an alias for tune.treesize.rhf
tuneRisk <- tune.rhf(f, d)  ## same as tune.treesize.rhf(f, d)

## For a controlled tree-size search, keep the trimming fraction fixed.
tuneRisk.fixed <- tune.rhf(f, d, adaptive = FALSE)

oldpar <- par(mfrow = c(1, 1))
plot(tuneRisk)
par(oldpar)

tuneRisk$best.size    # treesize minimizing mean OOB risk
tuneRisk$best.err

## Access the tuned forest and refit AUC diagnostics if desired
best.forest <- tuneRisk$forest
a.risk <- auct.rhf(best.forest)
print(a.risk$iAUC.uno)


## ------------------------------------------------------------
## Example: tuning treesize by iAUC
## ------------------------------------------------------------

set.seed(7)
tuneiAUC <- tune.iAUC(f, d)

oldpar <- par(mfrow = c(1, 1))
plot(tuneiAUC)
par(oldpar)

print(tuneiAUC$best.size)
print(1 - tuneiAUC$best.err)

best.forest.iauc <- tuneiAUC$forest
a.iauc <- auct.rhf(best.forest.iauc)
print(1 - a.iauc$iAUC.uno)


## ------------------------------------------------------------
## Example: tuning treesize by iAUC with bootstrap s.e.
## ------------------------------------------------------------

set.seed(7)
tuneiAUC.boot <- tune.iAUC(f, d, auct.args=list(bootstrap.rep=25))
oldpar <- par(mfrow = c(1, 1))
plot(tuneiAUC.boot)        
par(oldpar)

## ------------------------------------------------------------
## Example: tuning by hazard-based incident AUC
## increase min.cases for stability
## ------------------------------------------------------------

set.seed(7)
tuneHazInc <- tune.iAUC(
  f, d,
  auct.args = list(
    min.cases    = 100,
    marker       = "hazard",
    method       = "incident"
  )
)

oldpar <- par(mfrow = c(1, 1))
plot(tuneHazInc)
par(oldpar)

print(tuneHazInc$best.size)
print(tuneHazInc$best.err)

## Confirm the tuned incident AUC on the selected forest
best.forest.haz <- tuneHazInc$forest
aHazInc <- auct.rhf(
  best.forest.haz,
  min.cases = 100,
  marker  = "hazard",
  method  = "incident"
)
print(aHazInc)
print(1 - aHazInc$iAUC.uno)