28 Which variables contribute to prediction?

28.1 When to use it

Once a forest predicts adequately, a co-author usually asks which variables carry that prediction. Permutation variable importance (VIMP) answers by breaking one predictor’s information in the OOB data and measuring how much prediction error changes.

That is a model question. VIMP is not an effect size, does not give the direction of association, and does not show that a variable causes the outcome. Correlated predictors can share predictive information, so permuting one may redistribute or understate the importance we would see if it stood alone.

28.2 Fit a classification forest

We again use iris so the returned object exposes both global and class-specific VIMP. importance = TRUE asks rfsrc() to retain permutation importance. The bounded 300-tree fit is for teaching; a study analysis should verify both the error trajectory and the ranking’s stability across an adequate forest.

set.seed(20260828)
rf <- rfsrc(
  Species ~ .,
  data = iris,
  ntree = 300,
  importance = TRUE
)
rf

                         Sample size: 150
           Frequency of class labels: setosa=50, versicolor=50, virginica=50
                     Number of trees: 300
           Forest terminal node size: 1
       Average no. of terminal nodes: 9.5433
No. of variables tried at each split: 2
              Total no. of variables: 4
       Resampling used to grow trees: swor
    Resample size used to grow trees: 95
                            Analysis: RF-C
                              Family: class
                      Splitting rule: gini *random*
       Number of random split points: 10
                   (OOB) Brier score: 0.02448324
        (OOB) Normalized Brier score: 0.11017457
                           (OOB) AUC: 0.9937
                      (OOB) Log-loss: 0.12563704
   (OOB) Requested performance error: 0.05333333, 0.02, 0.06, 0.08

Confusion matrix:

            predicted
  observed   setosa versicolor virginica class.error
  setosa         49          1         0        0.02
  versicolor      0         47         3        0.06
  virginica       0          4        46        0.08

      (OOB) Misclassification rate: 0.05333333

Random-classifier baselines (uniform):
   Brier: 0.22222222   Normalized Brier: 1   Log-loss: 1.09861229

28.3 Build and inspect it

vimp_dta <- gg_vimp(rf)
vimp_dta

<gg_vimp>  from randomForestSRC  |  family: class  |  ntree: 300  |  n: 150  |  variables: 16

head(as.data.frame(vimp_dta), 8)

          vars        set      vimp positive
1  Petal.Width versicolor 0.7846667     TRUE
2 Petal.Length versicolor 0.6260000     TRUE
3  Petal.Width        all 0.6187549     TRUE
4  Petal.Width     setosa 0.5880000     TRUE
5 Petal.Length        all 0.5177924     TRUE
6 Petal.Length     setosa 0.4853333     TRUE
7  Petal.Width  virginica 0.4626667     TRUE
8 Petal.Length  virginica 0.4246667     TRUE

with(as.data.frame(vimp_dta),
     aggregate(vimp, list(outcome = set), max))

     outcome         x
1        all 0.6187549
2     setosa 0.5880000
3 versicolor 0.7846667
4  virginica 0.4626667

The vars column holds raw predictor identities. set = "all" is the global permutation score; the named class sets show the one-vs-rest contribution for each outcome. Positive VIMP means permutation worsened prediction. Negative VIMP means permutation improved it in this fit, often a sign of noise or an unstable small contribution. The sign is not the direction of a clinical association.

28.4 Use display labels at plot time

Human-readable labels belong in the display, not in the model object. The current plot method uses labels. We confirm that applying those labels leaves the returned variables and their order unchanged.

feature_labels <- c(
  Sepal.Length = "Sepal length",
  Sepal.Width = "Sepal width",
  Petal.Length = "Petal length",
  Petal.Width = "Petal width"
)

raw_ids <- as.character(vimp_dta$vars)
vimp_plot <- plot(vimp_dta, labels = feature_labels)

data.frame(
  raw_names_unchanged = identical(raw_ids, as.character(vimp_dta$vars)),
  rows = nrow(vimp_dta),
  outcome_sets = length(unique(vimp_dta$set))
)

  raw_names_unchanged rows outcome_sets
1                TRUE   16            4

vimp_plot +
  geom_hline(yintercept = 0, color = "grey45", linewidth = 0.35) +
  theme_hv_manuscript() +
  scale_y_continuous(
    breaks = function(x) pretty(x, n = 3),
    expand = expansion(mult = c(0.05, 0.05))
  ) +
  labs(x = "Predictor", y = "Permutation VIMP", fill = "VIMP > 0") +
  theme(axis.text.x = element_text(size = 7))

Figure 28.1: Permutation VIMP overall and by iris class; display labels do not alter the stored predictor names

28.5 Read it

Start with the all facet for the global ranking, then inspect whether the class-specific facets tell the same story. Petal width and petal length carry most of the prediction in this fit. A predictor can rank differently by class because the forest may need it to distinguish one outcome more than another. The grey zero line separates positive from negative VIMP, and the axis range is set from the returned scores so neither negative values nor larger positive values are clipped.

The object contains point estimates, not confidence intervals. Small gaps near the bottom should not be treated as settled ranks. If the order matters to a manuscript conclusion, refit across seeds or use a resampling analysis and report that uncertainty separately.

28.6 Adapt it

For a survival or regression forest there is usually one outcome set. For a classification forest, inspect unique(vimp_dta$set) before deciding whether to report the global ranking or a class-specific one. gg_vimp(rf, nvar = 10) can reduce a long display, but keep the full ranking in a supplement when the cutoff affects interpretation.

28.7 Pitfalls

Reading sign as clinical direction. Positive VIMP means worse prediction after permutation, not that larger predictor values increase risk.
Calling VIMP an effect size. Scores depend on this forest, this outcome, this predictor set, and the prediction-error scale.
Ignoring correlated predictors. Shared information can be split or redistributed across variables.
Treating labels as data changes. Use labels in plot(). Keep the raw predictor names intact for downstream code and audit trails.
Claiming causality. VIMP describes predictive contribution. It does not identify an intervention effect.