When including a linear interaction between two continuous predictors, should one generally also include quadratic predictors?

Question

Suppose I am fitting a linear model, and I have two continuous predictors x1 and x2. I think that they might interact, so I add the (linear) interaction term to my model (i.e., the pointwise product x1 * x2).

Now, as a rule of thumb, it seems to me that as soon as I do this, I should probably also include quadratic terms, x1^2 and x2^2. Otherwise, in case my predictors are correlated and have non-linear effects (and aren't everyone's predictors always correlated and having non-linear effects, at least to some extent?), then I might get a spuriously significant interaction. The correlation between predictors means that their linear interaction is confounded with the quadratic effect of either predictor alone, and so if I leave those quadratic predictors out of my model, then my model's best strategy to capture any single-variable quadratic trend is by using the two-variable interaction, making it look significant when it really isn't.

My question is:

• Is this reasoning valid?
• Is this a well known rule of thumb, e.g., is it stated as a recommendation in any standard books?

---

Another way of putting it: in R terms, suppose I find myself writing something like:

lm(y ~ (x1 + x2)^2, data)

Now the way R will interpret this is, it expands it to

lm(y ~ x1 + x2 + x1:x1 + x1:x2 + x2:x2, data)

but then it reduces the self-interactions x1:x1 and x2:x2 to just plain x1 and x2, so we end up with

lm(y ~ x1 + x2 + x1:x2, data)

So, if I find myself writing something like (x1 + x2)^2, should I instead make a habit of writing

lm(y ~ (x1 + x2)^2 + I(x1^2) + I(x2^2), data)

so as to include the quadratic "self-interactions"?

[Background: I'm wondering if I ought to alter patsy to interpret self-interactions for continuous variables like x:x as equivalent to I(x^2), instead of following R and interpreting them as equivalent to x. So I really am interested in the general case/"rule of thumb" rather than a specific instance.]

Answer 1

Is this reasoning valid?

Yes it is. See below.

Is this a well known rule of thumb, e.g., is it stated as a recommendation in any standard books?

I think it should be but I don't think it is, at least judging by the number of postgraduate students (and beyond) that haven't really considered it.

---

This question made me think about a section in the brilliant paper by W N Venables "Exegeses on Linear Models" (which is required reading for all my students) and I encourage anyone who has not read it, or hasn't read it recently, to do so. I will provide the full reference and link at the end of this answer. Most of what follows is taken from the paper, almost verbatim.

Let's start with a model that, on the face of it, is not very interesting:

$$ Y = f(x,Z)$$

where $x$ is a matrix of continuous explanatory variables and $Z$ is a random variable which we can think of as normally distributed around zero, but it does not have to be. If we take a first-order Taylor series approximation around $x_0$, then we have:

$$Y \approx f(x_0, 0) + \sum_{i = 1}^{p} f^{(i)}(x_0,0)(x_i-x_{i0}) + f^{(p+i)}(x_0,0)Z $$

or equivalently

$$Y \approx \beta_0 + \sum_{i = 1}^{p} \beta_{i}(x_i-x_{i0}) + \sigma Z $$

Note that it is common practice to subsume all the $x_{i0}$ into the intercept and then the model takes on a very familiar form. At this point we could naturally discuss whether a global intercept is a good idea, and whether centring the data could be of value.

If we continue with the Taylor series, the next approximation will be:

$$Y \approx \beta_0 + \sum_{i = 1}^{p} \beta_{i}(x_i-x_{i0}) + \sum_{i = 1}^{p}\sum_{j = 1}^{p} \beta_{ij}(x_i-x_{i0})(x_j-x_{j0}) + \left( \sigma + \sum_{i = 1}^{p} \gamma_i(x_i-x_{i0}) \right) Z + \sigma Z^2 $$

and so now we find:

• nonlinearty in the main effect (quadratic terms)
• a linear x linear interaction (cross product of two linear terms)
• heteroskedasticity (the $(x_i-x_{i0}) Z$ terms)
• skewness (the term in $Z^2$)

So this brings us back to the question about whether we should include quadratic terms when we include an interaction, and this approach tells us that we should. In general I think it is always a good idea to consider quadratic terms when fitting an interaction. Perhaps a good question to ask is why should we not do so ?

I always encourage students and colleagues to step back and look at what we are doing from a wider viewpoint. Statistical models are models, an abstraction of reality, which as George Box famously said, are all wrong, but some are useful. It is our job to make them as useful as possible, whether that be for prediction or inference.

It might very well be the case that in a particular context (eg a small region of the domain of $x$) that the nonlinear (and/or) interaction terms will not be needed, but at the very least it is a good idea to think about this, and the same goes for heteroskedasticity and skewness.

Source for the above:

Venables, W.N., 1998, October. Exegeses on linear models. In S-Plus User’s Conference, Washington DC.
http://www.stats.ox.ac.uk/pub/MASS3/Exegeses.pdf

Answer 2

I have an argument for, and an argument against. My recommendation is not to force the quadratic terms when adding interaction. If they fit your problem, then go ahead and add them, but not as a general rule of thumb.

Yes

Include quadratic terms, because 2nd order Taylor expansion would suggest so:

$$y=f(x,z)+e\approx f_x\Delta x+f_z\Delta z+2f_{xz}\Delta x\Delta z+f_{xx}\Delta x^2+f_{zz}\Delta z^2+e$$

From a point of view of mathematical elegance it would make a sense to throw in the quadratic terms whenever the interaction is added, if you like Taylor expansion interpretation of the regression coefficients.

No

This would make your approximation quadratic, which usually is not what you want:

$$y=\beta_0+\beta_xx+\beta_zz+\beta_{xz}xz+\beta_{xx}x^2+\beta_{zz}z^2+e$$

So, e.g. when $x\to-\infty$ then $y\sim |x|^2$, i.e. increasing at quadratic speed. This is rarely a desired behavior in linear regression because we usually expect monotonic effect of regressors. On the other hand, without quadratic terms the dependence on $x$ remains linear with an interaction $(\beta_x+\beta_{xz}z)x$, although the slope varies with $z$, which is a desired behavior.

This is the reason, for instance, why in LOESS type of local polynomial fits the odd order polynomials are used with p=1 or 3.