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.

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()