--- title: "Introduction to the ale package" author: "Chitu Okoli" date: 2026-10-05 date-format: "MMMM D, YYYY" format: html: toc: true toc-location: left toc-title: "Outline" number-sections: true code-fold: false smooth-scroll: true vignette: > %\VignetteIndexEntry{Introduction to the ale package} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} --- Accumulated Local Effects (ALE) were initially developed as a model-agnostic approach for global explanations of the results of black-box machine learning algorithms (Apley, Daniel W., and Jingyu Zhu. 'Visualizing the effects of predictor variables in black box supervised learning models.' Journal of the Royal Statistical Society Series B: Statistical Methodology 82.4 (2020): 1059-1086 ). ALE has at least two primary advantages over other approaches like partial dependency plots (PDP) and SHapley Additive exPlanations (SHAP): its values are not affected by the presence of interactions among variables in a model and its computation is relatively rapid. This package reimplements the original algorithm from the [`{ALEPlot}` package](https://CRAN.r-project.org/package=ALEPlot) and reimplements the plotting of ALE values. It also extends the original ALE concept to add bootstrap-based confidence intervals and ALE-based statistics that can be used for statistical inference. For more details, see Okoli, Chitu. 2023. “Statistical Inference Using Machine Learning and Classical Techniques Based on Accumulated Local Effects (ALE).” arXiv. . This vignette demonstrates the basic functionality of the `{ale}` package on standard large datasets used for machine learning. A separate vignette is devoted to its use on [small datasets](ale-small-datasets.html "ale package for small datasets"), as is often the case with statistical inference. (How small is small? That's a tough question, but as that vignette explains, most datasets of less than 2000 rows are probably "small" and even many datasets that are more than 2000 rows are nonetheless "small".) Other vignettes introduce [ALE-based statistics for statistical inference](ale-statistics.html) and show how the `{ale}` package handles [various datatypes of input variables](ale-x-datatypes.html). We begin by loading the necessary libraries. ```{r load libraries} library(ale) library(dplyr) ``` ## diamonds dataset For this introduction, we use the `diamonds` dataset, included with the `{ggplot2}` graphics system. We cleaned the original version by [removing duplicates](https://lorentzen.ch/index.php/2021/04/16/a-curious-fact-on-the-diamonds-dataset/ "errors in the diamonds dataset") and invalid entries where the length (x), width (y), or depth (z) is 0. ```{r diamonds_print} # Clean up some invalid entries diamonds <- ggplot2::diamonds |> filter(!(x == 0 | y == 0 | z == 0)) |> # https://lorentzen.ch/index.php/2021/04/16/a-curious-fact-on-the-diamonds-dataset/ distinct( price, carat, cut, color, clarity, .keep_all = TRUE ) |> rename( x_length = x, y_width = y, z_depth = z, depth_pct = depth ) summary(diamonds) ``` Here is the description of the modified dataset. | Variable | Description | |-----------------|-------------------------------------------------------| | price | price in US dollars (\$326--\$18,823) | | carat | weight of the diamond (0.2--5.01) | | cut | quality of the cut (Fair, Good, Very Good, Premium, Ideal) | | color | diamond color, from D (best) to J (worst) | | clarity | a measurement of how clear the diamond is (I1 (worst), SI2, SI1, VS2, VS1, VVS2, VVS1, IF (best)) | | x_length | length in mm (0--10.74) | | y_width | width in mm (0--58.9) | | z_depth | depth in mm (0--31.8) | | depth_pct | total depth percentage = z / mean(x, y) = 2 \* z / (x + y) (43--79) | | table | width of top of diamond relative to widest point (43--95) | ```{r diamonds_str} str(diamonds) ``` ```{r diamonds_price} summary(diamonds$price) ``` Interpretable machine learning (IML) techniques like ALE should be applied on the same dataset that was used to train the model. An explanation is an explanation of a trained model and a trained model is intrinsically linked to the dataset on which it is trained. (When a dataset is too small to feasibly split into training and test sets, then the ale package has tools to appropriately handle such [small datasets](ale-small-datasets.html "ale package for small datasets"). ## Modelling with a random forest ALE is a model-agnostic IML approach, that is, it works with any kind of machine learning model. As such, `{ale}` works with any R model with the only condition that it can predict numeric outcomes (such as raw estimates for regression and probabilities or odds ratios for classification). For this demonstration, we use a random forest from the `{ranger}` package. A random forest combines predictions from multiple decision trees and can capture nonlinear effects and interactions. The examples here will work with any statistical or machine learning algorithm that produces numeric predictions. We train a random forest to predict diamond prices: ```{r train_rf} rf_diamonds <- ranger::ranger( price ~ ., diamonds, # train for IML on full dataset mtry = 7, # number of variables per tree min.node.size = 4, # minimum leaf (terminal node) size sample.fraction = 0.465, # fraction of rows to sample # small forest still accurate for a large dataset; 10× faster than 500 trees num.trees = 50, seed = 0 # keep the random forest consistent ) rf_diamonds ``` ## Creating an `ALE` object The core object in the `{ale}` package is the `{S7}` `ALE` object. Its `effect` property stores the ALE data and, optionally, ALE statistics and bootstrap data for one or more outcome categories. The first argument to the `ALE()` constructor is a model object--any R model object that can generate numeric predictions is acceptable. By default, `ALE()` calculates 1D (or "first-order") ALE for every predictor in the model's training data. Ranger models do not retain their training data, so we supply `data = diamonds` explicitly. The outcome column can be detected from this model; for models where it cannot be detected, also supply `y_col`. ```{r parallel} # For faster processing, set the number of CPU cores available. See help(ALE). options(ale.parallel = 2) ``` For faster demonstrations, this vignette uses precreated ALE objects. For the full experience, you can uncomment the relevant lines in the code below. ```{r serialized_objects_site} # For speed, load pre-created objects from the current wip branch. serialized_objects_site <- "https://github.com/tripartio/ale/raw/wip/download" ``` ```{r ale_simple, eval=FALSE, echo=TRUE} # To run the slow code yourself, uncomment and execute the following call. # ale_rf_diamonds <- ALE( # rf_diamonds, # data = diamonds # ) ale_rf_diamonds <- serialized_objects_site |> file.path("ale_rf_diamonds.0.5.3.rds") |> url() |> readRDS() # Print a statistical summary of the ALE object summary(ale_rf_diamonds) ``` ```{r ale_simple-rds, eval=TRUE, echo=FALSE} ale_rf_diamonds <- retrieve_rds( # For speed, load a pre-created object by default. if (file.exists("../download/ale_rf_diamonds.0.5.3.rds")) { readRDS("../download/ale_rf_diamonds.0.5.3.rds") } else { c(serialized_objects_site, "ale_rf_diamonds.0.5.3.rds") }, { ALE( rf_diamonds, data = diamonds ) } ) # saveRDS(ale_rf_diamonds, "../download/ale_rf_diamonds.0.5.3.rds") # Print a statistical summary of the ALE object summary(ale_rf_diamonds) ``` The `summary()` method reports the model and data represented by the object, followed by its ALE statistics and statistically significant confidence regions when available. It returns the `ALE` object invisibly, so it is intended for inspection rather than data extraction. The default constructor calculates ALE statistics, although p-values and their confidence regions are not calculated unless p-values are supplied or generated during bootstrapping. See the [ALE statistics vignette](ale-statistics.html "ALE-based statistics") for their interpretation. Most core functions in the `{ale}` package support parallel processing. You can enable parallelization globally with `options(ale.parallel = x)`, where "x" is the number of CPU cores you want to use. It can also be enabled for specific function runs with the `parallel` argument. Parallelization is disabled by default with `parallel = 0`. See the documentation for `ALE()` for details. To access the plot for a specific variable, we must first create an `ALEPlots` object by calling the `plot()` method on the `ALE` object which internally generates `ggplot` objects with the full flexibility of {ggplot2}: ```{r create-plots} # Generate an ALEPlots object diamonds_plots <- plot(ale_rf_diamonds) ``` To retrieve a specific variable plot, you can use the `get()` method of the `ALEPlots` object. For example, to access and print the `carat` ALE plot, we can simply refer to `get(diamonds_plots, 'carat')`: ```{r print-carat, fig.width=4, fig.height=3} # Print a plot by entering its reference get(diamonds_plots, 'carat') ``` To display all the ALE plots in an `ALEPlots` object, we can simply call its `print()` or `plot()` methods. Behind the scenes, they use the `{patchwork}` package to arrange multiple plots in a common plot grid using `patchwork::wrap_plots()`, so we can pass arguments from that function. For example, we can specify that we want two plots per row with the `ncol` argument: ```{r print-ale_simple, fig.width=12, fig.height=22} # Print all plots plot(diamonds_plots, ncol = 2) ``` ## Selecting variables and interactions The `x_cols` argument controls the effects that `ALE()` calculates. Its list syntax can request 1D effects, 2D interactions, or both in the same object. For example, the following requests every 1D effect and only the interaction between `carat` and `clarity`: ```{r select-list, eval=FALSE, echo=TRUE} ALE( rf_diamonds, data = diamonds, x_cols = list(d1 = TRUE, d2 = "carat:clarity") ) ``` Formula syntax is convenient when selecting particular effects. Terms on their own request 1D ALE, while terms joined by `:` request 2D ALE: ```{r select-formula, eval=FALSE, echo=TRUE} ALE( rf_diamonds, data = diamonds, x_cols = ~ carat + cut + clarity + carat:clarity + cut:table ) ``` These formats can be combined with other selection options. See `help(resolve_x_cols)` for the complete set of formats, including ways to request all interactions involving a particular variable and to exclude variables. ## Bootstrapped ALE One of the key features of the ALE package is bootstrapping of the ALE results to ensure that the results are reliable, that is, generalizable to data beyond the sample on which the model was trained. As mentioned above, this assumes that IML analysis is carried out on a model whose hyperparameters were determined by cross-validation. When samples are too small for cross-validation, we provide a different approach by bootstrapping the entire model with a `ModelBoot` object, explained in the vignette for [small datasets](ale-small-datasets.html "ale package for small datasets"). Although ALE is faster than most other IML techniques for global explanation such as partial dependence plots (PDP) and SHAP, it still requires some time to run. Bootstrapping multiplies much of that time by the number of bootstrap iterations. We create bootstrapped ALE data and plots using the `boot_it` argument. ALE is a relatively stable IML algorithm (compared to others like PDP), so 100 bootstrap samples should be sufficient for relatively stable results, especially for model development. Final results could be confirmed with 1000 bootstrap samples or more, but there should not be much difference in the results beyond 100 iterations. ```{r ale_boot, eval=FALSE, echo=TRUE, fig.width=12, fig.height=22} # To run the slow code yourself, uncomment and execute the following call. # ale_rf_diamonds_boot <- ALE( # rf_diamonds, # data = diamonds, # # Request all 1D effects and two specific 2D interactions. # x_cols = list(d1 = TRUE, d2 = c("carat:clarity", "cut:table")), # boot_it = 100 # ) ale_rf_diamonds_boot <- serialized_objects_site |> file.path("ale_rf_diamonds_boot.0.5.3.rds") |> url() |> readRDS() # Bootstrapping produces confidence intervals for the 1D effects plot(ale_rf_diamonds_boot) |> subset(list(d1 = TRUE)) |> print(ncol = 2) ``` ```{r ale_boot-rds, eval=TRUE, echo=FALSE, fig.width=12, fig.height=22} ale_rf_diamonds_boot <- retrieve_rds( # For speed, load a pre-created object by default. if (file.exists("../download/ale_rf_diamonds_boot.0.5.3.rds")) { readRDS("../download/ale_rf_diamonds_boot.0.5.3.rds") } else { c(serialized_objects_site, "ale_rf_diamonds_boot.0.5.3.rds") }, { ALE( rf_diamonds, data = diamonds, # Request all 1D effects and two specific 2D interactions. x_cols = list(d1 = TRUE, d2 = c("carat:clarity", "cut:table")), boot_it = 100 ) } ) # saveRDS(ale_rf_diamonds_boot, "../download/ale_rf_diamonds_boot.0.5.3.rds") # Bootstrapping produces confidence intervals for the 1D effects plot(ale_rf_diamonds_boot) |> subset(list(d1 = TRUE)) |> print(ncol = 2) ``` Compare the bootstrapped confidence intervals with the single (non-bootstrapped) ALE curves above. The intervals show how much the estimated effects vary when the training data are resampled while keeping the fitted forest fixed. Use the bootstrapped results to assess the reliability of the patterns rather than relying on a single curve. By default, requesting bootstrapping while retaining ALE statistics also asks `ALE()` to generate fast surrogate p-values. These are useful for exploratory analysis but should not be used for definitive or published conclusions. See `help(ALEpDist)` and the [ALE statistics vignette](ale-statistics.html) for more rigorous options. ## ALE interactions Another advantage of ALE is that it provides data for 2D interactions between variables. Because the bootstrapped object contains the selected `carat:clarity` and `cut:table` interactions as well as all 1D effects, the same `ALEPlots` object provides both kinds of plots. We can use `subset()` to select just the interactions for a separate grid: ```{r print-interactions, fig.width=12, fig.height=5} ale_rf_diamonds_boot |> plot() |> subset(list(d2 = TRUE)) |> print(ncol = 2) ``` The `get()` method accepts either a character term or standard R formula notation: ```{r print-specific-ixn, fig.width=6, fig.height=4} ale_rf_diamonds_boot |> plot() |> get("carat:clarity") ``` The interaction plot shows the joint effect beyond the separate 1D effects. Grey values fall in the middle range of ALE values; colours away from grey indicate an additional joint effect. The sizes of the squares indicate the relative percentage of actual data in each interaction intersection, helping us assess where the data support the estimated interaction. Note that ALE interactions are very particular: an ALE interaction means that two variables have a composite effect over and above their separate independent effects. Two predictors can each have strong one-way effects without their combination having an additional interaction effect. ## Custom prediction functions `ALE()` normally calls a model's prediction method and handles ranger's prediction format automatically. To demonstrate a custom `pred_fun`, we explicitly call ranger's prediction method and extract its numeric predictions. The function must accept `object`, `newdata`, and `type`, and return a numeric vector or matrix. Explicitly invoking ranger's method also makes this custom function available during parallel processing: ```{r ale_custom, eval=FALSE, echo=TRUE, fig.width=12, fig.height=22} custom_predict <- function(object, newdata, type = 'response') { ranger:::predict.ranger( # must explicitly invoke ranger for parallelization to work object = object, data = newdata, type = type )$predictions } # To run the slow code yourself, uncomment and execute the following call. # ale_rf_diamonds_custom <- ALE( # rf_diamonds, # data = diamonds, # pred_fun = custom_predict # ) ale_rf_diamonds_custom <- serialized_objects_site |> file.path("ale_rf_diamonds_custom.0.5.3.rds") |> url() |> readRDS() plot(ale_rf_diamonds_custom) |> print(ncol = 2) ``` ```{r ale_custom-rds, eval=TRUE, echo=FALSE, fig.width=12, fig.height=22} custom_predict <- function(object, newdata, type = 'response') { ranger:::predict.ranger( # must explicitly invoke ranger for parallelization to work object = object, data = newdata, type = type )$predictions } ale_rf_diamonds_custom <- retrieve_rds( # For speed, load a pre-created object by default. if (file.exists("../download/ale_rf_diamonds_custom.0.5.3.rds")) { readRDS("../download/ale_rf_diamonds_custom.0.5.3.rds") } else { c(serialized_objects_site, "ale_rf_diamonds_custom.0.5.3.rds") }, { ALE( rf_diamonds, data = diamonds, pred_fun = custom_predict ) } ) # saveRDS(ale_rf_diamonds_custom, "../download/ale_rf_diamonds_custom.0.5.3.rds") plot(ale_rf_diamonds_custom) |> print(ncol = 2) ``` A custom function is especially useful for model classes whose prediction method has a non-standard signature or returns a structured object. ## Retrieving ALE data and statistics Plots are useful for interpretation, but the `get()` method retrieves the underlying data for further analysis. To demonstrate all of its principal outputs, we create a smaller selection of 1D and 2D effects, run ten bootstrap iterations, and retain the otherwise-discarded iteration-level data with `output_boot_data = TRUE`: ```{r ale_boot_data, eval=FALSE, echo=TRUE} # To run the slow code yourself, uncomment and execute the following call. # ale_diamonds_with_boot_data <- ALE( # rf_diamonds, # data = diamonds, # x_cols = ~ carat + cut + clarity + carat:clarity + cut:table, # output_boot_data = TRUE, # boot_it = 10 # ) ale_diamonds_with_boot_data <- serialized_objects_site |> file.path("ale_rf_diamonds_with_boot_data.0.5.3.rds") |> url() |> readRDS() ``` ```{r ale_boot_data-rds, eval=TRUE, echo=FALSE} ale_diamonds_with_boot_data <- retrieve_rds( # For speed, load a pre-created object by default. if (file.exists("../download/ale_rf_diamonds_with_boot_data.0.5.3.rds")) { readRDS("../download/ale_rf_diamonds_with_boot_data.0.5.3.rds") } else { c(serialized_objects_site, "ale_rf_diamonds_with_boot_data.0.5.3.rds") }, { ALE( rf_diamonds, data = diamonds, x_cols = ~ carat + cut + clarity + carat:clarity + cut:table, output_boot_data = TRUE, boot_it = 10 ) } ) # saveRDS(ale_diamonds_with_boot_data, "../download/ale_rf_diamonds_with_boot_data.0.5.3.rds") ``` By default, `get()` returns ALE data. A formula can limit the result to selected 1D and 2D terms: ```{r get_ale_data} get(ale_diamonds_with_boot_data, ~ carat + cut:table) ``` Because this object was constructed with `output_boot_data = TRUE`, the raw results from individual bootstrap iterations are also available: ```{r get_boot_data} get(ale_diamonds_with_boot_data, what = "boot_data") ``` The `stats` argument retrieves different summaries of the effects: ```{r get_ale_stats} # Estimates for every available statistic get(ale_diamonds_with_boot_data, stats = "estimate") # Complete results for selected statistics get(ale_diamonds_with_boot_data, stats = c("aled", "naler")) # Complete results for all available statistics get(ale_diamonds_with_boot_data, stats = "all") # All confidence regions, or only statistically significant regions get(ale_diamonds_with_boot_data, stats = "conf_regions") get(ale_diamonds_with_boot_data, stats = "conf_sig") ``` These calls return data rather than plots: `"estimate"` provides a compact overview, named statistics return their estimates and inferential details, and the confidence-region options identify portions of an ALE curve or surface that exceed the relevant threshold. See `help("get-ALE-method")` for the precise return structures and the [ALE statistics vignette](ale-statistics.html) before drawing inferential conclusions. ## Where to go next This introduction has covered the main `ALE()` workflow: default effects, selected 1D and 2D terms, bootstrapping, custom prediction functions, plotting, and data extraction. The other vignettes cover [ALE-based statistical inference](ale-statistics.html), [different predictor datatypes](ale-x-datatypes.html), and [small datasets and full-model bootstrapping](ale-small-datasets.html). For exhaustive argument and selection details, see `help(ALE)`, `help(resolve_x_cols)`, and `help("get-ALE-method")`.