PD and ICE curves with ORSF
Partial dependence (PD)
Partial dependence (PD) shows the expected prediction from a model as
a function of a single predictor or multiple predictors. The expectation
is marginalized over the values of all other predictors, giving
something like a multivariable adjusted estimate of the model’s
prediction.
Begin by fitting an ORSF ensemble. Set a prediction horizon of 5
years when we fit the ensemble so that any aorsf function
that we pass this ensemble to will assume we want to compute predictions
at 5 years.
library(aorsf)
pred_horizon <- 365.25 * 5
set.seed(329730)
index_train <- sample(nrow(pbc_orsf), 150)
pbc_orsf_train <- pbc_orsf[index_train, ]
pbc_orsf_test <- pbc_orsf[-index_train, ]
fit <- orsf(data = pbc_orsf_train,
formula = Surv(time, status) ~ . - id,
oobag_pred_horizon = pred_horizon)
fit
#> ---------- Oblique random survival forest
#>
#> Linear combinations: Accelerated
#> N observations: 150
#> N events: 52
#> N trees: 500
#> N predictors total: 17
#> N predictors per node: 5
#> Average leaves per tree: 12
#> Min observations in leaf: 5
#> Min events in leaf: 1
#> OOB stat value: 0.83
#> OOB stat type: Harrell's C-statistic
#> Variable importance: anova
#>
#> -----------------------------------------
Three ways to compute PD
You can compute PD three ways with aorsf:
using in-bag predictions for the training data
pd_inb <- orsf_pd_inb(fit, pred_spec = list(bili = 1:5))
pd_inb
#> pred_horizon bili mean lwr medn upr
#> 1: 1826.25 1 0.2054232 0.01599366 0.0929227 0.8077278
#> 2: 1826.25 2 0.2369077 0.02549869 0.1268457 0.8227315
#> 3: 1826.25 3 0.2808514 0.05027265 0.1720280 0.8457834
#> 4: 1826.25 4 0.3428065 0.09758988 0.2545869 0.8575243
#> 5: 1826.25 5 0.3992909 0.16392752 0.3232681 0.8634269
using out-of-bag predictions for the training data
pd_oob <- orsf_pd_oob(fit, pred_spec = list(bili = 1:5))
pd_oob
#> pred_horizon bili mean lwr medn upr
#> 1: 1826.25 1 0.2068300 0.01479443 0.08824123 0.8053317
#> 2: 1826.25 2 0.2377046 0.02469718 0.12623031 0.8258154
#> 3: 1826.25 3 0.2810546 0.04080813 0.18721220 0.8484846
#> 4: 1826.25 4 0.3417839 0.09076851 0.24968438 0.8611884
#> 5: 1826.25 5 0.3979925 0.16098228 0.32147532 0.8554402
using predictions for a new set of data
pd_test <- orsf_pd_new(fit,
new_data = pbc_orsf_test,
pred_spec = list(bili = 1:5))
pd_test
#> pred_horizon bili mean lwr medn upr
#> 1: 1826.25 1 0.2510900 0.01631318 0.1872414 0.8162621
#> 2: 1826.25 2 0.2807327 0.02903956 0.2269297 0.8332956
#> 3: 1826.25 3 0.3247386 0.05860235 0.2841853 0.8481825
#> 4: 1826.25 4 0.3850799 0.10741224 0.3405760 0.8588955
#> 5: 1826.25 5 0.4394952 0.17572657 0.4050864 0.8657886
in-bag PD indicates relationships that the model has learned
during training. This is helpful if your goal is to interpret the
model.
out-of-bag PD indicates relationships that the model has learned
during training but using the out-of-bag data simulates application of
the model to new data. if you want to test your model’s reliability or
fairness in new data but you don’t have access to a large testing
set.
new data PD shows how the model predicts outcomes for
observations it has not seen. This is helpful if you want to test your
model’s reliability or fairness.
Let’s re-fit our ORSF to all available data before proceeding to the
next sections.
set.seed(329730)
fit <- orsf(pbc_orsf,
Surv(time, status) ~ . -id,
oobag_pred_horizon = pred_horizon)
One variable, one horizon
Computing PD for a single variable is straightforward:
pd_sex <- orsf_pd_oob(fit, pred_spec = list(sex = c("m", "f")))
pd_sex
#> pred_horizon sex mean lwr medn upr
#> 1: 1826.25 m 0.3552858 0.03995954 0.2427792 0.9474016
#> 2: 1826.25 f 0.2951752 0.01132241 0.1435909 0.9625991
The output shows that the expected predicted mortality risk for men
is substantially higher than women at 5 years after baseline.
One variable, moving horizon
What if the effect of a predictor varies over time? PD can show
this.
pd_sex_tv <- orsf_pd_oob(fit, pred_spec = list(sex = c("m", "f")),
pred_horizon = seq(365, 365*5))
ggplot(pd_sex_tv, aes(x = pred_horizon, y = mean, color = sex)) +
geom_line() +
labs(x = 'Time since baseline',
y = 'Expected risk')
From inspection, we can see that males have higher risk than females
and the difference in that risk grows over time. This can also be seen
by viewing the ratio of expected risk over time:
library(data.table)
ratio_tv <- pd_sex_tv[
, .(ratio = mean[sex == 'm'] / mean[sex == 'f']), by = pred_horizon
]
ggplot(ratio_tv, aes(x = pred_horizon, y = ratio)) +
geom_line(color = 'grey') +
geom_smooth(color = 'black', se = FALSE) +
labs(x = 'time since baseline',
y = 'ratio in expected risk for males versus females')
#> `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
Multiple variables, marginally
If you want to compute PD marginally for multiple variables, just
list the variable values in pred_spec and specify
expand_grid = FALSE.
pd_two_vars <-
orsf_pd_oob(fit,
pred_spec = list(sex = c("m", "f"), bili = 1:5),
expand_grid = FALSE)
pd_two_vars
#> pred_horizon variable value level mean lwr medn upr
#> 1: 1826.25 sex NA m 0.3552858 0.03995954 0.2427792 0.9474016
#> 2: 1826.25 sex NA f 0.2951752 0.01132241 0.1435909 0.9625991
#> 3: 1826.25 bili 1 <NA> 0.2342067 0.01251527 0.1212992 0.8733544
#> 4: 1826.25 bili 2 <NA> 0.2895945 0.04965949 0.1920790 0.9097945
#> 5: 1826.25 bili 3 <NA> 0.3475991 0.07248908 0.2557485 0.9208106
#> 6: 1826.25 bili 4 <NA> 0.3959820 0.09116718 0.3278895 0.9334875
#> 7: 1826.25 bili 5 <NA> 0.4415943 0.12358330 0.3921014 0.9375022
Now would it be tedious if you wanted to do this for all the
variables? You bet. That’s why we made a function for that. As a bonus,
the printed output is sorted from most to least important variables.
pd_smry <- orsf_summarize_uni(fit)
pd_smry
#>
#> -- bili (VI Rank: 1) -------------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 0.80 0.2294737 0.1157904 0.04125087 0.3670285
#> 1.4 0.2508836 0.1351933 0.05464883 0.3875328
#> 3.5 0.3743323 0.2932788 0.16292560 0.5585283
#>
#> -- age (VI Rank: 2) --------------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 41 0.2731750 0.1319029 0.03878754 0.4605820
#> 50 0.3009128 0.1548903 0.04441111 0.5395029
#> 57 0.3324560 0.2082692 0.06778463 0.5768839
#>
#> -- copper (VI Rank: 3) -----------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 43 0.2659530 0.1288386 0.04283502 0.4644505
#> 74 0.2819925 0.1467904 0.05306841 0.4942437
#> 130 0.3326852 0.2110308 0.09796624 0.5333395
#>
#> -- protime (VI Rank: 4) ----------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 10 0.2819480 0.1472897 0.04428083 0.5111052
#> 11 0.2947989 0.1510190 0.04779299 0.5311810
#> 11 0.3175249 0.1843402 0.06216519 0.5486052
#>
#> -- ascites (VI Rank: 5) ----------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 0 0.2960603 0.1435909 0.04608949 0.5373328
#> 1 0.4681807 0.3841361 0.26174441 0.6476854
#>
#> -- stage (VI Rank: 6) ------------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 1 0.658076 0.6070399 0.5189756 0.7865766
#> 2 0.658076 0.6070399 0.5189756 0.7865766
#> 3 0.658076 0.6070399 0.5189756 0.7865766
#> 4 0.658076 0.6070399 0.5189756 0.7865766
#>
#> -- chol (VI Rank: 7) -------------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 249 0.2883441 0.1345402 0.03932400 0.5131917
#> 310 0.2958828 0.1468409 0.04305168 0.5290724
#> 402 0.3164329 0.1858660 0.07206203 0.5453219
#>
#> -- sex (VI Rank: 8) --------------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> m 0.3552858 0.2427792 0.11136359 0.5912710
#> f 0.2951752 0.1435909 0.04511121 0.5262298
#>
#> -- albumin (VI Rank: 9) ----------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 3.3 0.3174609 0.1792141 0.05618604 0.5850415
#> 3.5 0.2943950 0.1488761 0.04425814 0.5428160
#> 3.8 0.2789835 0.1439382 0.04309754 0.4997154
#>
#> -- spiders (VI Rank: 10) ---------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 0 0.2932967 0.1409483 0.04148682 0.5359432
#> 1 0.3333689 0.1950764 0.07789822 0.5655366
#>
#> -- edema (VI Rank: 11) -----------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 0 0.2916386 0.1412817 0.04456894 0.5346685
#> 0.5 0.3545577 0.2253930 0.09172178 0.6284023
#> 1 0.4477276 0.3450991 0.24036516 0.6662038
#>
#> -- ast (VI Rank: 12) -------------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 82 0.2859670 0.1445626 0.04299327 0.5266939
#> 117 0.2998852 0.1609152 0.04946382 0.5527030
#> 153 0.3193772 0.1721107 0.06368975 0.5863355
#>
#> -- platelet (VI Rank: 13) --------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 200 0.3086503 0.1609407 0.04633406 0.5939281
#> 257 0.3025122 0.1546509 0.04495160 0.5723953
#> 318 0.2995643 0.1456786 0.04767242 0.5555624
#>
#> -- hepato (VI Rank: 14) ----------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 0 0.2886564 0.1449674 0.04246348 0.5374039
#> 1 0.3186173 0.1843793 0.06072296 0.5538878
#>
#> -- trig (VI Rank: 15) ------------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 85 0.2957353 0.1428190 0.03972731 0.5358985
#> 108 0.3012346 0.1540859 0.04217322 0.5560627
#> 152 0.3134312 0.1796233 0.04994688 0.5564560
#>
#> -- trt (VI Rank: 16) -------------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> d_penicill_main 0.3085388 0.1756992 0.04966480 0.5487896
#> placebo 0.3009444 0.1552409 0.04556513 0.5611608
#>
#> -- alk.phos (VI Rank: 17) --------------------------------
#>
#> |---------------- risk ----------------|
#> Value Mean Median 25th % 75th %
#> 921 0.3009482 0.1433016 0.04407730 0.5562279
#> 1278 0.3032376 0.1476446 0.04521654 0.5564515
#> 2072 0.3084310 0.1540810 0.04858680 0.5649115
#>
#> Predicted risk at time t = 1826.25 for top 17 predictors
It’s easy enough to turn this ‘summary’ object into a
data.table for downstream plotting and tables.
head(as.data.table(pd_smry))
#> variable importance Value Mean Median 25th % 75th %
#> 1: bili 0.08991457 0.80 0.2294737 0.1157904 0.04125087 0.3670285
#> 2: bili 0.08991457 1.4 0.2508836 0.1351933 0.05464883 0.3875328
#> 3: bili 0.08991457 3.5 0.3743323 0.2932788 0.16292560 0.5585283
#> 4: age 0.02156699 41 0.2731750 0.1319029 0.03878754 0.4605820
#> 5: age 0.02156699 50 0.3009128 0.1548903 0.04441111 0.5395029
#> 6: age 0.02156699 57 0.3324560 0.2082692 0.06778463 0.5768839
#> pred_horizon level
#> 1: 1826.25 <NA>
#> 2: 1826.25 <NA>
#> 3: 1826.25 <NA>
#> 4: 1826.25 <NA>
#> 5: 1826.25 <NA>
#> 6: 1826.25 <NA>
Multiple variables, jointly
PD can show the expected value of a model’s predictions as a function
of a specific predictor, or as a function of multiple predictors. For
instance, we can estimate predicted risk as a joint function of
bili, edema, and trt:
pred_spec = list(bili = seq(1, 5, length.out = 20),
edema = levels(pbc_orsf_train$edema),
trt = levels(pbc_orsf$trt))
pd_bili_edema <- orsf_pd_oob(fit, pred_spec)
library(ggplot2)
ggplot(pd_bili_edema, aes(x = bili, y = medn, col = trt, linetype = edema)) +
geom_line() +
labs(y = 'Expected predicted risk')
From inspection,
the model’s predictions indicate slightly lower risk for the
placebo group, and these do not seem to change much at different values
of bili or edema.
There is a clear increase in predicted risk with higher levels of
edema and with higher levels of bili
the slope of predicted risk as a function of bili
appears highest among patients with edema of 0.5. Is the
effect of bili modified by edema being 0.5? A
quick sanity check with coxph suggests there is.
library(survival)
pbc_orsf$edema_05 <- ifelse(pbc_orsf$edema == '0.5', 'yes', 'no')
fit_cph <- coxph(Surv(time,status) ~ edema_05 * bili,
data = pbc_orsf)
anova(fit_cph)
#> Analysis of Deviance Table
#> Cox model: response is Surv(time, status)
#> Terms added sequentially (first to last)
#>
#> loglik Chisq Df Pr(>|Chi|)
#> NULL -550.19
#> edema_05 -546.83 6.7248 1 0.009508 **
#> bili -513.59 66.4689 1 3.555e-16 ***
#> edema_05:bili -510.54 6.1112 1 0.013433 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Individual conditional expectations (ICE)
Unlike partial dependence, which shows the expected prediction as a
function of one or multiple predictors, individual conditional
expectations (ICE) show the prediction for an individual observation as
a function of a predictor.
Just like PD, we can compute ICE using in-bag, out-of-bag, or testing
data, and the same principles apply. We’ll use out-of-bag estimates
here.
Visualizing ICE curves
Inspecting the ICE curves for each observation can help identify
whether there is heterogeneity in a model’s predictions. I.e., does the
effect of the variable follow the same pattern for all the data, or are
there groups where the variable impacts risk differently?
I am going to turn off boundary checking in orsf_ice_oob
by setting boundary_checks = FALSE, and this will allow me
to generate ICE curves that go beyond the 90th percentile of
bili.
pred_spec <- list(bili = seq(1, 10, length.out = 25))
ice_oob <- orsf_ice_oob(fit, pred_spec, boundary_checks = FALSE)
ice_oob
#> pred_horizon id_variable id_row bili pred
#> 1: 1826.25 1 1 1 0.9035317
#> 2: 1826.25 1 2 1 0.1012207
#> 3: 1826.25 1 3 1 0.6983324
#> 4: 1826.25 1 4 1 0.3368231
#> 5: 1826.25 1 5 1 0.1054981
#> ---
#> 6896: 1826.25 25 272 10 0.4456932
#> 6897: 1826.25 25 273 10 0.5718185
#> 6898: 1826.25 25 274 10 0.4973769
#> 6899: 1826.25 25 275 10 0.3913366
#> 6900: 1826.25 25 276 10 0.5289952
id_variable is an identifier for the current value
of the variable(s) that are in the data. It is redundant if you only
have one variable, but helpful if there are multiple variables.
id_row is an identifier for the observation in the
original data. It is used to group an observation’s predictions together
in plots.
For plots, it is helpful to scale the ICE data. I subtract the
initial value of predicted risk (i.e., when bili = 1) from
each observation’s conditional expectation values. So,
Now we can visualize the curves.
library(ggplot2)
ggplot(ice_oob, aes(x = bili,
y = pred,
group = id_row)) +
geom_line(alpha = 0.15) +
labs(y = 'Change in predicted risk') +
geom_smooth(se = FALSE, aes(group = 1))
#> `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
From inspection of the figure,
Limitations of PD
Partial dependence has a number of known
limitations and assumptions that users should be aware of (see
Hooker, 2021). In particular, partial dependence is less intuitive when
>2 predictors are examined jointly, and it is assumed that the
feature(s) for which the partial dependence is computed are not
correlated with other features (this is likely not true in many cases).
Accumulated local effect plots can be used (see here)
in the case where feature independence is not a valid assumption.
References
- Giles Hooker, Lucas Mentch, Siyu Zhou. Unrestricted Permutation
forces Extrapolation: Variable Importance Requires at least One More
Model, or There Is No Free Variable Importance. arXiv e-prints
2021 Oct; arXiv-1905. URL: https://doi.org/10.48550/arXiv.1905.03151