What test can I use to compare intercepts from two or more regression models when slopes might differ?

Question

I wish to test whether intercepts in linear regression models differ between two or more groups, when group-specific slopes might themselves differ (i.e., an interaction term may be present). Specifically, I want to compare intercepts between all pairwise combinations of groups.

The question entitled What test can I use to compare slopes from two or more regression models? shows how to test whether the slopes differ between all pairwise combinations of groups. However, I have not been able to find an equivalent way to test whether the intercepts differ between pairwise combinations of groups.

ANOVA (e.g., anova(lm(Sepal.Length ~ Petal.Width*Species, data = iris))) will tell you whether the intercepts of each group differ relative to the baseline group (i.e., the first level in the grouping factor). However, is there a convenient way to perform all pairwise comparisons between groups?

Answer 1

I will answer the technical question, then try to talk you out of doing this.

The intercept is the predicted value when the abscissa is equal to zero. Hence, the intercepts in the example are obtained via:

> mod = lm(Sepal.Length ~ Petal.Width*Species, data = iris)

> library("emmeans")
> (emm = emmeans(mod, "Species", at = list(Petal.Width = 0)))
NOTE: Results may be misleading due to involvement in interactions
 Species    emmean    SE  df lower.CL upper.CL
 setosa       4.78 0.173 144     4.43     5.12
 versicolor   4.04 0.464 144     3.13     4.96
 virginica    5.27 0.509 144     4.26     6.28

Confidence level used: 0.95 

... and the comparisons thereof can be tested this way:

> pairs(emm)
 contrast               estimate    SE  df t.ratio p.value
 setosa - versicolor       0.733 0.495 144  1.480  0.3037 
 setosa - virginica       -0.492 0.538 144 -0.915  0.6316 
 versicolor - virginica   -1.225 0.689 144 -1.779  0.1804 

P value adjustment: tukey method for comparing a family of 3 estimates 

That said, it is an unusual instance that the intercept is an interesting or meaningful quantity to want to do inferences about. In many datasets, the intercept is a severe extrapolation because a predictor value of zero is distant from its observed values. Models are only approximations to the truth, and it is highly questionable that one should believe that the straight line you have fitted actually represents the trend at a distant point.

Thus, I urge you to re-think what you are trying to do and decide on what meaningful question you are really trying to answer.

Answer 2

In principle, once you have the linear regression object generated by lm(), you can test for significant differences between any desired linear combinations of predictor values that you wish by applying the formula for the variance of a sum to the covariance matrix for the linear regression, the matrix provided provided by vcov(lm()). You should correct, of course, for the multiple comparisons. The emmeans package provides one convenient way to perform such tests.

That said, you need to be very careful about comparing intercepts when interaction terms are involved, particularly when one of the predictors is continuous with values far from 0. In your example, the intercept represents the Sepal.Length for each Species when the Petal.Width is exactly 0. Is that really a comparison you care about?

In the iris dataset, the average Petal.Width values are 0.246, 1.326, and 2.026 for the Species setosa, versicolor, and virginica respectively. With such a range of Petal.Width values, what would be the point of examining the intercepts for Sepal.Lengthin the fictitious case when the Petal.Width is 0?

I understand that you are using this dataset only as an example of a more general question, but the same cautions apply broadly to evaluations of intercepts when interactions are involved.