Difference in multiple regression with x=factor/coefficients matrix in R

Question

I have got 2 models. The first model uses a factor as the independent variable. The 2nd model uses a matrix with 0/1 coefficients as the independent variable.

# fit a model where an iris' petal width depends on species
# first up with a categorical variable
summary(lm(Petal.Width ~ 0+Species, data=iris))

# 2nd with a matrix indicating the species
smat = model.matrix(~ 0+Species, data=iris)
summary(lm(Petal.Width ~ 0+smat, data=iris))

The outcome of the two models above is the same. This is what I would expect:

Call:
lm(formula = Petal.Width ~ 0 + smat, data = iris)

Residuals:
   Min     1Q Median     3Q    Max 
-0.626 -0.126 -0.026  0.154  0.474 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
smatSpeciessetosa      0.24600    0.02894    8.50 1.96e-14 ***
smatSpeciesversicolor  1.32600    0.02894   45.82  < 2e-16 ***
smatSpeciesvirginica   2.02600    0.02894   70.00  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2047 on 147 degrees of freedom
Multiple R-squared:  0.9796,    Adjusted R-squared:  0.9792 
F-statistic:  2357 on 3 and 147 DF,  p-value: < 2.2e-16

However, when the model gets more complicated the summary is not the same anymore:

x = c(rep("x1",19),rep("x2",23),rep("x3",8),rep("x4",101),rep("x5",66),rep("x6",26),rep("x7",88))
y = c(0.74, 0.37, -0.34, 0.18, 0.18, -1.2, 0.04, 0.61, 0.85, 0, -0.03,  -0.18, -0.12, 1.37, -0.34, -0.4, -1.71, -0.13, 3.63, 1.6, 0.47,  0.68, 0.14, 0.86, 0.68, 0.07, -0.33, 4.14, 5.3, 0.94, -0.54,  -2.14, 0.12, 0.81, 1.51, -1.53, 3.37, -1.83, -6.99, -4.69, 1.59,  -1.01, -0.09, 0.04, 0.01, 20.35, -0.28, -1.67, -1.08, -5.33,  1.26, 4.35, 0.07, 1.55, 0.65, 0.15, -0.08, -0.31, -1.61, 0.57,  -0.36, 0.47, 0.72, -0.66, 0.92, 0.77, -1.93, -1.06, -0.93, 0.33,  0, 1.23, -1.53, -0.15, -3.5, 1.01, 0.63, 0.98, 1.06, 0.77, 0.33,  0.72, 1.56, 1.81, 0.74, 1.06, 0.09, 0.15, 0.75, 0.17, 0.77, 0.63,  0.54, 0.22, 0.49, -0.81, -0.03, 2.19, 0.76, 0.11, 0.94, -1.7,  -2.11, 0.32, -4.9, 0.39, 1.59, -5.66, 2.08, -10.43, -0.08, -0.9,  -0.16, -0.86, 0.43, -2.37, -2.16, -0.5, -0.6, 0.27, -0.34, 1.41,  -0.42, -1.43, -2.25, -0.69, -0.71, -1.54, -3.12, -0.93, -2.01,  0.26, -0.53, -1.09, -1.45, -2.75, -1.63, 0.16, -0.13, 0.08, -2.93,  -1.46, 0.18, -0.45, -0.88, 0.05, 0.68, -0.62, -2.55, -3.33, -2.35,  -3.56, -0.07, -0.29, 1.81, 1.07, 0.01, 0.85, -0.21, -1.68, 1.13,  1.58, -0.31, 1.14, -0.73, 2.04, 0.82, 0.78, 0.69, 1.18, 3.78,  0.97, 0.7, 13.63, -2.35, 4.52, 7.46, 4.09, 1.36, 0.77, -0.51,  0.09, 3.12, -3.36, -1.04, 2.79, 0.1, 2.83, 1.96, 0.15, -0.09,  0.39, 0.65, -1.7, -0.84, 1.86, -0.6, 0.18, 0.49, -0.66, -0.16,  2.16, 0.46, 1.66, 0.13, 2.12, -0.03, -0.31, 0.13, -0.45, 0.18,  -0.17, 0.25, -0.16, -0.06, 0.01, 0.83, -2.04, 1.06, 0.35, 0.86,  0.8, 2.01, 2.55, 0.07, 0.55, 1.59, -2.79, -2.77, 1.15, 0.46,  0.45, -0.4, 1.54, 0.29, 0.17, -0.36, 0.99, -9.59, 2.6, 1.79,  -0.56, -0.49, 0.28, 1.46, -1.19, -0.48, -2.52, 1.08, -4.1, -2.91,  -4.89, -0.12, -10.36, -4.97, -7.07, -3.81, -4.34, -12.82, -2.51,  -8.47, -6.12, -1.56, -4.48, -2.88, -1.94, -6.01, -1.35, 1.23,  3.08, -2.3, -0.99, 0.63, 2.4, 2.27, -0.79, 1.76, 1.21, -4.43,  -1.34, -1.94, -0.57, -2.77, -3.07, -1.15, -2.28, 0.68, 1.5, -1.78,  0.46, 0.39, -2.04, -1.22, 0.66, 1.22, 2.61, 5.21, 1.17, 2.73,  2.48, -0.95, 1.59, 0.14, -2.62, -1.73, -0.55, -1.25, 0.15, -0.11,  -0.23, -3.47, 1.46, -1.32, 0.24, 1.48, 2.68, 0.07, 0.24, -0.76,  -1.18, -1.57, -4.08, -2.74, 0.03, 0.34, 0.34, 0.84, 1.26, -0.81,  -0.79, -1.59)

with the following output:

> summary(lm(y ~ 0+x))

Call:
lm(formula = y ~ 0 + x)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.6615  -1.0202   0.0515   1.1498  18.8562 

Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
xx1  0.18526    0.59168   0.313    0.754    
xx2  0.14000    0.53778   0.260    0.795    
xx3  1.49375    0.91184   1.638    0.102    
xx4 -0.41149    0.25663  -1.603    0.110    
xx5  0.81182    0.31746   2.557    0.011 *  
xx6  0.01077    0.50580   0.021    0.983    
xx7 -1.15852    0.27493  -4.214 3.26e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.579 on 324 degrees of freedom
Multiple R-squared:  0.08401,   Adjusted R-squared:  0.06422 
F-statistic: 4.245 on 7 and 324 DF,  p-value: 0.0001676

> summary(lm(y ~ 0+model.matrix(~x)))

Call:
lm(formula = y ~ 0 + model.matrix(~x))

Residuals:
     Min       1Q   Median       3Q      Max 
-11.6615  -1.0202   0.0515   1.1498  18.8562 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)  
model.matrix(~x)(Intercept)  0.18526    0.59168   0.313   0.7544  
model.matrix(~x)xx2         -0.04526    0.79956  -0.057   0.9549  
model.matrix(~x)xx3          1.30849    1.08699   1.204   0.2296  
model.matrix(~x)xx4         -0.59675    0.64494  -0.925   0.3555  
model.matrix(~x)xx5          0.62656    0.67147   0.933   0.3515  
model.matrix(~x)xx6         -0.17449    0.77841  -0.224   0.8228  
model.matrix(~x)xx7         -1.34379    0.65244  -2.060   0.0402 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.579 on 324 degrees of freedom
Multiple R-squared:  0.08401,   Adjusted R-squared:  0.06422 
F-statistic: 4.245 on 7 and 324 DF,  p-value: 0.0001676

Why are the p-values of the coefficients different here?

If it's the order in which the coefficients are fitted, is there a way to get the same behaviour?

Answer 1

I forgot that the command

model.matrix(~x)

will fit coefficients with intercept, which I didn't want. The right syntax is

model.matrix(~0+x)

and then the results are the same.

So, to sum up, the following 2 models behave the same way:

summary(lm(y ~ 0+x))
summary(lm(y ~ 0+model.matrix(~0+x))) # notice that '0' appears twice

To get the same behaviour with intercept, the following syntax is needed:

summary(lm(y ~ x))
summary(lm(y ~ model.matrix(~x))) # notice that no '0' is in model.matrix()