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. 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:
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
withcyl=6
is different from that withcyl=4
(baseline) after controllingdrat
.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
withcyl=6
is the same as that withcyl=4
(baseline) after controllingdrat
.How to perform the multiple comparisons on each pair of
cyl
after controllingdrat
? In R,pairwise.t.test()
doesn't seem to handle this.
Thanks for any helps!