### 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][1]).

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][2] the reason for the impact of the order is due to the [position of the intercept][3]. 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.
 
(There must be some illustrative problems that show how the factors of $n$ levels are only fitted with $n-1$ coefficients and one term gets absorbed into the intercept. I have seen such problems several times, but I can not find them back. I found one, but it is a bit silly example: https://stats.stackexchange.com/questions/427475/)


  [1]: https://stats.stackexchange.com/search?q=type%20sums
  [2]: https://stats.stackexchange.com/questions/63189/order-of-variables-in-r-lm-model
  [3]: https://stats.stackexchange.com/a/308644/