View source: R/partial.rfsrc.R
| partial.rfsrc | R Documentation |
Compute predictions for partial dependence analysis using a fitted forest.
This advanced interface is the workhorse called by
plot.variable(..., partial = TRUE), which provides automatic grid
selection and plotting. Use partial() directly for custom predictor
grids, additional fixed predictors, several survival evaluation times,
or observation-level predictions.
partial() returns predictions for each training observation.
get.partial.plot.data() extracts their means or individual values
for custom plots. Supported outcomes include regression, classification,
multivariate and mixed outcomes, right-censored survival, and competing
risks. See plot.variable.rfsrc for the simpler plotting
interface.
partial.rfsrc(object, oob = TRUE,
partial.type = NULL, partial.xvar = NULL, partial.values = NULL,
partial.xvar2 = NULL, partial.values2 = NULL,
partial.time = NULL, get.tree = NULL, seed = NULL, do.trace = FALSE, ...)
get.partial.plot.data(o, target, m.target = NULL, granule = FALSE)
object |
An object of class |
oob |
Logical. If |
partial.type |
Prediction summary for survival outcomes. For
right-censored survival, choose |
partial.xvar |
Name of the predictor to vary. |
partial.values |
Nonempty vector of finite numeric values for
|
partial.xvar2 |
Optional character vector naming additional
predictors to hold fixed. Names must be distinct and different from
|
partial.values2 |
Numeric vector giving one finite fixed value per
predictor in |
partial.time |
Nonempty vector of finite evaluation times for
|
get.tree |
Optional vector of tree indices. By default all trees
are considered, with OOB eligibility applied when |
seed |
Negative integer specifying the random-number seed. |
do.trace |
Number of seconds between progress updates; |
... |
Additional internal prediction options. Supply graphical arguments to the subsequent plotting call. |
o |
An object returned by |
target |
For classification, a class label or its numeric position
in the selected response's levels. For competing risks, an event code
from |
m.target |
Response name to extract from a multivariate or mixed
|
granule |
Logical used by |
At each value in partial.values, partial() sets
partial.xvar to that value for every training observation and
predicts with the fitted forest. The other predictors stay at each
observation's training values, except for any additional predictors
fixed by partial.xvar2. No trees are refitted.
For a predictor X_j, with no additional predictors fixed, write
the prediction for observation i at value z as
P_i(z) = \widehat f(z, X_{i,-j}).
Averaging these predictions gives the partial dependence at z:
the model's mean prediction when X_j is set to z, averaged
over the training values of the other predictors.
partial() returns the individual P_i(z) values;
get.partial.plot.data() computes their means. With
oob = TRUE, each P_i(z) uses the trees for which
observation i was OOB. This changes the contributing trees,
while the training observations still provide the values over which
the mean is calculated.
Choose a grid for the range of interest, such as equally spaced values or selected training quantiles. Each grid point requires a prediction for every observation, so fewer points reduce computation and output size.
For factors, use the integer codes corresponding to the training
levels. For example, levels c("low", "middle", "high") have
codes 1, 2, and 3. Use match(labels, levels(training.factor))
to obtain codes for selected labels. The extraction helper returns
numeric codes, which can be labeled in the plotting call.
Use partial.xvar2 and partial.values2 to hold additional
predictors at specified values for every observation. For example,
vary wind over a grid while fixing temperature at one value. The
resulting curve averages predictions over the remaining predictors.
To compare curves at several temperatures, repeat the call with a different fixed temperature each time. This also provides the values for a two-dimensional display: one predictor varies along the grid, and the other changes between calls.
get.partial.plot.data() extracts predicted responses for
regression and the selected class probability for classification.
For multivariate and mixed outcomes, partial() computes all
responses in one call. Use m.target in the extraction helper
to choose a response, and target to choose a class when that
response is a factor. Several responses or classes can be extracted
from the same partial object without recalculating predictions.
For competing risks, target selects an event code from
o$event.info$event.type. This differs from plot.variable(),
where target selects the event's position in that vector.
For right-censored survival, "mort" returns mortality,
"surv" returns survival probabilities, and "chf" returns
cumulative hazards. In partial() and its extraction helper,
"rel.freq" returns the same unnormalized mortality as
"mort".
For competing risks, "years.lost" returns expected years lost,
"cif" returns cumulative incidence, and "chf" returns
cause-specific cumulative hazard. The extraction helper selects the
event of interest.
Supply partial.time for "surv", "cif", or
"chf". Predictions use the forest's time.interest grid.
Each request uses the last grid time less than or equal to it;
a request after the final grid time uses that final time. Using
values from time.interest gives exact matches. The returned
partial.time contains the requested times, which can differ
from the grid times used for prediction.
Mortality and years lost are scalar summaries and do not use
partial.time.
Pass the partial object to get.partial.plot.data() to obtain
predictor values x and corresponding predictions yhat.
The granule argument controls the summary:
granule = FALSEReturn the mean at each grid value.
Only finite predictions contribute to the mean; if none are
finite, the mean is NA.
granule = TRUEReturn the individual predictions,
grouped by grid value and then in training-observation order.
Each grid value is repeated once per observation. Nonfinite
predictions become NA.
For time-dependent survival summaries, yhat has one column
per requested time. The individual values are on the original
prediction scale. Unlike the boxplot values from
plot.variable(..., partial = TRUE), they are not shrunk toward
their mean.
For automatic plots and numerical summaries, use
plot.variable(..., partial = TRUE). Its partial.summary
component contains means, standard deviations, row-averaging standard
errors, and finite prediction counts. See
plot.variable.rfsrc for the interpretation of its
standard errors and display bands.
partial() returns an object of class (rfsrc, partial, family).
It includes call, family, partial.values (the grid
used), partial.time (the requested times), yvar.names, and
event.info, together with the prediction components below.
Let n be the number of training observations, K the number
of predictor grid values, T the number of requested times, and
J the number of events.
regrOutput |
Named list with one |
classOutput |
Named list with one
|
survOutput |
For right-censored mortality ( |
get.partial.plot.data() returns a list with:
x |
Predictor grid values. With |
yhat |
Mean predictions when |
partial.time |
The evaluation times supplied to |
Hemant Ishwaran and Udaya B. Kogalur
Ishwaran H., 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.
plot.variable.rfsrc
## ------------------------------------------------------------
##
## Regression: explicit grid and raw observation-level output.
##
## ------------------------------------------------------------
set.seed(19)
dta <- na.omit(airquality)
fit <- rfsrc(Ozone ~ ., data = dta, ntree = 100)
wind <- unique(as.numeric(quantile(fit$xvar$Wind, seq(0, 1, length.out = 11))))
p <- partial(fit, partial.xvar = "Wind", partial.values = wind)
pd <- get.partial.plot.data(p)
print(head(data.frame(Wind = pd$x, partial = pd$yhat)))
plot(pd$x, pd$yhat, type = "b", xlab = "Wind", ylab = "Partial prediction")
raw <- get.partial.plot.data(p, granule = TRUE)
boxplot(raw$yhat ~ raw$x, xlab = "Wind", ylab = "Individual partial predictions")
## ------------------------------------------------------------
##
## Vary wind while fixing temperature at three values.
##
## ------------------------------------------------------------
temperature <- as.numeric(quantile(fit$xvar$Temp, c(.25, .5, .75)))
slices <- lapply(temperature, function(tt) {
pp <- partial(fit, partial.xvar = "Wind", partial.values = wind,
partial.xvar2 = "Temp", partial.values2 = tt)
get.partial.plot.data(pp)$yhat
})
matplot(wind, do.call(cbind, slices), type = "l", lty = 1,
xlab = "Wind", ylab = "Partial prediction")
legend("topright", legend = paste("Temp =", temperature),
col = seq_along(temperature), lty = 1)
## ------------------------------------------------------------
##
## regression: partial effects for three variables simultaneously
## (can be slow, so modify accordingly)
##
## ------------------------------------------------------------
n <- 200
x <- matrix(rnorm(n * 3), ncol = 3)
y <- x[, 1] + x[, 1] * x[, 2] + x[, 1] * x[, 2] * x[, 3]
o <- rfsrc(y ~ ., data = data.frame(y = y, x))
## define target x values
x1 <- seq(-3, 3, length = 40)
x2 <- x3 <- seq(-3, 3, length = 10)
## extract second order partial effects
pdta <- do.call(rbind,
lapply(x3, function(x3v) {
cat("outer loop x3 = ", x3v, "\n")
do.call(rbind,lapply(x2, function(x2v) {
o <- partial(o,
partial.xvar = "X1",
partial.values = x1,
partial.xvar2 = c("X2", "X3"),
partial.values2 = c(x2v, x3v))
cbind(x1, x2v, x3v, get.partial.plot.data(o)$yhat)
}))
}))
pdta <- data.frame(pdta)
colnames(pdta) <- c("x1", "x2", "x3", "effectSize")
## coplot of partial effects
coplot(effectSize ~ x1|x2*x3, pdta, pch = 16, overlap = 0)
## ------------------------------------------------------------
##
## Classification: extract a named class from a univariate response.
##
## ------------------------------------------------------------
iris.fit <- rfsrc(Species ~ ., data = iris, ntree = 100)
sp <- partial(iris.fit, partial.xvar = "Sepal.Length",
partial.values = seq(4.5, 7.5, length.out = 11))
sp.data <- get.partial.plot.data(sp, target = "versicolor")
plot(sp.data$x, sp.data$yhat, type = "b", ylim = c(0, 1),
xlab = "Sepal length", ylab = "Partial probability of versicolor")
## ------------------------------------------------------------
##
## Multivariate regression: one partial call, two extracted responses.
##
## ------------------------------------------------------------
mv <- rfsrc(cbind(Ozone, Temp) ~ ., data = dta, ntree = 100)
mp <- partial(mv, partial.xvar = "Wind", partial.values = wind)
print(names(mp$regrOutput))
mp.ozone <- get.partial.plot.data(mp, m.target = "Ozone")
mp.temp <- get.partial.plot.data(mp, m.target = "Temp")
print(head(data.frame(Wind = mp.temp$x, Temp = mp.temp$yhat)))
## ------------------------------------------------------------
##
## Mixed outcomes: extract a named class from the selected response.
##
## ------------------------------------------------------------
car.mix <- mtcars
car.mix$am <- factor(car.mix$am)
mix <- rfsrc(Multivar(mpg, am) ~ ., data = car.mix, ntree = 100)
mix.p <- partial(mix, partial.xvar = "wt", partial.values = c(2, 3, 4))
mix.am <- get.partial.plot.data(mix.p, m.target = "am", target = "1")
print(data.frame(wt = mix.am$x, probability = mix.am$yhat))
## ------------------------------------------------------------
##
## Factors: keep the training level codes when choosing labels.
##
## ------------------------------------------------------------
car.data <- mtcars
car.data$cyl <- factor(car.data$cyl)
car.fit <- rfsrc(mpg ~ ., data = car.data, ntree = 100)
lev <- levels(car.fit$xvar$cyl)
code <- match(lev, levels(car.fit$xvar$cyl))
cp <- partial(car.fit, partial.xvar = "cyl", partial.values = code)
cd <- get.partial.plot.data(cp, granule = TRUE)
boxplot(cd$yhat ~ factor(cd$x, levels = code, labels = lev),
xlab = "Cylinders", ylab = "Individual partial predictions")
## ------------------------------------------------------------
##
## Right-censored survival: several stored evaluation times.
##
## ------------------------------------------------------------
data(veteran, package = "randomForestSRC")
vfit <- rfsrc(Surv(time, status) ~ ., data = veteran, ntree = 100)
grid <- vfit$time.interest
when <- grid[unique(round(seq(1, length(grid), length.out = 3)))]
karno <- unique(as.numeric(quantile(vfit$xvar$karno)))
vp <- partial(vfit, partial.type = "surv", partial.xvar = "karno",
partial.values = karno, partial.time = when)
vd <- get.partial.plot.data(vp)
print(vd$partial.time)
matplot(vd$x, vd$yhat, type = "l", lty = 1,
xlab = "Karnofsky score", ylab = "Partial survival probability")
legend("bottomright", legend = paste("Time =", vd$partial.time),
col = seq_along(vd$partial.time), lty = 1)
## ------------------------------------------------------------
##
## Survival curves at selected quantiles of Karnofsky score.
## Reuse vfit from the preceding right-censored survival example.
##
## ------------------------------------------------------------
karno.probs <- c(.10, .50, .90)
karno.values <- quantile(vfit$xvar$karno, probs = karno.probs,
na.rm = TRUE)
print(karno.values)
## Set Karnofsky score to each requested value for every reference row.
## Keep the other predictors at their observed values and use OOB trees.
surv.partial <- partial(vfit,
oob = TRUE,
partial.type = "surv",
partial.xvar = "karno",
partial.values = as.numeric(karno.values),
partial.time = vfit$time.interest)
surv.data <- get.partial.plot.data(surv.partial)
## yhat has one row per distinct Karnofsky value and one column per time.
## Transpose it so that each plotted column is a survival curve over time.
## Each curve averages over the same reference rows.
curve.id <- seq_along(surv.data$x)
matplot(surv.data$partial.time, t(surv.data$yhat),
type = "s", col = curve.id, lty = curve.id, lwd = 2,
ylim = c(0, 1), xlab = "Time", ylab = "Survival probability",
main = "Partial survival curves by Karnofsky score")
## Use the returned grid for labels: partial() removes tied quantile values.
legend("topright",
legend = paste("Karnofsky =", format(surv.data$x, trim = TRUE)),
col = curve.id, lty = curve.id, lwd = 2, bty = "n")
## ------------------------------------------------------------
##
## Competing risks: select an event code from the saved event map.
##
## ------------------------------------------------------------
data(follic, package = "randomForestSRC")
crfit <- rfsrc(Surv(time, status) ~ ., data = follic, ntree = 100, nsplit = 3)
age <- unique(as.numeric(quantile(crfit$xvar$age)))
crp <- partial(crfit, partial.type = "years.lost",
partial.xvar = "age", partial.values = age)
event <- crp$event.info$event.type[1]
crd <- get.partial.plot.data(crp, target = event)
plot(crd$x, crd$yhat, type = "b", xlab = "Age", ylab = "Partial years lost")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.