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