Order of variables in R's lm

Question

I dont quite understand the answer given in Order of variables in R lm model in lm() function of R (and generally formulas) why changing the order of variable matters? My own guess is that the model first calculates the effect of first variable, and then uses the second variable for remaining variation in dependent variable and so on.

set.seed(1)
a = seq(1:100) + rnorm(100, sd=5)
    b = seq(0.01:0.99, by=0.01) + rnorm(100,sd=3)/100
    c = seq(1:100) + rnorm(100,sd=3)
    d = seq(1:100) + rnorm(100,sd=3)
    summary(lm(a~c+b+d))
    summary(lm(a~b+c+d))

Answer 1

Differences due to type I/II/III sums

The order is not important for the summary of the linear model (which is based on t-tests that don't change). You can see this in your output which is the same.

However when you do an ANOVA then you might get different results depending on the order (this happens for type I sums)

> anova(lm(a~c+b+d))
Analysis of Variance Table

Response: a
          Df Sum Sq Mean Sq   F value    Pr(>F)    
c          1  82067   82067 3412.9019 < 2.2e-16 ***
b          1    494     494   20.5397 1.683e-05 ***
d          1     77      77    3.1872   0.07738 .  
Residuals 96   2308      24                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(lm(a~b+c+d))
Analysis of Variance Table

Response: a
          Df Sum Sq Mean Sq   F value    Pr(>F)    
b          1  82146   82146 3416.2075 < 2.2e-16 ***
c          1    414     414   17.2341 7.155e-05 ***
d          1     77      77    3.1872   0.07738 .  
Residuals 96   2308      24                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 

Note the different p-values for the factors b and c.

The reason is that ANOVA is a comparison of models and there are different ways to interpret this comparison (see type I/II/III sums).

The standard anova function is performing the models in a cascading way, dropping terms one by one starting from the back. Those are type I sums

It goes a bit like this (but slightly different F-scores because the degrees of freedom are computed differently)

anova(lm(a~1+b+c), lm(a~1+b+c+d))  # testing the effect of d
anova(lm(a~1+b  ), lm(a~1+b+c  ))  # testing the effect of c
anova(lm(a~1    ), lm(a~1+b    ))  # testing the effect of b

The t-scores and the related p-values (from the summary of the lm function) relate to the F-test/ANOVA in the case of type III sums, which is dropping terms relative to the full model (and that is why the order doesn't matter for the t-test)

anova(lm(a~1+b+c), lm(a~1+b+c+d))  # testing the effect of d
anova(lm(a~1+b+d), lm(a~1+b+c+d))  # testing the effect of c
anova(lm(a~1+c+d), lm(a~1+b+c+d))  # testing the effect of b

This can also be done with the drop1 function

> drop1(lm(a~b+c+d), test = "F")
Single term deletions

Model:
a ~ b + c + d
       Df Sum of Sq    RSS    AIC F value   Pr(>F)   
<none>              2308.4 321.92                    
b       1   232.725 2541.2 329.52  9.6783 0.002456 **
c       1   147.721 2456.2 326.12  6.1433 0.014937 * 
d       1    76.639 2385.1 323.18  3.1872 0.077377 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Differences due to the position of the intercept

In the referenced question the reason for the impact of the order is due to the position of the intercept. In that question, the intercept is explicitly excluded from the model, but indirectly it is still part of the model because the categorical variables often add up to one. In that case, the intercept is placed for whichever variable and factor level is first in the order.

An illustrative problems that shows how the factors of $n$ levels are only fitted with $n-1$ coefficients and one term gets absorbed into the intercept is: Fitting a Logistic Regression Without an Intercept . In that problem you see that the person who ask's the questions tries to get rid of this 'dropping of one level for each factor' by not using an intercept. But this only works for one factor. The factor for which this works is the one which is the beginning of the model.

A more silly example is: Why do output coefficients not resemble true coefficients in a linear model? In that particular example pay especially attention to the nls model where the dropping of the first level of each factor must be done explicitly

modelnls2 <- nls(Y ~ exp(a + c(0,b1,b2)[Home] + c(0,c1)[Gender] + c(0,d1,d2)[Rank]), 
                start = c(a=1,b1=1,b2=1,c1=1,d1=1,d2=1), data=dune)

For each of the factors (Home, Gender, Rank), the related coefficients are set explicitly at 0 for one of the levels. If you would take away the intercept coefficient a then you could add it to one of the others. For instance:

modelnls2 <- nls(Y ~ exp(c(a,b1+a,b2+a)[Home] + c(0,c1)[Gender] + c(0,d1,d2)[Rank]), 
                start = c(a=1,b1=1,b2=1,c1=1,d1=1,d2=1), data=dune)
### or equivalent
modelnls2 <- nls(Y ~ exp(c(a,b1,b2)[Home] + c(0,c1)[Gender] + c(0,d1,d2)[Rank]), 
                start = c(a=1,b1=1,b2=1,c1=1,d1=1,d2=1), data=dune)

This is what happens with the lm function for the referenced question and is the reason why the order matters.

Answer 2

Here is a more formal answer (more elegant proofs that start with something like "consider the space spanned by the columns of $X$..." are surely possible) to the question of why just changing the order of the regressors does not matter. As I sketch at the end, the question you link to deals with a case where more than just permuting the columns is happening (as the accepted answer there, I believe, also explains quite well).

Consider changing the position of the variables in the $(n\times p)$ regressor matrix $$ X=(X_1,\ldots,X_p), $$ where $X_j=(x_{1j},\ldots,x_{nj})'$, $j=1,\ldots,p$, amounts to postmultiplying $X$ with a $(p\times p)$ permutation matrix $P$ that has a single entry 1 in each column $j$ that indicates the new column position of that regressor $X_j$.

For example, if the new columns are to be the old columns 2, 1 and 3, we have $$P=\begin{pmatrix} 0&1&0\\ 1&0&0\\ 0&0&1 \end{pmatrix}$$

This matrix $P$ is invertible, being just a permuted version of the identity matrix.

Also, as it is easy to check that $P'P=I$, $P^{-1}$ is equal to the transpose of $P$, $P^{-1}=P'$.

Thus, the OLS coefficient of the regression of $y$ on the transformed regressors, call it $\hat\beta_t$, is \begin{eqnarray*} \hat\beta_t&=&((XP)'XP)^{-1}(XP)'y\\ &\stackrel{(AB)'=B'A'}{=}&(P'X'XP)^{-1}P'X'y\\ &\stackrel{(ABC)^{-1}=C^{-1}B^{-1}A^{-1}}{=}&P^{-1}(X'X)^{-1}\underbrace{(P')^{-1}P'}_{=I}X'y\\ &=&P^{-1}(X'X)^{-1}X'y\\ &=&P'(X'X)^{-1}X'y\\ &=&P'\hat\beta\\ \end{eqnarray*} Here, $P'$ is a matrix that permutes the row elements of $\hat\beta$, and hence permutes the entries of the original coeffient estimator $\hat\beta$ according to the permutation of the columns.

To see why the question you link to addresses a slightly different situation in which something else happens than just permuting the columns, I suggest to inspect model.matrix(out_1) and model.matrix(out_2) in that code, which gives you the different regressor matrices in the two models.

Answer 3

The two models you fitted are exactly the same. The results are shown in the order you have entered them into R. The only 'difference' as such is that the rows for b and c have been swapped over, but the p-values, estimates and so on are identical

> summary(lm(a~c+b+d))

Call:
lm(formula = a ~ c + b + d)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.6989  -3.4114  -0.0175   3.4746  12.1792 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   0.8080     0.9928   0.814  0.41778   
c             0.3388     0.1367   2.479  0.01494 * 
b            40.5137    13.0228   3.111  0.00246 **
d             0.2504     0.1402   1.785  0.07738 . 

> summary(lm(a~b+c+d))

Call:
lm(formula = a ~ b + c + d)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.6989  -3.4114  -0.0175   3.4746  12.1792 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   0.8080     0.9928   0.814  0.41778   
b            40.5137    13.0228   3.111  0.00246 **
c             0.3388     0.1367   2.479  0.01494 * 
d             0.2504     0.1402   1.785  0.07738 . 
```

Answer 4

Couple of additional points.

1. your understanding of the calculation process is incorrect.

My own guess is that the model first calculates the effect of first variable, and then uses the second variable for remaining variation in dependent variable and so on.

No, Linear Regression calculating all the coefficients at the same time, but one by one. In R, lm is using QR decomposition. And changing the order is just switching columns in the matrix.

2. When people say, order of the variable it may mean "Coding Systems for Categorical Variables in Regression" and reference level for categorical Variables.

For example if you are fitting a regression model, and one of the variable is education level. (high school, bachelor, master), If the reference level is different. Then the coefficient will be different.