how to interpret classes dependence that are not the reference class in a linear model

Question

If we run the three following codes:

n <- 250
df <- rbind(
  data.frame(cat=1, x=rnorm(n, 0), y=rnorm(n, 1)),
  data.frame(cat=2, x=rnorm(n, 0), y=rnorm(n, 0)),
  data.frame(cat=3, x=rnorm(n, 0), y=rnorm(n, 0))
)
df$cat <- as.factor(df$cat)
df$cat <- relevel(df$cat, ref = "1")

summary(lm(y ~ cat, df))

It will say here that intercept, cat2 and cat3 are statistically significant in trying to predict y. But doesn't indicate that cat=2 and cat=3 are actually redundant.

n <- 250
df <- rbind(
  data.frame(cat=1, x=rnorm(n, 0), y=rnorm(n, 1)),
  data.frame(cat=2, x=rnorm(n, 0), y=rnorm(n, 0)),
  data.frame(cat=3, x=rnorm(n, 0), y=rnorm(n, 0))
)
df$cat <- as.factor(df$cat)
df$cat <- relevel(df$cat, ref = "2")

summary(lm(y ~ cat, df))

Here only cat1 is significant.

n <- 250
df <- rbind(
  data.frame(cat=1, x=rnorm(n, 0), y=rnorm(n, 1)),
  data.frame(cat=2, x=rnorm(n, 0), y=rnorm(n, 0)),
  data.frame(cat=3, x=rnorm(n, 0), y=rnorm(n, 0))
)
df$cat <- as.factor(df$cat)
df$cat <- relevel(df$cat, ref = "3")

summary(lm(y ~ cat, df))

Same here (which makes sense).

Do we have to run several times the linear model to see that cat2 and cat3 don't need to be both used to predict y? What if only the intercept is significant, what does that mean? I'm not able to reproduce a case where that happens and it happens in my dataset. Does that mean we don't need any of the 3 categorical variables to predict y? Why wouldn't be all p values not significant in that case?

Answer 1

This is a frequent confusion arising from the way that coefficient estimates and p-values are typically displayed for multi-level categorical predictors. With this default treatment coding, the intercept is the outcome at the reference level and the coefficients are the differences from that outcome associated with the other levels.

The choice of reference level matters in terms of display, but not for the fundamental model. If you change the reference level, both the intercept and the other coefficient estimates will change.

In this particular case you can get a clearer display if you omit the intercept. (I set the seed to 2552 before continuing with your code to generate df, and kept the reference level at "1".)

lm2 <- lm(y ~ -1 + cat, df)

Instead of focusing on the individual coefficients and p-values, examine the overall significance of the categorical predictor. With these data, standard anova() is a good choice (particularly as this is just a balanced ANOVA in lm form). The predictor is highly significant, as intended.

anova(lm2)
# Analysis of Variance Table
# 
# Response: y
#            Df Sum Sq Mean Sq F value    Pr(>F)    
# cat         3 269.47  89.824  86.183 < 2.2e-16
# Residuals 747 778.56   1.042                      

With the intercept omitted, each of the levels is shown with its own mean value and the significance of its difference from 0.

summary(lm2)

## some output lines omitted

# Coefficients:
#      Estimate Std. Error t value Pr(>|t|)    
# cat1  1.03630    0.06457  16.050   <2e-16
# cat2 -0.05519    0.06457  -0.855    0.393    
# cat3  0.03040    0.06457   0.471    0.638    

That trick, however, only works with a single categorical predictor.

In general with treatment coding, the intercept represents a situation when all categorical predictors are at their reference levels and continuous predictors are at 0. So the "significance" of the intercept is just whether the outcome is significantly different from 0 for that particular choice of reference levels and centering of continuous predictors. In many circumstances that "significance" isn't very useful on its own.

Similarly the "significance" of a coefficient for a non-reference level is whether its difference from the outcome at the reference level can be distinguished from 0. So that also depends on the choice of reference level.

Nevertheless, outcome predictions for any level of a categorical predictor will be the same regardless of reference-level or predictor-centering choices. It's often wise just to ignore the coefficient p-values in this situation and look directly at predictions or differences between predictions.

The emmeans package provides a general way to display more useful model summaries of that type. It works with the original model to display particular estimates or comparisons of interest.

library(emmeans)
emmeans(lm2,pairwise~cat)
# $emmeans
#  cat  emmean     SE  df lower.CL upper.CL
#  1    1.0363 0.0646 747   0.9095   1.1631
#  2   -0.0552 0.0646 747  -0.1819   0.0716
#  3    0.0304 0.0646 747  -0.0964   0.1572
# 
# Confidence level used: 0.95 
# 
# $contrasts
#  contrast estimate     SE  df t.ratio p.value
#  1 - 2      1.0915 0.0913 747  11.953  <.0001
#  1 - 3      1.0059 0.0913 747  11.016  <.0001
#  2 - 3     -0.0856 0.0913 747  -0.937  0.6168
# 
# P value adjustment: tukey method for comparing a family of 3 estimates 
# 

The $emmeans are the outcome estimates for each level, the $contrasts are the estimates of differences between each pair of predictor levels, with a p-value correction for multiple comparisons. You would get the same results if you used modeled the same data while including the intercept, or with any other choice of reference level.