Linear models: estimating the b vector on R

Question

Starting from the basic linear model problem: $$ y=Xb+e $$ And the least squares method of $b$ estimation $$ \hat{b}=(X'X)^{-1}X'y $$ And also considering a model not of full rank (an one-way design, for example), I am having a big issue with R. I can get the SS and residuals estimations fine, but I need the $\hat{b}$ vector. I am aware that the R's lm function uses the qr-decomposition for fitting these models, and I need the pure $\hat{b}$ vector, with the factor effects, such as in:

$$ \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \\ y_7 \\ y_8 \\ \end{bmatrix} = \begin{bmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} \mu \\ \tau_1 \\ \tau_2 \\ \tau_3 \\ \end{bmatrix} + \begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \\ \epsilon_4 \\ \epsilon_5 \\ \epsilon_6 \\ \epsilon_7 \\ \epsilon_8 \\ \end{bmatrix} $$ Therefore, the $effects and the $coefficients attributes of the response lm object are not the desired vector, since some levels get missing ($\tau_1$ for example) and the values do not match the ones calculated manually using the least squares method.

For example: given a 2-factor completely randomized design with 4 levels of A, 2 levels of B and respective interactions, the $\hat{b}$ vector estimated by the described method (simple linear algebra, method described above) is

> b
             [,1]
 [1,]  26.1666667
 [2,]  11.5000000
 [3,]   0.2777778
 [4,]   2.1666667
 [5,]  12.2222222
 [6,]  16.2666667
 [7,]   9.9000000
 [8,]  -3.9333333
 [9,]  15.4333333
[10,] -11.0444444
[11,]  11.3222222
[12,]   5.4000000
[13,]  -3.2333333
[14,]  25.8444444
[15,] -13.6222222

But R's lm()$coefficients are

> mf$coefficients
  (Intercept)           X12           X13           X14           X22 
 5.000000e+01 -1.833333e+01 -1.751311e-14  3.050000e+01  1.300000e+01 
      X12:X22       X13:X22       X14:X22 
 3.000000e+00 -2.800000e+01 -5.883333e+01

And R's lm()$effects are

> mf$effects
  (Intercept)           X12           X13           X14           X22 
-2.147264e+02 -2.083998e+01 -2.391722e+01 -3.478505e+00  1.478977e+01 
      X12:X22       X13:X22       X14:X22                             
 3.287056e+01  1.267105e+00 -4.557210e+01 -9.818222e+00 -1.818222e+00 

-3.818222e+00 -7.380306e+00  2.619694e+00 -6.007701e+00 -7.700896e-03 

-7.700896e-03 -4.365427e-01  4.563457e+00 -7.572086e-01 -2.757209e+00 

-2.757209e+00 

Could anyone please help me obtaining this $\hat{b}$ vector on R?

R code for reproducing the calculations:

dt <- structure(list(X1 = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L), .Label = c("1", 
"2", "3", "4"), class = "factor"), X2 = structure(c(1L, 1L, 1L, 
2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 
2L, 2L), .Label = c("1", "2"), class = "factor"), X3 = structure(c(58, 
45, 47, 61, 65, 31, 35, 29, 43, 51, 49, 45, 55, 31, 37, 37, 78, 
83, 36, 34, 34), .Dim = c(21L, 1L))), .Names = c("X1", "X2", 
"X3"), row.names = c(NA, -21L), class = "data.frame")

X <- structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 
1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 
1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 
0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 1, 1, 1), .Dim = c(21L, 15L))

y <- structure(c(58, 45, 47, 61, 65, 31, 35, 29, 43, 51, 49, 45, 55, 
31, 37, 37, 78, 83, 36, 34, 34), .Dim = c(21L, 1L))

XtX <- structure(c(21, 5, 6, 5, 5, 10, 11, 3, 2, 3, 3, 2, 3, 2, 3, 5, 
5, 0, 0, 0, 3, 2, 3, 2, 0, 0, 0, 0, 0, 0, 6, 0, 6, 0, 0, 3, 3, 
0, 0, 3, 3, 0, 0, 0, 0, 5, 0, 0, 5, 0, 2, 3, 0, 0, 0, 0, 2, 3, 
0, 0, 5, 0, 0, 0, 5, 2, 3, 0, 0, 0, 0, 0, 0, 2, 3, 10, 3, 3, 
2, 2, 10, 0, 3, 0, 3, 0, 2, 0, 2, 0, 11, 2, 3, 3, 3, 0, 11, 0, 
2, 0, 3, 0, 3, 0, 3, 3, 3, 0, 0, 0, 3, 0, 3, 0, 0, 0, 0, 0, 0, 
0, 2, 2, 0, 0, 0, 0, 2, 0, 2, 0, 0, 0, 0, 0, 0, 3, 0, 3, 0, 0, 
3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 3, 0, 3, 0, 0, 0, 3, 0, 0, 0, 3, 
0, 0, 0, 0, 2, 0, 0, 2, 0, 2, 0, 0, 0, 0, 0, 2, 0, 0, 0, 3, 0, 
0, 3, 0, 0, 3, 0, 0, 0, 0, 0, 3, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, 
0, 0, 0, 0, 0, 2, 0, 3, 0, 0, 0, 3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 
3), .Dim = c(15L, 15L))

Xty <- structure(c(984, 276, 238, 205, 265, 506, 478, 150, 126, 95, 
143, 100, 105, 161, 104), .Dim = c(15L, 1L))

myfit <- lm(X3 ~ X1 * X2, data=dt)

Answer 1

After receiving clarification on what you are trying to do through the chat linked hereto, here is the consolidated answer:

myfit <- lm(X3 ~ X1 * X2, data=dt)
newx <- model.matrix(myfit)
solve(t(newx) %*% newx) %*% t(newx) %*% y

Essentially, you need to be sure that you have specified the correct parametrization for the model and correctly specified the model. You can obtain the same parameterization as R using the "set first to zero" parameterization. Take a look at the model matrix produced and stored in newx above. This is how your model matrix must be specified in order to solve for $\hat{b}$ and replicate R's coefficients output.