Detecting interactions with random forests: a comment on Gries’ words of caution and suggestions for improvement

2.1 From trees to ensemble methods

A problematic property of trees, that is often mentioned in the literature and has led to the development of ensemble methods like random forests, is their instability to small changes in the data. While Gries (2020) mentions that these changes can affect both the prediction and the tree structure, Philipp et al. (2016, 2018 point out that quite different-looking tree structures can actually lead to essentially the same predictions. This highlights that often the tree structure is over- or misinterpreted. For example, while the variable used for the first split in a tree is the one that showed the strongest main effect, it is the entire pattern of the predictions in the end nodes that shows whether and how several splits in the same variable approximate the functional form (as illustrated for a linear effect approximated by several splits in a tree by Gries 2020, Fig. 3) or whether the first variable works together with the variables below to form an interaction effect (as illustrated by Strobl et al. 2009b, Fig. 4).

Ensemble methods like random forests (Breiman 2001) and their predecessor method bagging (Breiman 1996) have overcome the instability issue by means of averaging predictions over several hundreds or thousands of trees. This makes the predictions of ensembles much more smooth and typically more accurate than those of single trees, but comes at the price that one looses the interpretability that is inherent to individual trees (as long as they are not too large, as Gries 2020, points out). This lack of interpretability refers to the fact that random forests do not offer direct insights into the relationship between the individual predictor variables and the response variable. It remains unclear whether and how the individual predictor variables contribute to the prediction of the response variable alone and/or in interactions. For this reason, in the machine learning literature ensemble methods are often termed black box methods (a metaphor for the fact that one enters the predictor variables and out comes the prediction – but what happened in between is obscure).

2.2 Overfitting and related issues

The way that ensembles of bagged trees and random forests are constructed – by fitting individual trees on bootstrap samples with replacement or (preferably for unbiased variable selection, as shown by Strobl et al. 2007b) subsamples without replacement from the original data – also has the advantage that each tree comes with its own, built-in test data set, namely those observations that were not used for training the respective tree (termed out-of-bag, or OOB, data). This is a very useful property of random forests, which Gries (2020) rightly draws attention to. However, since the terminology used by Gries (2020) is not entirely in line with the statistical and machine learning literature, we would like to elaborate on this topic.

To explain why it is important to assess the quality of machine learning methods on data that were not used for fitting the model, let’s first review how models are compared and assessed in classical, parametric statistics. For example, imagine that you want to use a linear regression model to predict a response variable from several potential predictor variables, but do not know in advance which predictors are and are not informative for predicting the response (or whether they work in linear main effects, nonlinear or interaction effects). Imagine you estimated two regression models for predicting the response variable: model 1, containing only predictors A and B, and model 2, containing predictors A, B, C, D, E and F. How can you decide which model is better? While, for example, the simple R² statistic for linear regression tends to increase with the number of added variables, even if spurious effects are added, the adjusted R² (but also F-tests or likelihood-ratio-tests for model comparisons in nested models) will account for the model complexity.^[3] So in the world of parametric statistics, there are clear-cut criteria for deciding whether model 2 is better than model 1, which balance the explained variance against the model complexity. But these clear-cut criteria are based on certain mathematical assumptions.

Machine learning methods come with far fewer assumptions than classical parametrical statistical approaches. This also means that we no longer have the entire toolbox of statistical significance testing and confidence intervals available for assessing machine learning solutions. Without this toolbox, the only information we have about how good one machine learning model performed compared to another one is their prediction accuracy (which for categorical response variables is also termed classification accuracy).

Now, similarly to what would happen if you would compare the un-adjusted R² of two regression models, if you compared the prediction accuracies of two machine learning models (again model 1 containing only predictors A and B, model 2 containing A, B, C, D, E and F) on the original data, the more complex model will typically perform better, just because it is more flexible. You can imagine this when you think about predicting the scores of your students in an exam based on a few truly relevant predictors (how long did they study for the exam – probably helpful for getting a high score; did they attend a party the night before the exam – probably not so much). But your data set may also contain other potential predictor variables (do your students have a cat or dog; what is their shoe size; what is the shoe size of their grandmother; … – you get our point, these variables are not informative for the students’ exam performance). If you would use a regression or machine learning model containing only the two truly influential predictors (study time and party attendance) the model would show a decent prediction accuracy, but not a perfect one. Remaining unexplained differences between the students’ exam scores could be due to other, unobserved characteristics, such as their IQ or previous experience with the course topic, but also how they feel that particular day, if they made any careless errors etc. Still, this model would be very useful, because it does not only decently predict the exam scores of this year’s students, but it will also decently predict the grades of next year’s students, as long as the general mechanisms underlying your course and exam don’t change.

If, however, in addition to the truly informative variables you included several irrelevant variables in your model (such as pet ownership and shoe sizes), what would happen is the following: The model would show a higher prediction accuracy for this year’s students. The reason for this is that by means of the additional irrelevant variables, you can make more fine-grained predictions. In the extreme case you would include so many irrelevant variables that only one particular student is described by each combination. For example there might be only one student who studied dozens of hours, did not attend a party the night before the exam, has two cats, shoe size 38, a grandmother with very large feet, etc. A model that is extremely flexible because it contains so many predictors will learn what exam score this particular student achieved, and predict exactly that same exam score for all students with the same values on these predictors. This is actually a very good strategy for predicting the exam score of this particular student in the data set from this particular year, and the model will achieve a very high prediction accuracy on this year’s data this way. However, we would not expect this prediction to work very well for next year’s students – or generally speaking, another sample from the same population. When a model is too flexible and thus adheres too strongly to random variations in the original sample it does not generalize well to other samples from the same population. This is termed overfitting.

Due to the stabilizing effect of averaging over the individual trees, random forests tend to be less affected by overfitting than other machine learning methods. In Gries (2020) this point is correctly made. However, presumably for didactic reasons, the terminology about predictions used by Gries (2020) is not the one common in the statistics and machine learning literature. To avoid confusion, we would like to point out that: (i) In the statistics and machine learning literature, the term prediction describes the process of entering the predictor variable values of a real or hypothetical observation into a fitted model and receiving the model’s prediction of the response variable for that observation. This can be done for individual or all observations in a data set, and for the same data set that the model was fitted on (termed training data or learning data) or new data (termed test data). (ii) The prediction accuracy is computed by comparing the predicted and true values of the response variable and aggregating over all observations. When the response variable is categorical, the prediction accuracy is also termed classification accuracy. It is typically described by the percentage of observations for which the predicted class is equal to the true class.^[4]

In contrast to these conventions, Gries (2020) calls the accuracy of out-of-bag predictions of a random forest (which are based on different observations than those the forest was fitted to, similar to a fresh test sample) “prediction accuracy”, while he calls the prediction accuracy of a tree on the original learning data “classification accuracy”. We fully agree that the prediction accuracy of out-of-bag predictions of a random forest gives a better estimate of the prediction accuracy that is to be expected for other samples from the same population. At the same time, in order to avoid confusion for readers consulting articles or textbooks from statistics or machine learning, we would like to point out that the contrasting use of the terms “prediction” versus “classification” accuracy in Gries (2020) is not in line with how these terms are used in the statistics and machine learning literature.

2.3 Interpretability and stability

Coming back to the interpretability (or lack thereof) of random forests, we agree with Gries (2021) that when the true pattern in the data is more complex, trying to interpret a random forest by means of fitting an additional single tree to the data will produce an oversimplification. Below, we will discuss more adequate means of interpretable machine learning. However, we would also like to point out that when the true pattern in the data is simple, a single tree may be all that it takes to capture this pattern – and in this case, the interpretability of a single tree is a big advantage over the black box property of a random forest. Whether a specific data set may contain a pattern that is sufficiently described by a single tree can be explored by means of the stability diagnostics available in the stablelearner package in R (Philipp et al. 2016, 2018, 2023). In this framework, the variable and cutpoint selection of the tree fit on the original data is compared to that of trees fit on hundreds of samples from the same data (similar to a random forest, but here the original tree has a special role). If the same variables and similar cutpoint locations as in the original tree are selected by the majority of resampling-based trees, the results can be considered as stable and it is safe to interpret the original tree. If not, the complexity of the pattern in the data will be better captured by an ensemble method, and in that case no single tree is suitable for describing the pattern.

2.4 Tuning and runtime

Random forests are often said to “work well off the shelf”, i.e., with the default settings for the tuning parameters. Still, the user should always double check the default settings of each random forest implementation they are using. In particular, while in the randomForest (Liaw and Wiener 2022) function the default value of the mtry argument is set to the default values suggested by Breiman for the original implementation (for classification to the square root of the number of predictor variables, for regression to one third of the number of predictor variables), in party the mtry argument is an arbitrary fixed number. As correctly highlighted by Gries (2021) for randomForest, it is often worthwhile to tune (i.e., to optimally select) the value of mtry. The same holds for cforest. For cforest from the party package, tuning based, e.g., on cross validation is available in the caret package (Kuhn 2008; Kuhn et al. 2023).

Cross validation means that the sample is randomly split into k parts, of which k − 1 parts are used for fitting the model and the remaining part is held back to assess the prediction accuracy on fresh data. This process is repeated k times until all observations have been used in turn for fitting the model and for assessing the prediction accuracy. The entire procedure is conducted for several different values of the to-be-tuned parameter, and the value with the best cross-validated prediction accuracy is chosen. The number of trees in the forest, ntree, is typically not tuned this way, because, while for mtry values smaller or larger than the optimal value can lead to suboptimal performance, increasing the number of trees will only make the results more accurate and more stable – only the extent of the increase in accuracy flattens after a certain point. However, when the number of trees is chosen too small, one might not yet have reached the point where the results are sufficiently stable to, e.g., interpret random forest variable importance scores. To check whether the number of trees is sufficiently large, Strobl et al. (2009a) recommend to re-fit a random forest and re-compute the variable importances using a different seed (which determines the status of the random number generator in R). If the results vary notably, the number of trees was not yet sufficiently high.

Since for larger data sets runtime can be an issue, it may be worth mentioning that in the partykit package the functions for fitting conditional inference trees and forests were written entirely in R (while in party parts of their implementation were outsourced to C). This makes partykit very flexible, but means that, at the moment, fitting conditional inference trees and forests in the older package party is actually faster than in the newer package partykit. For computing the conditional permutation importance (see also below), the permimp package (Debeer et al. 2021) provides a faster implementation.

2.5 Variable importance and functional form

Random forest variable importance measures can give a first impression which variables were relevant for the prediction. They have been motivated rather heuristically and their absolute values depend on various factors beyond the true effect of a variable, so that they are not directly comparable between different data sets. Variable importance scores should therefore merely be interpreted in a qualitative fashion. For example, when a few variables show much higher importance scores than the rest, it makes sense to concentrate on these variables in future studies.^[5] However, researchers often desire clear-cut decisions which variables are “significantly important”. Rothacher and Strobl (2023) review different heuristics for significance tests for random forest variable importance measures. In particular, they show that the rule of thumb suggested by Strobl et al. (2009b), that has been erroneously communicated by some researchers as if it was a significance test, does not possess the properties of a formal statistical test. This rule of thumb is based on the idea that – in a sparse setting, where the number of potential predictors is large, but only few are truly informative – the importance scores of the non-informative predictors will randomly fluctuate around zero. In this kind of setting, which is common, e.g., in genetics, but probably not so common in linguistics, the absolute value of the strongest negative importance can serve as a rule-of-thumb lower bound, in the sense that variables not even exceeding this bound can be considered as noise variables. It is, however, not a formal significance test, and was not presented as one by Strobl et al. (2009b). When nevertheless used as such, Rothacher and Strobl (2023) show that it has an inflated type I error rate. Interestingly, some other, computationally much more intensive approaches show even worse statistical properties.

In addition to the original permutation variable importance, a conditional permutation variable importance is available (Debeer and Strobl 2020; Strobl et al. 2008). As is very well explained by Gries (2021), the conditional permutation scheme takes correlations between predictor variables into account and thus avoids “exaggerating the importance of predictors that are correlated with other [informative] predictors” (Gries 2021, p. 470, with our addition in square brackets). However, the way the conditioning is implemented at the moment leads to increased computation times and also overcorrects in the sense that the conditional permutation importance of all correlated predictors is lowered compared to that of uncorrelated predictors (Henninger et al. 2023a). While work is being done currently to further investigate and improve the behavior of the conditional permutation importance, we would like to point out that neither the unconditional nor the conditional permutation importance should be considered as “the truth”, but rather as different points on the marginal-to-partial continuum (Debeer and Strobl 2020) and that the most information lies in their comparison rather than either result alone.

Variable importance scores indicate descriptively, which variables – alone or in interactions – contributed to the prediction of the random forest. For a more thorough interpretation it would be desirable to learn more about which variables contribute in main effects versus interactions and also in which functional form. Here we would like to clarify two misconceptions, but first need to review a few important terms and concepts. When you think of a linear or logistic regression model, everything is very tidy, in the sense that there are individual coefficients quantifying the linear, quadratic, or other curvilinear effect of a single variable or the interaction effect of two or more variables. As was mentioned above, in a single tree, the tree structure together with the predictions in the end nodes describe these effects and the approximation of the functional form captured by the tree – at least to the trained eye. In a random forest with hundreds or thousands of trees it is not possible to visually captures these effects. There are, however, additional approaches from the field of interpretable machine learning, that provide additional hints. Partial dependence plots (Apley and Zhu 2020; Henninger et al. 2023a, for an introduction) or other plots of the predicted responses for different values of the predictors, such as those suggested by Szmrecsanyi et al. (2016), Gries (2020) and Hundt et al. (2020), are one possibility to illustrate the functional form of the association between a single predictor variable and the prediction for the response variable, typically averaging over the values of all other predictors. It is important to be aware that the information about the functional form of the association between predictor and response, that is approximated rather smoothly by a random forest and then graphically displayed in a partial dependence plot, is a valuable information, because in exploratory modelling the functional form of the association between predictors and response is learned from the data rather than pre-specified by the user.

Still, a univariate partial dependence plot is an approximation of the functional form of the main effect of only one predictor variable at a time. So while Gries (2021) argues correctly that partial dependence plots contain more information than, say, bivariate χ²-tests between each individual predictor and the response – namely the information on the functional form – his statement on page 467 (“…it also makes no sense to summarize a multifactorial forest with monofactorial tables – this is what we use partial dependence plots for.”) should not be misinterpreted in the sense that univariate partial dependence plots could capture interaction effects.

We would like to illustrate this with a simple example: Figure 1 shows partial dependence plots for two predictors from a random forest. We see that the predicted response is essentially the same over the range of the hypothetical numeric (left) and the two values of the binary (right) predictor variable. If there was a strong main effect, the predicted response values should show a visible monotone or non-monotone pattern, but these two predictors seem to be irrelevant – when looking at one of them at a time, marginalizing over all other predictors. We will come back to this example soon, but first would like to provide a little more detail about how partial dependence plots are constructed.

Figure 1:

Partial dependence plots for two predictor variables, x₁ (left) and x₂ (right).

Partial dependence plots display the average predicted value of the response variable (on the y-axis) as a function of the value of a predictor variable (on the x-axis). In order to create a partial dependence plot, the average predicted value of the response variable needs to be computed for different values of the predictor variable.^[6] In order to compute the average predicted value, first a prediction for each single observation in the dataset needs to made by the fitted random forest (or another machine learning model) for each of the different values of the predictor variable. The predictions for the single observations are then averaged over all observations in a second step.

One observation here corresponds to one row in the dataset. In many machine learning applications, every row in the dataset represents one person, and every column represents one variable, such as age, gender, economic status or wellbeing. Imagine that a random forest is used to predict the value of the response variable wellbeing from the values of the predictor variables age, gender and economic status. The partial dependence plot then has the following rationale: For making the prediction for a single observation, the observed values of all other predictor variables for this observation are used – only for the variable of interest the value is varied. For example, in order to create a partial dependence plot for the predictor variable age, the prediction for a single person will use that person’s observed values for gender and economic status, but try out different values for the person’s age. The random forest’s predicted value of the response variable, wellbeing, is then computed for each value of age. These predictions show the levels of wellbeing predicted by the model for a certain person if we could hold all their other properties constant but vary their age to make them younger or older.

For every value of age, the predictions are then averaged over all persons in the second step to create the average prediction. In this sense, partial dependence plots average or marginalize over all other predictor variables. The average prediction (on the y-axis) is displayed in the partial dependence plot as a function of the different values of age that were tried out for each person (on the x-axis, like for the numeric variable at the x-axis of the left panel in Figure 1). If, for example, the predicted wellbeing averaged over all persons decreased with age, we would see a decreasing shape in the partial dependence plot. This would correspond to a monotone main effect of the variable age. Due to the flexibility of random forests to approximate functional forms, it would also be possible that we find a non-monotone effect, such as a u-shaped effect where wellbeing is higher for middle aged and lower for younger and elderly persons.

In corpus linguistics, the rows of a data set typically do not correspond to persons but to linguistic instances, such as occurrences of a verb in a corpus of written or transcribed text. The partial dependence plot is then created accordingly by varying one property of the instance while leaving the other properties as they were observed, computing the predictions for each instance and then averaging over all instances.

Returning to Figure 1, we saw that the two predictors from this toy example seem to be irrelevant when looking at one predictor at a time, like it is done in a univariate partial dependence plot that marginalizes over all other predictor variables. There are, however, other plots from interpretable machine learning, that can capture the effects of two variables at a time: bivariate partial dependence plots and individual conditional expectation or ICE plots (Goldstein et al. 2015). We will further discuss bivariate partial dependence plots below, but in the following concentrate on introducing ICE plots, because they appear to be less widely known in linguistics.

The main difference between univariate partial dependence plots and ICE plots is that in ICE plots the predictions are not averaged over all observations, but displayed separately for each observation. Considering the above description of how partial dependence plots are constructed, this means that the aggregation in the second step of the procedure is left out, leaving one individual line of predicted values for each observation (i.e., each person or each instance) over the range of the predictor variable. These lines are displayed in Figure 2 for the toy example data already used in Figure 1.

Figure 2:

ICE plot for x₁ (left) and ICE plot for x₁ colored w.r.t. x₂ (right).

ICE plots thus display the predicted value of the response variable (on the y-axis) as a function of the value of a predictor variable (on the x-axis) for each observation separately. Accordingly, the partial dependence plots in Figure 1 are the averages over the individual ICE curves in Figure 2. Through this, ICE plots can display more detailed information than partial dependence plots. For more details on partial dependence and ICE plots see Henninger et al. (2023a). We will now describe how ICE plots may help detect interaction effects.

When we look at the left panel of Figure 2, there seem to be observations for which the effect of the numeric predictor is positive, and others for which it is negative. Moreover, if we color the lines for the individual observations according to the second, binary predictor (see Figure 2, right), we see that the pattern actually corresponds to a perfect interaction between the two predictors. This kind of interaction cannot be detected by looking at one variable at a time (it is the XOR problem mentioned earlier). So univariate partial dependence plots or plots displaying the predicted response only for one predictor at a time are also monofactorial in the sense that they risk to miss interactions.

It is of course also possible to manually include the interaction between two predictor variables as an additional predictor variable in the machine learning model like Gries (2020), and create a partial dependence plot for this additional predictor. The partial dependence plot will then capture the effect of the specified interaction. However, all these approaches do not allow to detect any unexpected and thus unspecified higher order interaction effects: Just like a univariate partial dependence plot cannot capture a 2-factor interaction, a partial dependence plot for a 2-factor interaction, a bivariate partial dependence plot or an ICE plot colored w.r.t. a second variable cannot capture a 3-factor interaction, etc. While we as humans tend to think only in interactions of order two or three at the most, the real pattern in the data may be more complex and contain interactions of higher order.

2.6 Detecting interactions

Gries (2020) suggests an approach based on random forest variable importances that has raised the hope to be able to better identify interactions and has already been picked up in the linguistics literature (Bernaisch and Funke 2024; Deshors 2021; Deshors and Gries 2020; Schmidt and Funke 2024): The first step of this approach is to explicitly add interactions between the predictor variables as additional predictor variables, termed interaction predictors in parts of the literature. This approach is perfectly legitimate and is also mentioned in Strobl et al. (2009b). It can increase the prediction accuracy in settings like the XOR case, where predictors have small or no main effects but may exhibit informative interactions, which are easier to detect when explicitly included. However, Gries (2020) suggestion goes beyond this first step and seems to imply that the resulting variable importance scores of the interaction predictors could help identify whether a predictor variable contributes to the prediction through a 2-factor interaction rather than through a main effect.

To illustrate this, Gries (2020) presents a case study with a toy data set. The example has been formulated such that an interaction predictor alone is able to perfectly predict the outcome. The interaction predictors are added and the finding is described this way (Gries 2020: 639): “Variable importance: every single forest chooses P2:P3 as the by far most important predictor, but it is worth noting that, because of the sampling, all other predictors’ variable importance scores are also not 0.” This indicates that this result is what the author considers as the correct behavior of the variable importance and to our understanding implies that the approach of including interaction predictors and computing their importance is promoted as a method for detecting truly relevant interaction effects.

Other papers seem to have interpreted Gries’ suggestion in the same way. For example, Deshors and Gries (2020: 225) write “[…] we follow Gries’s (forthcoming) recommendations: […] the first step of our statistical analysis consisted of manually creating a number of new predictors that essentially represent all two-way interactions […]. These were then added as predictors to a forest […]. […] Second, we computed the version of variable importance scores proposed in Janitza et al. (2013) […].” Deshors and Gries (2020) use a linear surrogate model in addition to the random forest analysis, but concentrate the interpretation on the four interaction predictors that have shown the highest AUC^[7] variable importance scores (besides the single variable variety; the variable importance scores for these four interaction predictors are reported on p. 226, the importance scores for all predictors in Appendix B of Deshors and Gries 2020), indicating that the authors consider these interaction effects to be most relevant. A similar approach was followed by Deshors (2021).

Schmidt and Funke (2024, p. 8 of online version) write “Following Gries (2020), we explicitly included interaction variables of the two sociolinguistic variables under investigation, namely GENDER and TIME, with each of the other variables. […]” and conclude (p. 9 of online version) “Figure 1 shows the variable importance scores measured in mean decrease in accuracy. While the main effects of GENDER and TIME rank among the least important variables, the plot highlights the high importance of many of the interaction effects for the model’s accuracy. Especially the interaction effects of GENDERxTRIGGER.LEMMA and TIMExTRIGGER.LEMMA turn out to be the most important variables in the random forest followed by the overall effect of TRIGGER.LEMMA and the interaction of TIMExNEWSPAPER.” The authors then go on to focus the interpretation on these effects. A similar approach was followed by Bernaisch and Funke (2024).

While it is correct to say that any (interaction) predictors with a high variable importance have highly contributed to predicting the response variable in a random forest, we would like to caution the readers against the idea that including interaction predictors provided a reliable way of detecting whether interaction effects are truly relevant. The evidence from our simulation study, presented below, shows that this is not the case. In particular, we will see that the importance of interaction predictors can be higher than that of individual predictors even when the corresponding interaction effect does not exist while the main effect of the individual predictor does. This means that the practice of concluding from a high variable importance ranking of interaction predictors that the corresponding interaction effect was truly relevant is not justified.

Unlike our simulation study, the case study of Gries (2020) only looks at the ability of the suggested approach to identify an interaction that is actually in the data (true positive finding).^[8] In order to validate a new methodological approach, however, it is important to look at different settings, also including ones that allow us to check whether the approach may erroneously identify interactions that are not in the data (false positive findings).