# Controlling covariates in linear regression in R

I have made these two models:

(model1 <- summary(lm(mpg ~ drat + wt + cyl, mtcars)))

Call:
lm(formula = mpg ~ drat + wt + cyl, data = mtcars)

Residuals:
Min      1Q  Median      3Q     Max
-4.2944 -1.5576 -0.4667  1.5678  6.1014

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  39.7677     6.8729   5.786 3.26e-06 ***
drat         -0.0162     1.3231  -0.012 0.990317
wt           -3.1947     0.8293  -3.852 0.000624 ***
cyl          -1.5096     0.4464  -3.382 0.002142 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.613 on 28 degrees of freedom
Multiple R-squared:  0.8302,    Adjusted R-squared:  0.812
F-statistic: 45.64 on 3 and 28 DF,  p-value: 6.569e-11

(model2 <- summary(lm(mpg ~  wt + cyl + drat, mtcars)))

Call:
lm(formula = mpg ~ wt + cyl + drat, data = mtcars)

Residuals:
Min      1Q  Median      3Q     Max
-4.2944 -1.5576 -0.4667  1.5678  6.1014

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  39.7677     6.8729   5.786 3.26e-06 ***
wt           -3.1947     0.8293  -3.852 0.000624 ***
cyl          -1.5096     0.4464  -3.382 0.002142 **
drat         -0.0162     1.3231  -0.012 0.990317
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.613 on 28 degrees of freedom
Multiple R-squared:  0.8302,    Adjusted R-squared:  0.812
F-statistic: 45.64 on 3 and 28 DF,  p-value: 6.569e-11


My understanding is that R uses "sequential" partitioning for the variance in mpg. So in model1, drat should be unadjusted, wt should be adjusted for drat and cyl should be adjusted for drat and wt. In model2, wt should be unadjusted, cyl should be adjusted for wt and drat should be adjusted for wt and cyl.

However, the coefficients in each model appear to be exactly the same, suggesting coefficients are not being adjusted at all. Are the coefficients not being adjusted at all?

-

The question, as phrased, is slightly ambiguous. It states that "the coefficients in each model appear to be exactly the same". There are two ways that statement could be interpreted, with respect to: (1) the Estimates of the coefficients, or (2) the tests of the coefficients.

1. Regarding the Estimates of the coefficients, they are being adjusted for the other variables in the model, but you don't see any difference because you have the same variables in both model1 and model2; the order they are listed in doesn't matter. The parameter estimates will only differ if the variables are correlated and the sets of variables that are included in the models differ. Consider:

model1 <- lm(mpg ~ drat + wt + cyl,  mtcars)
model2 <- lm(mpg ~  wt + cyl + drat, mtcars)
model3 <- lm(mpg ~  wt       + drat, mtcars)
cor(mtcars$wt, mtcars$cyl)
# [1] 0.7824958
summary(model2)
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)
#   (Intercept)  39.7677     6.8729   5.786 3.26e-06 ***
#   wt           -3.1947     0.8293  -3.852 0.000624 ***
#   cyl          -1.5096     0.4464  -3.382 0.002142 **
#   drat         -0.0162     1.3231  -0.012 0.990317
summary(model3)
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)
#   (Intercept)   30.290      7.318   4.139 0.000274 ***
#   wt            -4.783      0.797  -6.001 1.59e-06 ***
#   drat           1.442      1.459   0.989 0.330854


Notice that the Estimate for wt is -3.1974 in model2 and -4.783 in model3. To learn more about how the parameter estimates for variables can change depending on whether a correlated variable is included or not, it may help you to read by answer here.

2. Regarding the tests of the coefficients, it depends on which function you use to get your output. You used summary(). What is reported then are the $t$-tests associated with the parameter estimates. Those are not computed by partitioning sums of squares. However, they are equivalent to $F$-tests using type III SS. The order variables are listed in doesn't matter for $t$-tests or $F$-tests that are based on type III SS. You can also use anova() to get significance tests of your parameter estimates. That is where R uses type I SS. And because type I SS are sequential, the order the variables are listed in does matter (although again, only if the variables are correlated). Consider:

summary(model1)
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)
#   (Intercept)  39.7677     6.8729   5.786 3.26e-06 ***
#   drat         -0.0162     1.3231  -0.012 0.990317
#   wt           -3.1947     0.8293  -3.852 0.000624 ***
#   cyl          -1.5096     0.4464  -3.382 0.002142 **
summary(model2)
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)
#   (Intercept)  39.7677     6.8729   5.786 3.26e-06 ***
#   wt           -3.1947     0.8293  -3.852 0.000624 ***
#   cyl          -1.5096     0.4464  -3.382 0.002142 **
#   drat         -0.0162     1.3231  -0.012 0.990317
anova(model1)
# Analysis of Variance Table
#             Df Sum Sq Mean Sq F value    Pr(>F)
#   drat       1 522.48  522.48  76.525 1.691e-09 ***
#   wt         1 334.33  334.33  48.967 1.308e-07 ***
#   cyl        1  78.07   78.07  11.435  0.002142 **
#   Residuals 28 191.17    6.83
anova(model2)
# Analysis of Variance Table
#             Df Sum Sq Mean Sq  F value    Pr(>F)
#   wt         1 847.73  847.73 124.1627 8.382e-12 ***
#   cyl        1  87.15   87.15  12.7645  0.001304 **
#   drat       1   0.00    0.00   0.0001  0.990317
# Residuals 28 191.17    6.83


Notice that the $p$-value for wt in model1 is 0.000624 in both summary() outputs, but is 1.308e-07 in anova(model1) and it 8.382e-12 in anova(model2). To learn more about sums of squares in general, it may help to read my answer here. Lastly, note that you can get an ANOVA table in R that uses other types of SS, such as II and III, by using Anova() in the car package.

-
Hmmm…regarding your first point R in Action (2011, Kabacoff) states something different. He says order matters and in general, models should be ordered as covariates, main effects, interactions. In the model y~a+b+a:b, a will be unadjusted, b will adjust for a and a:b will adjust for a and b. He states adjustment should occur if (a) unequal number of observations in any predictor OR (b) covariates are present –  luciano Mar 6 '14 at 16:24
I'm not sure what that is supposed to mean, @luciano. If that is what he says, he's wrong, although there may be some ambiguity in the way he put something & he didn't quite mean that. You can see in the examples what changes & what doesn't, & under what circumstances. All those examples were simply run in R & cut&pasted here. You can run them yourself & get the same results. –  gung Mar 6 '14 at 16:28
Having re-read Kabacoff, I think he's just refering to sum of squares in an ANOVA table, not the coefficients. My mistake –  luciano Mar 6 '14 at 16:58
I suspected something like that might be true, @luciano, it's very easy to miss something subtle like that. Let me know if you need anything else. –  gung Mar 6 '14 at 18:23