Jeff Gill writes:
For some reason the misinterpretations about interactions in regression models just won’t go away. I teach the point that mathematically and statistically one doesn’t have to include the main effects along with the multiplicative component, but if you leave them out it should be because you have a strong theory supporting this decision (i.e. GDP = Price * Quantity, in rough terms). Yet I got this email from a grad student yesterday:
As I was reading the book, “Introduction to Statistical Learning,” I came across the following passage. This book is used in some of our machine learning courses, so perhaps this is where the idea of leaving the main effects in the model originates. Maybe you can send these academics a heartfelt note of disagreement.
“The hierarchical principle states that if we include an interaction in a model, we should also include the main effects, even if the p-values associated with their coefficients are not significant. In other words, if the interaction between X1 and X2 seems important, then we should include both X1 and X2 in the model even if their coefficient estimates have large p-values. The rationale for this principle is that if X1 × X2 is related to the response, then whether or not the coefficients of X1 or X2 are exactly zero is of little interest. Also X1 × X2 is typically correlated with X1 and X2, and so leaving them out tends to alter the meaning of the interaction.”
(Bousquet, O., Boucheron, S. and Lugosi, G., 2004. Introduction to statistical learning theory. Advanced Lectures on Machine Learning: ML Summer Schools 2003, Canberra, Australia, February 2-14, 2003, Tübingen, Germany, August 4-16, 2003, Revised Lectures, pp.169-207.)
There are actually two errors here. It turns out that the most cited article in the history of the journal Political Analysis was about interpreting interactions in regression models, and there are seemingly many other articles across various disciplines. I still routinely hear the “rule of thumb” in the quote above.
To put it another way, suppose you start with the model with all the main effects and interactions, and then you consider the model including the interactions but excluding one or more main effects. You can think of this smaller model in two ways:
1. You could consider it as the full model with certain coefficients set to zero, which in a Bayesian sense could be considered as very strong priors on these main effects, or in a frequentist sense could be considered as a way to lower variance and get more stable inferences by not trying to estimate certain parameters.
2. You could consider it as a different model of the world. This relates to Jeff’s reference to having a strong theory. A familiar example is a model of the form, y = a + b*t + error, with a randomly assigned treatment z that occurs right after time 0. A natural model is then, y = a + b*t + c*z*t + error. You’d not want to fit the model, y = a + b*t + c*z * d*z*t + error—except maybe as some sort of diagnostic test—because, by design, the treatment cannot effect y at time 0.
I have three problems with the above-quoted passage. The first is the “even if the p-values” bit. There’s no good reason, theoretically or practically, that p-values should determine what is in your model. So it seems weird to refer to them in this context. My second problem is where they say, “whether or not the coefficients of X1 or X2 are exactly zero is of little interest.” In all my decades of experience, whether or not certain coefficients are exactly zero is never of interest! I think the problem here is that they’re trying to turn an estimation problem (fitting a model with interactions) into a hypothesis testing problem, and I think this happened because they’re working within an old-fashioned-but-still-dominant framework in theoretical statistics in which null hypothesis significance testing is fundamental. Finally, calling it a “hierarchical principle” seems to be going too far. “Hierarchical heuristic,” perhaps?
That all said, usually I agree with the advice that, if you include an interaction in your model, you should include the corresponding main effects too. Hmmm . . . let’s see what we say in Regression and Other Stories . . . section 10.3 is called Interactions, and here’s what we’ve got . . .
We introduce the concept of interactions in the context of a linear model with a continuous predictor and a subgroup indicator:
Figure 10.3 suggests that the slopes differ substantially. A remedy for this is to include an interaction . . . that is, a new predictor defined as the product of these two variables. . . . Care must be taken in interpreting the coefficients in this model. We derive meaning from the fitted model by examining average or predicted test scores within and across specific subgroups. Some coefficients are interpretable only for certain subgroups. . . .
An equivalent way to understand the model is to look at the separate regression lines for [the two subgroups] . . .
Interactions can be important, and the first place we typically look for them is with predictors that have large coefficients when not interacted. For a familiar example, smoking is strongly associated with cancer. In epidemiological studies of other carcinogens, it is crucial to adjust for smoking both as an uninteracted predictor and as an interaction, because the strength of association between other risk factors and cancer can depend on whether the individual is a smoker. . . . Including interactions is a way to allow a model to be fit differently to different subsets of data. . . . Models with interactions can often be more easily interpreted if we preprocess the data by centering each input variable about its mean or some other convenient reference point.
We never actually get around to giving the advice that, if you include the interaction, you should usually be including the main effects, unless you have a good theoretical reason not to. I guess we don’t say that because we present interactions as flowing from the main effects, so it’s kind of implied that the main effects are already there. And we don’t have much in Regression and Other Stories about theoretically-motivated models. I guess that’s a weakness of our book!
“In all my decades of experience, whether or not certain coefficients are exactly zero is never of interest!”
Does this sentence mean that variable selection is of no interest at all, and people working on it should find a nice retirement home?
Anon:
Variable selection (that is, setting some coefficients to be exactly zero) can be useful for various reasons, including:
– It’s a simple form of regularization.
– It can reduce costs in future data collection.
Variable selection can be fine as a means to an end. Problems can arise if it’s taken too seriously, for example as an attempt to discover a purported parsimonious true model.
Include the main effect if it improves the out of sample predictive skill enough to outweigh the costs of increased memory and computation time.
Why is it not so basic and simple?
I think this is another case someone wracking their brain trying to interpret arbitrary model coefficients.
Arbitrary numbers have no meaning, and none of them really exactly equal zero because everything that happens is collectively caused by everything in its own past light cone. But it is fine to leave out the negligible ones (ie, set them to zero).
If you are using variable selection to find a useful prediction model (the typical scenario), then it is the simplified model that is of actual interest, not whether or not particular coefficients are zero. Setting some coefficients to zero via some procedure (whether statistical or first principles) is not the same as being interested in testing whether they are different from zero as the inferential goal of one’s research.
(this was supposed to be a reply to Anon above)
I think there’s a big misunderstanding about what “interaction effects” are. And this is because there’s a big misunderstanding about what regression is.
Suppose you have an experimental outcome Y that is influenced by the status of two measurements, A and B.
Let’s assume it’s not a dynamic process but rather something where it’s appropriate to write
Y = f(A,B) + err
the problem is we don’t know f
One way to address this is to use a universal approximation method to approximate f. Suppose we hold B constant, then we could write f as just a function of A, and if it’s a continuous function on a finite domain, which the vast majority of real world situations are, we could use a polynomial in A (note, for various reasons you should usually use chebyshev polynomials but here I use the A^n basis for educational purposes)
Y = a +bA +cA^2 +dA^3…
In the same sense B is symmetric we could use a polynomial in B holding A constant.
If both A and B can vary, what do we do? Simple. We write the polynomial in one variable with coefficients which are each functions of the other variable.
Y = P1(A)+P2(A)B+P3(A)B^2…
Where each Pn is itself a polynomial. Suppose we go up to 2nd order in both directions, then there are 9 unknowns.
The shape of the B polynomial is now a function of A, so for any cross section for fixed A, the curve is just some shape in B but if we change the A the B curve is a different shape.
If we use a polynomial of degree M and K then there are (M+1)(K+1) coefficients. If we were to use the Chebyshev basis these coefficients could be related to the values of the surface at special chebyshev points. We can think of an interpolator surface determined by those (M+1)(K+1) “grid points”
So far nothing we have done is special to the problem at hand. It’s just a universal way to approximate continuous functions on closed bounded domains.
You can make this be about your particular science problem by imposing prior information about the problem. Basically prior information will constrain the shape of this universal function. For example you could say that the P1(A),P2(A) etc are all the same… Or you could set some boundary conditions that on the edge of the domain the function goes to zero, or to a value that lies on a plane, or you could say that interior chebyshev points are all higher than the boundary points, or that the function is nearly flat around the edges, or whatever. Each kind of prior information specifies something about the behavior of the function which itself specifies something about the coefficients. So you can place prior info on the coeffs or on the behavior, they are dual to each other.
The stuff that is usually done in stats is to take this beautiful universal approximation method and botch it. People use low order polynomials, then they impose their prior info by either saying a coefficient could be literally any number (a flat nonstandard prior on the real line) or has to be exactly 0 (a delta function prior on 0).
Both are wrong headed. The useful way to specify prior information is via probability placed on the function BEHAVIOR.
I got tired of writing stuff on my phone, so I posted that. now I’m on my computer so I wanted to add some info about what kinds of behavior you could describe.
For example, you could say that the function isn’t very “wiggly”. Wigglyness could mean changes in slope. So let’s calculate the average of the Laplacian… Lavg = integrate(d^2Y/dA^2 + d^2Y/dB^2, dA,dB) over the domain. This quantity is a scalar that is a function of **all the coefficients** so we can utilize it as part of our prior by something like
Lavg ~ normal(0.0,scale)
Think of this as a component of the prior multiplying the pdf of the coefficients. To normalize it would be a pain in the ass but with MCMC we don’t need to worry about that.
If you make scale a small number, then the universal function will be relatively smooth, whereas if it’s a bigger number, the functions will be allowed to be more wiggly.
There are all sorts of variations on this. Do you think that once B gets above a certain region, the function is very nearly constant? then integrate((TheFunction(A,B)-Constant)^2,dA,dB) in the region you were thinking of, and constrain it to be near 0.
Etc.
We need to stop thinking about setting coefficients to particular values or “excluding interactions” or etc and just think about “what does the function do?”
I have a suspicion that how statistical software commonly works explains some of this. For example in R, if you fit a model with lm and then call summary, it will happily output a list of p-values for each coefficient, where the meanings of coefficients is super unclear unless you know exactly what you are doing. A user might reasonably want to know what to do with this important looking output (p-values!). For example, values of coefficients seemingly about the main effects and interactions of categorical variables have a non-obvious interpretation. This might be behind the idea “we should also include the main effects, even if the p-values associated with their coefficients are not significant”.
I think the step that is often missed is using something like Estimated Marginal Means to interrogate a model. Get estimates of quantities actually of interest, with confidence intervals. This is the right sort of output to be looking at.
I see many people using statistical software with quite a fuzzy idea of what is going on. Yes they should be taught better, but better software could help create a pit of success too. (And then maybe we can talk sensibly about what they should include in their model.)
The confusion is real. Stuff like figuring out which “contrasts” are being used and how those contrasts are encoded and what they mean used to drive me absolutely crazy. These days I just use Bayesian model fitting software and so I know exactly what mathematical calculation is being done. I’ll do a quick least squares linear model here and there but for the most part I really hate doing anything even slightly complicated via stuff like “lm” in R or Julia
99+% of the time a model is used because of whatever data is available, ie the model is arbitrary. Since the coefficients are conditional on the model, that makes them also arbitrary.
So far no one has been able to explain why they hope to discover “meaning” in these arbitrary numbers. My guess is the supposed people who find this meaning because they “know exactly what they are doing” are either mythical or committing some kind of fallacy.
Ah, I am mostly involved in designed experiments. My collaborators designed the experiment to contain the experimental samples with the treatments they are interested in, with a fully-crossed design. I know the causal graph because I can talk to the people who did the causing. I know what causal effects they are trying to estimate, and those are the meaningful numbers I want.
I would combine “strong theory supporting this decision (i.e. GDP = Price * Quantity, in rough terms)” and treatment given at the start of experiment example of Prof. Gelman into “unless 0 has some special significance in your model either for treatments or for the outcome”. Because if 0 is just another number we can shift everything to Z-Z0= a(x-x0)(y-y0) which is 4 parameters. But if you can think of why any of them should be exactly 0, you probably should incorporate it in some way.
Interesting discussion. I agree with the comments RE how interaction effects are often misunderstood. Years ago I wrote an overview article in the Journal of Statistics Education to help clarify many of the concepts, “Key Results of Interaction Models with Centering” (2011),
https://www.tandfonline.com/doi/abs/10.1080/10691898.2011.11889620
Jonas Haslbeck later wrote a summary blog post about the article on R-bloggers (2018):
https://www.r-bloggers.com/2018/05/regression-with-interaction-terms-how-centering-predictors-influences-main-effects/
Also see
Schielzeth, Holger. “Simple Means to Improve the Interpretability of Regression Coefficients: Interpretation of Regression Coefficients.” Methods in Ecology and Evolution 1, no. 2 (February 10, 2010): 103–13. https://doi.org/10.1111/j.2041-210X.2010.00012.x
which I recommend to people approximately once a week …
PS the “ff” and “ffi” glyphs seem to have disappeared from your quotation from “Regression and Other Stories” in the process of cutting and pasting text …
Fixed; thanks.
The “hierarchy principle” is defined in the classical 1975 text “Discrete Multivariate Analysis” by Bishop, Fienberg, and Holland. I don’t know if that’s where the “hierarchical principle” above originates, but Bishop et al. seem to be claiming to be the creators of this principle. They define the hierarchy principle for log linear modeling of entries of contingency tables: “The family of hierarchical models is defined as the family such that if any u-term is set equal to zero, all its higher-order relatives must also be set equal to zero. Conversely, if any u-term is not zero, its lower-order relatives must be present in the log-linear model. Thus if u_{12} = 0, we must have u_{123} = 0”.
They go on, “In larger tables with nonhierarchical structure a possible strategy is to partition and look at smaller sections of data. If u_{123} is not equal to 0 but u_{12} = 0 when we partition according to the categories of variable 3, then we have some tables with interaction between variables 1 and 2 in one direction and others with interaction in the opposite direction. We can either consider the structure of each set of tables separately or rearrange the cells to form new compound variables.
I have found that when researchers exclude the main effect, it is usually because the result runs counter to what everyone expects. Following from the example in the post, if someone is researching the association between a risk factor and cancer and the interaction of smoking * risk factor is positive and significant, but the main effect of smoking is negative and significant, then the researcher excludes the main effect and hopes nobody notices. In some of these studies, the total effect of smoking is often negative, even though the interaction is as predicted.
Would it be okay to exclude the main effects if they are collinear with the fixed effects and thus fully absorbed? The software is dropping them. I found that writing the model with or without the main effects has identical outcomes.