How to systematically choose which interactions to include in a multiple regression model?

Question

In an answer to this post a user suggests, based on Chapter 3 of the book "The Elements of Statistical Learning" by Hastie et. al, the following means of choosing which interaction effects to include in a model:

1. Trying out all possible subsets of variables and pick the one that gives a regression with the smallest Bayesian information criterion (BIC) value
2. Forward or backward stepwise selection

In the comments associated with that answer, both of these approaches are described as being bad.

So, if we shouldn't use method 1) or 2) above, how exactly do we decide what variables/interactions to use in the model? I have seen 'domain knowledge' suggested in a few places, but this seems like a bit of a cop out. Domain matter knowledge is not going to help in the very common situation in which we have no pre-existing knowledge of whether a particular interaction effect is present in nature and we are relying on the information in the data itself.

For the sake of example, suppose we have the predictors - age, gender, height, weight, experience, IQ - and the response variable salary. How do we decide what interaction effects to include/not include?

This example is probably the simplest possible scenario, as we understand all of these variables very well, and even still it is not clear how to decide which interactions to include or exclude. In other situations, we will be dealing with predictor variables for which we have no pre-existing intuition on whether interactions between them could affect the response variable.

So I am looking for a systematic method of choosing which interactions to include in a multiple regression model. How does an experienced statistician choose which interactions to include in the case when domain knowledge is not available or of no use?

Answer 1

I think a lot depends on what the purpose of the model is. Inference or prediction ?

If it is inference then you really need to incorporate some domain knowledge into the process, otherwise you risk identifying completely spurious associations, where an interaction may appear to be meaningful but in reality is either an artifact of the sample, or is masking some other issues such as non-linearity in one of more of the variables.

However, if the purpose is prediction then there are a various approaches you can take. One approach would be to fit all possible models and use a train / validate / test approach to find the model that gives the best predictions.

---

Edit : A simple simulation can show what can go wrong with inference without domain knowledge:

set.seed(50)
N <- 50

X1 <- runif(N, 1, 15)
X2 <- rnorm(N)
  
Y <-  X1 + X2^2 + rnorm(N)

So, here we posit an actual data generation process of $Y = X_1 + {X_2}^2$

If we had some domain / expert knowledge that suggested some nonlinearities could be involved, we might fit the model:

> lm(Y ~ X1 + I(X1^2) + X2 + I(X2^2) ) %>% summary()

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.89041    0.65047  -1.369    0.178    
X1           1.21915    0.19631   6.210 1.52e-07 ***
I(X1^2)     -0.01462    0.01304  -1.122    0.268    
X2          -0.19150    0.15530  -1.233    0.224    
I(X2^2)      1.07849    0.08945  12.058 1.08e-15 ***

which provides inferences consistent with the "true" data generating procress.

On the other hand, if we had no knowledge and instead thought about a model with just first order terms and the interaction we would obtain:

> lm(Y ~ X1*X2) %>% summary()

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.01396    0.58267  -0.024    0.981    
X1           1.09098    0.07064  15.443  < 2e-16 ***
X2          -3.39998    0.54363  -6.254 1.20e-07 ***
X1:X2        0.35850    0.06726   5.330 2.88e-06 ***

which is clearly spurious.

---

Further edit : However, when we look at predictive accuracy using root mean squared error we find that the interaction model performs slightly better:

> lm(Y ~ X1*X2) %>% predict() %>% `^`(2) %>% sum() %>% sqrt()
[1] 64.23458
> lm(Y ~ X1 + I(X1^2) + X2 + I(X2^2) ) %>% predict() %>% `^`(2) %>% sum() %>% sqrt()
[1] 64.87996

which underlines my central point that a lot depends on the purpose of the model.