rhf.RdRandom Hazard Forest (RHF) is a tree-ensemble survival method that
estimates case-specific hazards and cumulative hazard functions on a
working time grid. It accepts ordinary right-censored data through
Surv(time, event) and time-dependent covariates through
counting-process Surv(id, start, stop, event) input.
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,
...)Model formula specifying the response and predictors.
The response may be ordinary right-censored survival data,
Surv(time, event), or RHF counting-process data,
Surv(id, start, stop, event).
Data frame containing the variables referenced in
formula. Ordinary right-censored input has one row per subject;
counting-process input may have multiple rows per subject.
Number of trees to grow in the ensemble.
Non-negative integer specifying the number of random split points to consider for each variable. When set to 0, all possible split points are evaluated (deterministic splitting), which may be slower. The default is 10.
Specifies the desired tree size, defined as the number of terminal nodes. If unspecified, an internal default is used.
Minimum number of cases required in a terminal node. If not specified, an internal default is used.
Controls the granularity of cumulative error rate
reporting. Setting this to an integer between 1 and ntree
yields the cumulative error every block.size trees. Currently
ignored.
Bootstrap strategy for generating inbag samples. The
default is "by.root", which bootstraps by sampling with or
without replacement (default is without; see samptype). If set
to "none", no bootstrapping is performed (OOB predictions and
errors are then unavailable). The option "by.user" uses a
user-defined bootstrap specified by samp.
Sampling type used when bootstrap =
"by.root". Options are "swor" (sampling without replacement,
default) and "swr" (sampling with replacement).
User-specified bootstrap when bootstrap =
"by.user". Should be an n by ntree array, where each
entry gives the number of times a case appears in the inbag sample for
a given tree.
Non-negative vector of case weights (not required to sum to 1). Higher weights increase the probability of an observation being selected during bootstrapping or subsampling. Using real-valued weights is preferred over integer weights.
Logical flag indicating whether to return terminal node membership and inbag information.
Specifies the size of the bootstrap sample when
bootstrap = "by.root". For sampling without replacement, it
represents the requested subsample size (default is 0.632 times the
sample size). For sampling with replacement, it equals the sample
size. Can be supplied as a numeric value or a function.
Non-negative vector of variable selection probabilities for splitting. Values do not need to sum to 1. Defaults to uniform selection probability.
Controls the working time grid used for all ensemble calculations. Can be:
ntime >= 1 (integer):Request approximately ntime grid
points chosen from the observed event times. When used with
min.events.per.gap, grid points are selected adaptively so that each
interval contains at least min.events.per.gap events (when possible);
therefore the resulting grid may contain fewer than ntime points when
events are sparse (especially in the tail).
A user-supplied set of time points. Each value is aligned (snapped) to the nearest observed event time and duplicates are removed.
0 or NULL:Use all observed (unique) event times.
Minimum number of observed events required in
each time interval when ntime is specified as an integer.
Together with ntime, this provides an automatic event-balanced grid
that allocates more resolution where events are dense and avoids
sparse tail intervals.
Logical value controlling the default trimmed-mean
aggregation protocol. When TRUE, the default, RHF evaluates a
candidate vector of trimming fractions and retains the value selected
by the OOB-risk criterion. When FALSE, the default candidate
vector is replaced by the single value coe.trim = 0.05. This can be
useful when tuning other parameters, such as treesize, without
simultaneously optimizing the trimming fraction, otherwise leave at
the default.
Negative integer setting the random seed for reproducibility.
Time in seconds between progress updates printed to the console.
Additional expert-level controls passed to internal RHF methods.
Overview
rhf() is the main fitting function for Random Hazard Forests
(RHF). Ordinary right-censored data can be supplied directly as
Surv(time, event) with one row per subject. RHF converts this
internally to a one-interval predictable process with start = 0,
stop = time, and time-static predictors. Time-dependent data are
supplied explicitly in counting-process form as
Surv(id, start, stop, event). Multiple rows for the same
id then 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
For counting-process input, 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. For ordinary Surv(time, event)
input, each subject's supplied predictors are treated as time-static
over the interval from zero to the observed time.
Supported response formats
RHF accepts either of the following response forms.
Surv(time, event)Ordinary terminal right-censored data with one row per subject. The
time and event arguments must each be column names in
data. Times must be non-negative, and event must be
coded 0 for censoring and 1 for an event. RHF creates
a private subject identifier and internally represents each row as
one counting-process interval (0, time].
Surv(id, start, stop, event)Counting-process data for time-dependent covariates. The subject
identifier id is repeated across all rows for that subject.
Each row must satisfy stop > start. Intervals may be
contiguous or separated by gaps, but they may not overlap within a
subject. The row-level event increment must be coded
0/1; terminal-event data may contain at most one event per
subject, on that subject's final row.
Times are supplied in their original non-negative units. RHF maps them
to its internal working scale during data processing and returns
time.interest and time-based predictions in the original units.
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
Manual conversion is not required for ordinary right-censored data:
rhf(Surv(time, event) ~ ., data) performs the one-interval
conversion internally. The helper convert.counting() remains
available when an explicit counting-process data frame is desired.
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). The component
input.info records whether the forest was grown from ordinary
right-censored or counting-process input and retains the response names
needed for prediction. Use predict.rhf to obtain estimates
for new data.
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.
## ------------------------------------------------------------
## time-static pbc: parameters set for fast CRAN run
## ------------------------------------------------------------
## load the data
data(pbc, package = "randomForestSRC")
pbc.raw <- na.omit(pbc)
## ordinary right-censored data can be supplied directly
f <- Surv(days, status) ~ .
print((o <- rhf(f, pbc.raw, 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)
# \donttest{
## ------------------------------------------------------------
## time-static pbc: compare RHF to RSF
## ------------------------------------------------------------
data(pbc, package = "randomForestSRC")
pbc.raw <- na.omit(pbc)
## first we run rhf directly on the ordinary survival data
o.rhf <- rhf(Surv(days, status) ~ ., pbc.raw)
h.rhf <- o.rhf$hazard.oob
time <- o.rhf$time.interest
H.rhf <- o.rhf$chf.oob
S.rhf <- exp(-H.rhf)
## next we run rsf, using the same cases and time.interest grid
o.rsf <- randomForestSRC::rfsrc(
Surv(days, status) ~ ., pbc.raw, ntime = time, use.uno = FALSE
)
H.rsf <- o.rsf$chf.oob
S.rsf <- o.rsf$survival.oob
## compare C-stat of RHF to RSF
cat("RHF C-index:",
randomForestSRC::get.cindex(o.rhf$yvar[,2],o.rhf$yvar[,3],rowSums(H.rhf, na.rm = TRUE)),"\n")
cat("RSF C-index:",
randomForestSRC::get.cindex(o.rsf$yvar[,1],o.rsf$yvar[,2],rowSums(H.rsf, na.rm = TRUE)),"\n")
cat("RSF C-index (from rsf object):",tail(o.rsf$err.rate, 1),"\n")
## graphical parameters
oldpar <- par(mfrow=c(2,2))
## plot survival results
matplot(time,t(S.rhf),type="l",col=gray(.6),ylab="RHF survival")
matplot(time,t(S.rsf),type="l",col=gray(.6),ylab="RSF survival")
boxplot(t(S.rhf-S.rsf),ylab="survival diff.",xlab="cases",xaxt="n",outline=FALSE,range=1e-10)
abline(h=0,lwd=3,col=2)
## plot hazard direct approach - see plot.rhf() for user friendly function
finite <- is.finite(h.rhf)
matplot(time,t(h.rhf),ylab="RHF hazard",
pch=16,type="n",
ylim=c(0,quantile(h.rhf[finite],.99)))
nO <- lapply(which(rowSums(finite) >= 2L), function(i) {
lines(supsmu(time[finite[i,]],h.rhf[i,finite[i,]],bass=1),
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")
f <- Surv(ttodead, died) ~ .
## build the forest directly from ordinary right-censored data
o1 <- rhf(f, peakVO2, treesize=10)
o2 <- rhf(f, peakVO2, treesize=30)
## old-school C-stat
cstat.f <- function(o) {randomForestSRC::get.cindex(o$yvar[, 2],
o$yvar[, 3],rowSums(o$chf.oob, na.rm =TRUE))}
cat("C-stat treesize 10:", cstat.f(o1), "\n")
cat("C-stat treesize 30:", cstat.f(o2), "\n")
## 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)
# }