ALE-based statistics for statistical inference and effect sizes

Create p-value distribution objects

Before we bootstrap the model to create ALE and other data, there is an important preliminary step when our goal is to analyze ALE statistics. As with any statistics calculated on a dataset, there is some randomness to the statistic values that the procedure will give us. To quantify this randomness, we want to obtain p-values for these statistics. P-values for most standard statistics are based on the assumption that these statistics fit some distribution or another (e.g., Student’s t, χ², etc.). With these distributional assumptions, p-values can be calculated very quickly. However, a key characteristic of ALE is that there are no distributional assumptions: ALE data is a description of a model’s characterization of the data given to it. Accordingly, ALE statistics do not assume any distribution, either.

The implication of this for p-values is that the distribution of the data must be discovered by simulation rather than calculated based on distributional assumptions. The procedure for calculating p-values is the following:

• A random variable is added to the dataset.
• The model is retrained with all variables including this new random variable.
• ALE statistics are calculated for the random variable.
• This procedure is repeated 1,000 times or so to get 1,000 statistic values from the 1,000 random variables.
• p-values are calculated based on the frequency of times that the random variables obtain specific statistic values.

As you can imagine, this procedure is very slow: it involves retraining the entire model on the full dataset 1,000 times. The {ale} package speeds of the process significantly through parallel processing (implemented by default), but it still involves the speed of retraining the model hundreds of times.

To avoid having to repeat this procedure several times (as would be the case when you are doing exploratory analyses), the create_p_dist() function generates a p_dist object that can be run once for a given model-dataset pair. The p_dist object contains functions that can generate p-values based on statistics for any variable from the same model-dataset pair. It generates these p-values when passed to the ale() or model_bootstrap() functions. For very large datasets, the process of generating the p_dist object could be sped up by only using a subset of the data and by running fewer than 1,000 random iterations by setting the rand_it argument. However, the create_p_dist() function will not allow fewer than 100 iterations, otherwise the p-values thus generated would be meaningless.)

We now demonstrate how to create the p_dist object for our case.

# # To generate the code, uncomment the following lines.
# # But it is slow because it retrains the model 1000 times, so this vignette loads a pre-created ale_p object.
# gam_math_p_dist <- create_p_dist(
#   math,
#   gam_math
# )
# saveRDS(gam_math_p_dist, file.choose())
gam_math_p_dist <- url('https://github.com/tripartio/ale/raw/main/download/gam_math_p_dist.rds') |> 
  readRDS()

We can now proceed to bootstrap the model for ALE analysis.

Bootstrap the model with p-values

By default, model_bootstrap() runs 100 bootstrap iterations; this can be controlled with the boot_it argument. Bootstrapping is usually rather slow, even on small datasets, since the entire process is repeated that many times. The model_bootstrap() function speeds of the process significantly through parallel processing (implemented by default), but it still involves retraining the entire model dozens of times. The default of 100 should be sufficiently stable for model building, when you would want to run the bootstrapped algorithm several times and you do not want it to be too slow each time. For definitive conclusions, you could run 1,000 bootstraps or more to confirm the results of 100 bootstraps.

# # To generate the code, uncomment the following lines.
# # But bootstrapping is slow because it retrains the model, so this vignette loads a pre-created ale_boot object.
# mb_gam_math <- model_bootstrap(
#   math,
#   gam_math,
#   # Pass the p_dist object so that p-values will be generated
#   ale_options = list(p_values = gam_math_p_dist),
#   # For the GAM model coefficients, show details of all variables, parametric or not
#   tidy_options = list(parametric = TRUE),
#   # tidy_options = list(parametric = NULL),
#   boot_it = 100  # default
# )
# saveRDS(mb_gam_math, file.choose())
mb_gam_math <- url('https://github.com/tripartio/ale/raw/main/download/mb_gam_math_stats_vignette.rds') |> 
  readRDS()

We can see the bootstrapped values of various overall model statistics by printing the model_stats element of the model bootstrap object:

mb_gam_math$model_stats
#> # A tibble: 9 × 7
#>   name          boot_valid conf.low  median    mean conf.high     sd
#>   <chr>              <dbl>    <dbl>   <dbl>   <dbl>     <dbl>  <dbl>
#> 1 df                NA       29.5    42.3    42.6      58.0   7.78  
#> 2 df.residual       NA      102.    118.    117.      131.    7.78  
#> 3 nobs              NA      160     160     160       160     0     
#> 4 adj.r.squared     NA        0.844   0.896   0.895     0.938 0.0249
#> 5 npar              NA       66      66      66        66     0     
#> 6 mae                1.35     1.25   NA      NA         2.13  0.211 
#> 7 sa_mae_mad         0.726    0.541  NA      NA         0.754 0.0512
#> 8 rmse               1.71     1.58   NA      NA         2.77  0.294 
#> 9 sa_rmse_sd         0.724    0.574  NA      NA         0.748 0.0528

The names of the columns follow the broom package conventions:

• name is the specific overall model statistic described in the row.
• estimate is the bootstrapped estimate of the statistic. It is the same as the bootstrap mean by default, though it can be set to the median with the boot_centre argument of model_bootstrap(). Regardless, both the mean and median estimates are always returned. The estimate column is provided for convenience since that is a standard name in the broom package.
• conf.low and conf.high are the lower and upper confidence intervals respectively. model_bootstrap() defaults to a 95% confidence interval; this can be changed by setting the boot_alpha argument (the default is 0.05 for a 95% confidence interval).
• sd is the standard deviation of the bootstrapped estimate.

Our focus, however, in this vignette is on the effects of individual variables. These are available in the model_coefs element of the model bootstrap object:

mb_gam_math$model_coefs
#> # A tibble: 3 × 6
#>   term              conf.low median   mean conf.high std.error
#>   <chr>                <dbl>  <dbl>  <dbl>     <dbl>     <dbl>
#> 1 (Intercept)         11.7   12.7   12.7      13.6       0.484
#> 2 publicTRUE          -2.02  -0.652 -0.689     0.415     0.637
#> 3 high_minorityTRUE   -0.318  1.05   1.03      2.32      0.676

In this vignette, we cannot go into the details of how GAM models work (you can learn more with Noam Ross’s excellent tutorial). However, for our model illustration here, the estimates for the parametric variables (the non-numeric ones in our model) are interpreted as regular statistical regression coefficients whereas the estimates for the non-parametric smoothed variables (those whose variable names are encapsulated by the smooth s() function) are actually estimates for expected degrees of freedom (EDF in GAM). The smooth function s() lets GAM model these numeric variables as flexible curves that fit the data better than a straight line. The estimate values for the smooth variables above are not so straightforward to interpret, but suffice it to say that they are completely different from regular regression coefficients.

The ale package uses bootstrap-based confidence intervals, not p-values that assume predetermined distributions, to determine statistical significance. Although they are not quite as simple to interpret as counting the number of stars next to a p-value, they are not that complicated, either. Based on the default 95% confidence intervals, a coefficient is statistically significant if conf.low and conf.high are both positive or both negative. We can filter the results on this criterion:

mb_gam_math$model_coefs |> 
  # filter is TRUE if conf.low and conf.high are both positive or both negative because
  # multiplying two numbers of the same sign results in a positive number.
  filter((conf.low * conf.high) > 0)
#> # A tibble: 1 × 6
#>   term        conf.low median  mean conf.high std.error
#>   <chr>          <dbl>  <dbl> <dbl>     <dbl>     <dbl>
#> 1 (Intercept)     11.7   12.7  12.7      13.6     0.484

The statistical significance of the estimate (EDF) of the smooth terms is meaningless here because EDF cannot go below 1.0. Thus, even the random term s(rand_norm) appears to be “statistically significant”. Only the values for the non-smooth (parametric terms) public and high_minority should be considered here. So, we find that neither of the coefficient estimates of public nor of high_minority has an effect that is statistically significantly different from zero. (The intercept is not conceptually meaningful here; it is a statistical artifact.)

This initial analysis highlights two limitations of classical hypothesis-testing analysis. First, it might work suitably well when we use models that have traditional linear regression coefficients. But once we use more advanced models like GAM that flexibly fit the data, we cannot interpret coefficients meaningfully and so it is not so clear how to reach inferential conclusions. Second, a basic challenge with models that are based on the general linear model (including GAM and almost all other statistical analyses) is that their coefficient significance compares the estimates with the null hypothesis that there is no effect. However, even if there is an effect, it might not be practically meaningful. As we will see, ALE-based statistics are explicitly tailored to emphasize practical implications beyond the notion of “statistical significance”.

ALE plots with p-values

mb_gam_plots <- plot(mb_gam_math)
mb_gam_1D_plots <- mb_gam_plots$distinct$math_avg$plots[[1]]
patchwork::wrap_plots(mb_gam_1D_plots, ncol = 2)

We can see that most variables seem to have some sort of mean effect across various values. However, for statistical inference, our focus must be on the bootstrap intervals. Crucial to our interpretation is the middle grey band that indicates the median ± 5% of random values. Below, we will explain what exactly the ALE range (ALER) means, but for now, we can say this:

• The approximate middle of the grey band is the median of the y outcome variables in the dataset (math_avg, in our case). The middle tick on the right y axis indicates the exact median. (The plot() function lets you centre the data on the mean or on zero if you prefer with the relative_y argument.)
• We call this grey band the “ALER band”. 95% of random variables had ALE values that fully lay within the ALER band.
• The dashed lines above and below the ALER band expand the boundaries to where 99% of the random variables were constrained. These boundaries could be considered as demarcating an extended or outward ALER band.

The idea is that if the ALE values of any predictor variable falls fully within the ALER band, then it has no greater effect than 95% of purely random variables. Moreover, to consider any effect in the ALE plot to be statistically significant (that is, non-random), there should be no overlap between the bootstrapped confidence regions of a predictor variable and the ALER band. (For the threshold p-values, We use the conventional defaults of 0.05 for 95% confidence and 0.01 for 99% confidence, but the value can be changed with the p_alpha argument.)

For categorical variables (public and high_minority above), the confidence interval bars for all categories overlap the ALER band. The confidence interval bars indicate two useful pieces of information to us. When we compare them to the ALER band, their overlap or lack thereof tells us about the practical significance of the category. When we compare the confidence bars of one category with those of others, it allows us to assess if the category has a statistically significant effect that is different from that of the other categories; this is equivalent to the regular interpretation of coefficients for GAM and other GLM models. In both cases, the confidence interval bars of the TRUE and FALSE categories overlap each other, indicating that there is no statistically significant difference between categories. Whereas the coefficient table above based on classic statistics indicated this conclusion for public, it indicated that high_minority had a statistically significant effect; our ALE analysis indicates that high_minority does not. In addition, each confidence interval band overlaps the ALER band, indicating that none of the effects is meaningfully different from random results, either.

For numeric variables, the confidence regions overlap the ALER band for most of the domains of the predictor variables except for some regions that we will examine. The extreme points of each variable (except for discrim and female_ratio) are usually either slightly below or slightly above the ALER band, indicating that extreme values have the most extreme effects: math achievement increases with increasing school size, academic track ratio, and mean socioeconomic status, whereas it decreases with increasing minority ratio. The ratio of females and the discrimination climate both overlap the ALER band for the entirety of their domains, so any apparent trends are not supported by the data.

Of particular interest is the random variable rand_norm, whose average ALE appears to show some sort of pattern. However, we note that the 95% confidence intervals we use mean that if we were to retry the analysis for twenty different random seeds, we would expect at least one of the random variables to partially escape the bounds of the ALER band. We will return below to the implications of random variables in ALE analysis.

ALE plots without p-values

Before we continue, let us take a brief detour to see what we get if we run model_bootstrap() without passing it a p_dist object. This might be because we forget to do so, or because we want to see quick results without the slow process of first generating a p_dist object. Let us run model_bootstrap() again, but this time, without p-values.

# # To generate the code, uncomment the following lines.
# # But bootstrapping is slow because it retrains the model, so this vignette loads a pre-created ale_boot object.
# mb_gam_no_p <- model_bootstrap(
#   math, 
#   gam_math,
#   # For the GAM model coefficients, show details of all variables, parametric or not
#   tidy_options = list(parametric = TRUE),
#   # tidy_options = list(parametric = NULL),
#   boot_it = 40  # 100 by default but reduced here for a faster demonstration
# )
# saveRDS(mb_gam_no_p, file.choose())
mb_gam_no_p <- url('https://github.com/tripartio/ale/raw/main/download/mb_gam_no_p_stats_vignette.rds') |> 
  readRDS()

mb_no_p_plots <- plot(mb_gam_no_p)
mb_no_p_1D_plots <- mb_no_p_plots$distinct$math_avg$plots[[1]]
patchwork::wrap_plots(mb_no_p_1D_plots, ncol = 2)

In the absence of p-values, the {ale} packages uses alternate visualizations to offer meaningful results, but with somewhat different interpretations of the middle grey band. Without p-values, we do not have any point of reference for the ALER statistics, so we use percentiles around the median as the reference. The middle grey band here indicates the median ± 2.5%, that is, the middle 5% of all average mathematics achievement scores (math_avg) values in the dataset. We call this the “median band”. The idea is that if any predictor can do no better than influencing math_avg to fall within this middle median band, then it only has a minimal effect. For an effect to be considered statistically significant, there should be no overlap between the confidence regions of a predictor variable and the median band. (We use 5% around the median by default, but the value can be changed with the median_band_pct argument.) For further reference, the outer dashed lines indicate the interquartile range of the outcome values, that is, the 25th and 75th percentiles.

We can see that in this case, the 5% median band is much narrower than the 5% ALER band when p-values are calculated, though they might be more similar for a different dataset. This should give us pause to skipping the calculation of p-values, since we might be overly lax in interpreting apparent relationships as meaningful whereas the ALER band indicates that they might not be that different from what random variables might produce.

For most of the rest of this article, we will only analyze the results with ALER bands generated from p-values, though we will briefly revisit median bands without p-values.

ALE effect size measures on the scale of the y outcome variable

# Create object for convenient access to the relevant stats
mb_gam_math_stats <- mb_gam_math$ale$boot$distinct$math_avg$stats[[1]]

Although ALE plots allow rapid and intuitive conclusions for statistical inference, it is often helpful to have summary numbers that quantify the average strengths of the effects of a variable. Thus, we have developed a collection of effect size measures based on ALE tailored for intuitive interpretation. To understand the intuition underlying the various ALE effect size measures, it is useful to first examine the ALE effects plot, which graphically summarizes the effect sizes of all the variables in the ALE analysis. This is generated when ale is executed and both statistics and plots are requested (which is the case by default) and is accessible with the To focus on all the measures for a specific variable, we can access the ale$stats$effects_plot element:

# mb_gam_plots$distinct$math_avg$stats$effects
mb_gam_math |> 
  plot(type = 'effects')
#> $math_avg

#> 
#> attr(,"class")
#> [1] "ale_eff_plot"

This plot is unusual, so it requires some explanation:

• The y (vertical) axis displays the x variables, rather than the x axis. This is consistent with most effect size plots because they list the full names of variables. It is more readable to list them as labels on the y axis than the other way around.
• The x (horizontal) axis thus displays the y (outcome) variable. But there are two representations of this same axis, one at the bottom and one at the top.
• On the bottom is a more typical axis of the outcome variable, in our case, math_avg. It is scaled as expected. In our case, the axis breaks default to five units each from 5 to 20, evenly spaced.
• On the top, the outcome variable is expressed as percentiles ranging from 0 (the minimum outcome value in the dataset) to 100 (the maximum). It is divided into 10 deciles of 10% each. Because percentiles are usually not evenly distributed in a dataset, the decile breaks are not evenly spaced.
• Thus, this plot has two x axes, the lower one in units of the outcome variable and the upper one in percentiles of the outcome variable. To reduce the confusion, the major vertical gridlines that are slightly darker align with the units of the outcome (lower axis) and the minor vertical gridlines that are slightly lighter align with the percentiles (upper axis).
• The vertical grey band in the middle is the NALED band. Its width is the 0.05 p_value of the NALED (explained below). That is, 95% of random variables had a NALED equal or smaller than that width.
• The variables on the horizontal axis are sorted by decreasing ALED and then NALED value (explained below).

Although it is somewhat confusing to have two axes, the percentiles are a direct transformation of the raw outcome values. The first two base ALE effect size measures below are in units of the outcome variable while their normalized versions are in percentiles of the outcome. Thus, the same plot can display the two kinds of measures simultaneously. Referring to this plot can help understand each of the measures, which we proceed to explain in detail.

Before we explain these measures in detail, we must reiterate the timeless reminder that correlation is not causation. So, none of the scores necessarily means that an x variable causes a certain effect on the y outcome; we can only say that the ALE effect size measures indicate associated or related variations between the two variables.

mb_gam_1D_plots$public

mb_gam_math_stats$by_term$public
#> # A tibble: 6 × 7
#>   statistic estimate p.value conf.low median   mean conf.high
#>   <chr>        <dbl>   <dbl>    <dbl>  <dbl>  <dbl>     <dbl>
#> 1 aled         0.375     0     0.0184  0.332  0.375    0.968 
#> 2 aler_min    -0.344     0.2  -0.846  -0.320 -0.344   -0.0201
#> 3 aler_max     0.421     0     0.0174  0.369  0.421    1.20  
#> 4 naled        5.03      0     0       4.72   5.03    11.7   
#> 5 naler_min   -4.26      0.4 -11.2    -3.41  -4.26     0     
#> 6 naler_max    6.07      0.4   0       5.62   6.07    17.8

mb_gam_1D_plots$academic_ratio

mb_gam_math_stats$by_term$academic_ratio
#> # A tibble: 6 × 7
#>   statistic estimate p.value conf.low  median    mean conf.high
#>   <chr>        <dbl>   <dbl>    <dbl>   <dbl>   <dbl>     <dbl>
#> 1 aled         0.708       0    0.350   0.693   0.708     1.21 
#> 2 aler_min    -4.18        0   -8.10   -4.51   -4.18     -0.581
#> 3 aler_max     1.99        0    0.862   1.97    1.99      3.32 
#> 4 naled        8.93        0    4.18    8.71    8.93     14.9  
#> 5 naler_min  -33.4         0  -48.8   -39.4   -33.4      -5.23 
#> 6 naler_max   27.0         0   10.9    27.2    27.0      39.6

Normalized ALE effect size measures

Since ALER and ALED scores are scaled on the range of y for a given dataset, these scores cannot be compared across datasets. Thus, we present normalized versions of each with intuitive, comparable values. For intuitive interpretation, we normalize the scores on the minimum, median, and maximum of any dataset. In principle, we divide the zero-centred y values in a dataset into two halves: the lower half from the 0th to the 50th percentile (the median) and the upper half from the 50th to the 100th percentile. (Note that the median is included in both halves). With zero-centred ALE y values, all negative and zero values are converted to their percentile score relative to the lower half of the original y values while all positive ALE y values are converted to their percentile score relative to the upper half. (Technically, this percentile assignment is called the empirical cumulative distribution function (ECDF) of each half.) Each half is then divided by two to scale them from 0 to 50 so that together they can represent 100 percentiles. (Note: when a centred ALE y value of exactly 0 occurs, we choose to include the score of zero ALE y in the lower half because it is analogous to the 50th percentile of all values, which more intuitively belongs in the lower half of 100 percentiles.) The transformed maximum ALE y is then scaled as a percentile from 0 to 100%.

There is a notable complication. This normalization smoothly distributes ALE y values when there are many distinct values, but when there are only a few distinct ALE y values, then even a minimal ALE y deviation can have a relatively large percentile difference. If any ALE y value is less than the difference between the median in the data and the value either immediately below or above the median, we consider that it has virtually no effect. Thus, the normalization sets such minimal ALE y values as zero.

Its formula is:

$$ norm\_ale\_y = 100 \times \begin{cases} 0 & \text{if } \max(centred\_y < 0) \leq ale\_y \leq \min(centred\_y > 0), \\ \frac{-ECDF_{y_{\leq 0}}(ale\_y)}{2} & \text{if }ale\_y < 0 \\ \frac{ECDF_{y_{\geq 0}}(ale\_y)}{2} & \text{if }ale\_y > 0 \\ \end{cases} $$ where - centred_y is the vector of y values centred on the median (that is, the median is subtracted from all values). - ECDF_y ≥ 0 is the ECDF of the non-negative values in y. - −ECDF_y ≤ 0 is the ECDF of the negative values in y after they have been inverted (multiplied by -1).

Of course, the formula could be simplified by multiplying by 50 instead of by 100 and not dividing the ECDFs by two each. But we prefer the form we have given because it is explicit that each ECDF represents only half the percentile range and that the result is scored to 100 percentiles.

The median band and random variables

When we do not have p-values, the NALED is particularly helpful in comparing the practical relevance of variables against our threshold for the median band by which we consider that a variable needs to shift the outcome on average by more than 5% of the median values. This threshold is the same scale as the NALED. So, we can tell that public school status with its NALED of 5 just barely crosses our threshold.

It is particularly striking to note the ALE effect size measures for the random rand_norm:

mb_gam_1D_plots$rand_norm

mb_gam_math_stats$by_term$rand_norm
#> # A tibble: 6 × 7
#>   statistic estimate p.value conf.low  median    mean conf.high
#>   <chr>        <dbl>   <dbl>    <dbl>   <dbl>   <dbl>     <dbl>
#> 1 aled         0.358       0    0.108   0.351   0.358     0.577
#> 2 aler_min    -1.32        0   -3.33   -1.18   -1.32     -0.167
#> 3 aler_max     1.62        0    0.315   1.64    1.62      2.81 
#> 4 naled        4.80        0    1.43    4.58    4.80      8.20 
#> 5 naler_min  -14.3         0  -33.4   -13.8   -14.3      -3.75 
#> 6 naler_max   22.8         0    3.12   24.5    22.8      37.8

rand_norm has a NALED of 4.8. It might be surprising that a purely random value has any “effect size” to speak of, but statistically, it must have some numeric value or the other. However, by setting our default value for the median band at 5%, we effectively exclude rand_norm from serious consideration. In informal tests with several different random seeds, the random variables never exceeded this 5% threshold. Setting the median band too low at a value like 1% would not have excluded the random variable, but 5% seems like a nice balance. Thus, the effect of a variable like the discrimination climate score (discrim, 5) should probably not be considered practically meaningful.

We realize that 5% as a threshold for the median band is rather arbitrary, inspired by traditional α = 0.05 for statistical significance and confidence intervals. A proper analysis should use p-values, as most of this article does. However, our initial analyses here show that 5% seems to be an effective choice for excluding a purely random variable from consideration, even for quick initial analyses.

We return to using p-values for the rest of this article.

Interpretation of normalized ALE effect sizes

Here we summarize some general principles for interpreting normalized ALE effect sizes.

• Normalized ALE deviation (NALED): this is the average variation of ALE effect of the input variable.
  • Values range from 0 to 100%:
    • 0% means no effect at all.
    • 100% means the maximum possible effect any variable could have: for a binary variable, one value (50% of the data) sets the outcome at its minimum value and the other value (the other 50% of the data) sets the outcome at its maximum value.
  • Larger NALED means stronger effects.
• Normalized ALE range (NALER): these are the minimum and maximum effects of any value of the input variable.
  • NALER minimum ranges from –50% to 0%; NALER maximum ranges from 0% to +50%:
  • 0% means no effect at all. It indicates that the only effect of the input variable is to keep the outcome at the median of its range of values.
  • NALER minimum of n means that, regardless of the effect size of NALED, the minimum effect of any input value shifts the outcome by n percentile points of the outcome range. Lower values (closer to –50%) mean a stronger extreme effect.
  • NALER maximum of x means that, regardless of the effect size of NALED, the maximum effect of any input value shifts the outcome by x percentile points of the outcome range. Greater values (closer to +50%) mean a stronger extreme effect.

In general, regardless of the values of ALE statistics, we should always visually inspect the ALE plots to identify and interpret patterns of relationships between inputs and the outcome.

A common question for interpreting effect sizes is, “How strong does an effect need to be to be considered ‘strong’ or ‘weak’?” On one hand, we refuse to offer general guidelines for how “strong” is “strong”. The simple answer is that it depends entirely on the applied context. It is not meaningful to try to propose numerical values for statistics that are supposed to be useful for all applied contexts.

On the other hand, we do consider it very important to delineate the threshold between random effects and non-random effects. It is always important to distinguish between a weak but real effect from one that is just a statistical artifact due to random chance. For that, we can offer some general guidelines based on whether or not we have p-values.

When we have p-values for ALE statistics, then the boundaries of ALER should generally be used to determine the acceptable risk of considering a statistic to be meaningful. Statistically significant ALE effects are those that are less than the 0.05 p_value ALER minimum of a random variable and greater than the 0.05 p_value maximum of a random variable. As we explained above when introducing the ALER band, this is precisely what the {ale} package does, especially in the plots that highlight the ALER band and the confidence region tables that use the specified ALER p_value threshold.

In the absence of p-values, we suggest that NALED can be a general guide for non-random values. In our informal tests, we find that NALED values below 5% have the same average effect as a random variable. That is, the average effect is not reliable; it might be random. However, regardless of the average effect indicated by NALED, large NALER effects indicate that the ALE plot should be inspected to interpret the exceptional cases. This caveat is very important; unlike GLM coefficients, ALE analysis is sensitive to exceptions to the overall trend. This is precisely what makes it valuable for detecting non-linear effects.

In general, if NALED < 5%, NALER minimum > –5%, and NALER maximum < +5%, the input variable has no meaningful effect. All other cases are worth inspecting the ALE plots for careful interpretation: - NALED > 5% means a meaningful average effect. - NALER minimum < –5% means that there might be at least one input value that significantly lowers the outcome values. - NALER maximum > +5% means that there might be at least one input value that significantly increases the outcome values.

ALE data structures for categorical and numeric variables

To understand how bootstrapped ALE can be used for statistical inference, we must understand the structure of ALE data. Let’s begin simple with a binary variable with just two categories, public:

mb_gam_math$ale$boot$distinct$math_avg$ale[[1]]$public
#> # A tibble: 2 × 7
#>   public.bin    .n     .y  .y_lo .y_mean .y_median .y_hi
#>   <ord>      <int>  <dbl>  <dbl>   <dbl>     <dbl> <dbl>
#> 1 FALSE         70  0.383 -0.219   0.383     0.356 1.20 
#> 2 TRUE          90 -0.306 -0.846  -0.306    -0.308 0.196

Here is the meaning of each column of ale$data for a categorical variable:

• .bin: the different categories that exist in the categorical variable.
• .n: the number of rows for that category in the dataset provided to the function.
• ale_y: the ALE function value calculated for that category. For bootstrapped ALE, this is the same as ale_y_mean by default or ale_y_median if the boot_centre = 'median' argument is specified.
• ale_y_lo and ale_y_hi: the lower and upper confidence intervals for the bootstrapped ale_y value.

By default, the ale package centres ALE values on the median of the outcome variable; in our dataset, the median of all the schools’ average mathematics achievement scores is 12.9. With ALE centred on the median, the weighted sum of ALE y values (weighted on .n) above the median is approximately equal to the weighted sum of those below the median. So, in the ALE plots above, when you consider the number of instances indicated by the rug plots and category percentages, the average weighted ALE y approximately equals the median.

Here is the ALE data structure for a numeric variable, academic_ratio:

mb_gam_math$ale$boot$distinct$math_avg$ale[[1]]$academic_ratio
#> # A tibble: 63 × 7
#>    academic_ratio.ceil    .n     .y .y_lo .y_mean .y_median   .y_hi
#>                  <dbl> <int>  <dbl> <dbl>   <dbl>     <dbl>   <dbl>
#>  1                0        1 -5.45  -8.13  -5.45     -5.37  -1.82  
#>  2                0.05     2 -2.46  -3.96  -2.46     -2.60  -0.549 
#>  3                0.09     1 -1.40  -2.68  -1.40     -1.38  -0.435 
#>  4                0.1      2 -1.28  -2.34  -1.28     -1.29  -0.0767
#>  5                0.13     1 -0.613 -1.61  -0.613    -0.600  0.387 
#>  6                0.14     2 -0.755 -1.58  -0.755    -0.852  0.702 
#>  7                0.17     1 -0.700 -1.48  -0.700    -0.746  0.311 
#>  8                0.18     4 -0.573 -1.36  -0.573    -0.592  0.790 
#>  9                0.19     3 -0.513 -1.26  -0.513    -0.502  0.649 
#> 10                0.2      3 -0.530 -1.34  -0.530    -0.539  0.543 
#> # ℹ 53 more rows

The columns are the same as with a categorical variable, but instead of .bin, we have .ceil since there are no categories. To calculate ALE for numeric variables, the range of x values is divided into bins (by default 100, customizable with the max_num_bins argument). Because numeric variables often have multiple values in a bin, the ALE data stores the ceilings (upper bounds) of the bins. If the x values have fewer than 100 distinct values in the data, then each distinct value becomes its own bin and we record its value as the ceiling of the bin. (This is often the case with smaller datasets like ours; here academic_ratio has only distinct values.) If there are more than 100 distinct values, then the range is divided into 100 percentile groups. The other columns mean the same thing as with categorical variables: .n is the number of rows of data in each bin and .y is the calculated ALE for the bin whose ceiling is .ceil.

Bootstrap-based inference with ALE

In a bootstrapped ALE plot, values within the confidence intervals are statistically significant; values outside of the ALER band can be considered at least somewhat meaningful. Thus, the essence of ALE-based statistical inference is that only effects that are simultaneously within the confidence intervals AND outside of the ALER band should be considered conceptually meaningful.

We can see this, for example, with the plot of mean_ses:

mb_gam_1D_plots$mean_ses

It might not always be easy to tell from a plot which regions are relevant, so the results of statistical significance are summarized with the ale$conf_regions$by_term element, which can be accessed for each variable from its by_term element:

mb_gam_math_stats$by_term$mean_ses
#> # A tibble: 6 × 7
#>   statistic estimate p.value conf.low median   mean conf.high
#>   <chr>        <dbl>   <dbl>    <dbl>  <dbl>  <dbl>     <dbl>
#> 1 aled          1.15       0    0.760   1.15   1.15      1.55
#> 2 aler_min     -6.98       0  -10.2    -7.22  -6.98     -3.34
#> 3 aler_max      2.79       0    1.46    2.74   2.79      4.86
#> 4 naled        13.6        0    9.47   13.7   13.6      17.9 
#> 5 naler_min   -45.9        0  -50     -47.5  -45.9     -33.1 
#> 6 naler_max    34.2        0   18.0    35.6   34.2      44.4

For numeric variables, the confidence regions summary has one row for each consecutive sequence of x values that have the same status: all values in the region are below the middle irrelevance band, they overlap the band, or they are all above the band. Here are the summary components:

• start_x is the first and end_x is the last x value in the sequence. start_y is the y value that corresponds to start_x while end_y corresponds to end_x.
• n is the number of data elements in the sequence; pct is the percentage of total data elements out of the total number.
• x_span is the length of x of the sequence that has the same confidence status. However, so that it may be comparable across variables with different units of x, x_span is expressed as a percentage of the full domain of x values.
• trend is the average slope from the point (start_x, start_y) to (end_x, end_y). Because only the start and end points are used to calculate trend, it does not reflect any ups and downs that might occur between those two points. Since the various x values in a dataset are on different scales, the scales of the x and y values in calculating the trend are normalized on a scale of 100 each so that the trends for all variables are directly comparable. A positive trend means that, on average, y increases with x; a negative trend means that, on average, y decreases with x; a zero trend means that y has the same value at its start and end points–this is always the case if there is only one point in the indicated sequence.
• mid_bar is the key information here. It indicates if all the values in sequence from start_x to end_x are below, overlapping, or above the middle ALER band:
  • below: the higher limit of the confidence interval of ALE y (ale_y_hi) is below the lower limit of the ALER band.
  • above: the lower limit of the confidence interval of ALE y (ale_y_lo) is above the higher limit of the ALER band.
  • overlap: neither of the first two conditions holds; that is, the confidence region from ale_y_lo to ale_y_hi at least partially overlaps the ALER band.

These results tell us that, for mean_ses, from -1.19 to -1.04, ALE is below the median band from 6.1 to 7.6. From -0.792 to -0.792, ALE overlaps the median band from 10.2 to 10.2. From -0.756 to -0.674, ALE is below the median band from 10.2 to 10.8. From -0.663 to -0.663, ALE overlaps the median band from 10.9 to 10.9. From -0.643 to -0.484, ALE is below the median band from 10.8 to 11.2. From -0.467 to -0.467, ALE overlaps the median band from 11.5 to 11.5. From -0.46 to -0.46, ALE is below the median band from 11.4 to 11.4. A few other regions briefly exceeded the ALER band.-

Interestingly, most of the text of the previous paragraph was generated automatically by an internal (unexported function) ale:::summarize_conf_regions_1D_in_words. (Since the function is not exported, you must use ale::: with three colons, not just two, if you want to access it.)

mb_gam_math_stats$conf_regions$by_term |> 
  filter(term == 'mean_ses') |> 
  ale:::summarize_conf_regions_1D_in_words()
#> [1] "From -1.19 to -0.368, ALE is below the median band from -7.83 to -1.06. From -0.347 to 0.163, ALE overlaps the median band from -0.829 to 0.972. From 0.179 to 0.179, ALE is above the median band from 1.11 to 1.11. From 0.188 to 0.188, ALE overlaps the median band from 1.07 to 1.07. From 0.218 to 0.312, ALE is above the median band from 1.23 to 0.981. From 0.316 to 0.316, ALE overlaps the median band from 0.852 to 0.852. From 0.333 to 0.333, ALE is above the median band from 0.876 to 0.876. From 0.334 to 0.535, ALE overlaps the median band from 0.796 to 1.12. From 0.569 to 0.759, ALE is above the median band from 1.51 to 2.85. From 0.831 to 0.831, ALE overlaps the median band from 1.86 to 1.86."

While the wording is rather mechanical, it nonetheless illustrates the potential value of being able to summarize the inferentially relevant conclusions in tabular form.

Confidence region summary tables are available not only for numeric but also for categorical variables, as we see with public. Here is its ALE plot again:

mb_gam_1D_plots$public

And here is its confidence regions summary table:

mb_gam_math_stats$conf_regions$by_term |> 
  filter(term == 'public')
#> # A tibble: 2 × 12
#>   term   x     start_x end_x x_span_pct     n   pct      y start_y end_y trend
#>   <chr>  <chr>   <dbl> <dbl>      <dbl> <int> <dbl>  <dbl>   <dbl> <dbl> <dbl>
#> 1 public FALSE      NA    NA         NA    70  43.8  0.383      NA    NA    NA
#> 2 public TRUE       NA    NA         NA    90  56.2 -0.306      NA    NA    NA
#> # ℹ 1 more variable: mid_bar <ord>

Since we have categories here, there is no start or end positions and there is no trend. We instead have each x category and its single ALE y value, with the n and pct of the respective category and mid_bar as before to indicate whether the indicated category is below, overlaps with, or is above the ALER band.

Again with the help of ale:::summarize_conf_regions_1D_in_words(), these results tell us that, for public, for FALSE, the ALE of 13.3 overlaps the ALER band. For TRUE, the ALE of 12.6 overlaps the ALER band.

Again, our random variable rand_norm is particularly interesting. Here is its ALE plot:

mb_gam_1D_plots$rand_norm

And here is its confidence regions summary table:

mb_gam_math_stats$conf_regions$by_term |> 
  filter(term == 'rand_norm')
#> # A tibble: 1 × 12
#>   term    x     start_x end_x x_span_pct     n   pct     y start_y  end_y  trend
#>   <chr>   <chr>   <dbl> <dbl>      <dbl> <int> <dbl> <dbl>   <dbl>  <dbl>  <dbl>
#> 1 rand_n… <NA>    -2.40  2.61        100   160   100    NA   -1.35 -0.416 0.0626
#> # ℹ 1 more variable: mid_bar <ord>

Despite any apparent pattern, we see that from -2.4 to 2.61, ALE overlaps the median band from 12 to 12.8. So, despite the random highs and lows in the bootstrap confidence interval, there is no reason to suppose that the random variable has any effect anywhere in its domain.

We can conveniently summarize all the confidence regions from all variables that are statistically significant or meaningful by accessing the conf_regions$significant element:

mb_gam_math_stats$conf_regions$significant
#> # A tibble: 12 × 12
#>    term      x      start_x    end_x x_span_pct     n   pct     y start_y  end_y
#>    <chr>     <chr>    <dbl>    <dbl>      <dbl> <int> <dbl> <dbl>   <dbl>  <dbl>
#>  1 size      <NA>  2403      2.40e+3       0        2  1.25    NA   1.45   1.45 
#>  2 size      <NA>  2650      2.65e+3       0        2  1.25    NA   1.87   1.87 
#>  3 academic… <NA>     0      9   e-2       9        4  2.5     NA  -5.45  -1.40 
#>  4 academic… <NA>     0.96   1   e+0       4.00    14  8.75    NA   1.54   1.93 
#>  5 mean_ses  <NA>    -1.19  -3.68e-1      40.6     33 20.6     NA  -7.83  -1.06 
#>  6 mean_ses  <NA>     0.179  1.79e-1       0        2  1.25    NA   1.11   1.11 
#>  7 mean_ses  <NA>     0.218  3.12e-1       4.66    11  6.88    NA   1.23   0.981
#>  8 mean_ses  <NA>     0.333  3.33e-1       0        2  1.25    NA   0.876  0.876
#>  9 mean_ses  <NA>     0.569  7.59e-1       9.41    13  8.12    NA   1.51   2.85 
#> 10 minority… <NA>     0      5.71e-2       5.71    52 32.5     NA   1.63   0.972
#> 11 minority… <NA>     0.409  5.61e-1      15.2     11  6.88    NA  -1.15  -2.11 
#> 12 minority… <NA>     0.955  1   e+0       4.48     9  5.62    NA  -2.20  -3.82 
#> # ℹ 2 more variables: trend <dbl>, mid_bar <ord>

This summary focuses only on the x variables that have meaningful ALE regions anywhere in their domain. We can also conveniently isolate which variables have any such meaningful region by extracting the unique values in the term column:

mb_gam_math_stats$conf_regions$significant$term |> 
  unique()
#> [1] "size"           "academic_ratio" "mean_ses"       "minority_ratio"

This is especially useful for analyses with dozens of variables; we can thus quickly isolate and focus on the most meaningful ones.