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

```r
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. 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.$$

```r
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:

```r
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_{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.

---

Thanks for any helps!