# Perform the multiple comparisons after controlling other covariates

To describe this issue, we consider the following 2 models:

m1 <- lm(disp ~ drat, mtcars)
m2 <- lm(disp ~ drat + factor(cyl), mtcars)


drat is a continuous variable. cyl is a categorical variable with levels '4', '6', '8', and I treat '4' as the baseline, so there will be 2 dummy variables $$D_1$$ for cyl=6 & $$D_2$$ for cyl=8 respectively. The regression model is:

$$disp = \beta_0 + \beta_1 drat + \beta_{cyl=6} D_1 + \beta_{cyl=8} D_2 + \epsilon$$

I want to know if cyl is correlated to disp after controlling drat. I perform an F-test where the null hypothesis is

$$H_0: \beta_{cyl=6} = \beta_{cyl=8} = 0.$$

anova(m1, m2)

# Model 1: disp ~ drat
# Model 2: disp ~ drat + factor(cyl)
#   Res.Df    RSS Df Sum of Sq      F    Pr(>F)
# 1     30 235995
# 2     28  70246  2    165750 33.034 4.286e-08 ***


The p-value ($$\approx$$ 0) implies that we can reject the $$H_0$$ statement. Now I wonder which pairs of cyl have significant difference in disp. The following is the model summary table:

summary(m2)

# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)
# (Intercept)    266.28      97.32   2.736   0.0107 *
# drat           -39.58      23.62  -1.676   0.1048
# factor(cyl)6    58.97      26.79   2.201   0.0361 *
# factor(cyl)8   214.65      28.33   7.578 2.96e-08 ***


I have 2 main questions:

1. From the table above, the p-value of $$\beta_\rm{cyl=6}$$ is $$0.036 < 0.05$$. Could I interpret this value directly and conclude:

The mean disp with cyl=6 is different from that with cyl=4 (baseline) after controlling drat.

Or otherwise, should I adjust the p-value to avoid inflation of the type I error? E.g. multiply it by 3 (using Bonferroni correction) to get p-value $$= 0.036 \times 3 = 0.108 > 0.05$$, and hence conclude:

The mean disp with cyl=6 is the same as that with cyl=4 (baseline) after controlling drat.

2. How to perform the multiple comparisons on each pair of cyl after controlling drat? In R, pairwise.t.test() doesn't seem to handle this.

Regarding both pairwise comparisons and adjustment for multiple comparisons, doesn't emmeans solution work?

library(emmeans)
em1<-emmeans(m2, ~cyl) #this gives you estimated marginal means of disp at different levels of cyl from the model m2
em1

From the table above, the p-value of $$βcyl=6$$ is 0.036<0.05 Could I interpret this value directly and conclude: "The mean disp with cyl=6 is different from that with cyl=4 (baseline) after controlling drat"