Here is an introductory course about linear regression, starting with simple linear regression and then moving on to multiple linear regression (menu on the right side). It explains the basics, also mathematical, in more detail: https://online.stat.psu.edu/stat501/
In principle, linear regression tries to create a line through all given data points, where every data point is as close to the line as possible.
Mathematically, the $\beta$-vector of the equation
\begin{align*}
Y &= X \beta + \epsilon \\
&= \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \epsilon
\end{align*}
should minimize the sum of squared residuals RSS (residuals = distance between $Y$ and $\hat{Y}$, where $Y$ are the original Y-values and $\hat{Y}$ the fitted values (the results of the equation $\hat{Y} = X\beta$ where $\beta$ is replaced by the estimated $\hat{\beta}$)). Plotting all data points with the estimated $\hat{Y}$ results in a straight line closest to ALL original points in total (not just optimized for 1 point).
This minimization problem is also called Ordinary-Least-Squares (OLS), a derivation of the equation that finally needs to be solved, can be found here for example:
How to derive the least square estimator for multiple linear regression?
The final equation then is: $\beta = (X^\top X)^{-1}X^\top Y$
X <- model.matrix(model)
head(X)
X is, as pointed out earlier, the model matrix. It is constructed in such a way, that the first column contains a vector of $1$'s followed by the covariate column vectors. In case of categorial variables, for all levels except the reference category, a column vector is created. Using dummy coding for categorial variables $X1$, $X2$ either a $1$ or $0$ is made depending on the level (a, b, c) / (d, e). If all column vectors of a variable are set to 0 within an observation, the value of the variable for this specific observation equals the reference category. Each row contains 1 observation, thus X has 100 rows in this example.
Solving the equation from above with matrix multiplication gets you all $\beta$ estimates at once:
solve(t(X)%*%X)%*%t(X)%*%Y
There are several possibilities how OLS estimates can be derived. This post explains what the function lm() uses:
Algorithm for minimization of sum of squares in regression packages
One way to calculate the estimates "by hand" in detail can be found here (as an example for $\beta_0, \beta_1, \beta_2)$:
http://qed.econ.queensu.ca/pub/faculty/abbott/econ351/351note12.pdf
Adapted to your example this means:
\begin{equation}
\begin{array}{c}
\mathbf{X} =
\begin{pmatrix}
1 & X_{1b,1} & X_{1c,1} & X_{2e,1} & X_{3,1} \\
1 & X_{1b,2} & X_{1c,2} & X_{2e,2} & X_{3,2} \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
1 & X_{1b,n} & X_{1c,n} & X_{2e,n} & X_{3,n} \\
\end{pmatrix}
\end{array}
\quad
\begin{array}{c}
\mathbf{y} =
\begin{pmatrix}
y_1 \\
y_2 \\
\vdots \\
y_n \\
\end{pmatrix}
\end{array}
\\
\end{equation}
\begin{equation}
\mathbf{X}^\top \mathbf{X} =
\begin{pmatrix}
n & \sum_{i=1}^{n} X_{1b,i} & \sum_{i=1}^{n} X_{1c,i} & \sum_{i=1}^{n} X_{2e,i} & \sum_{i=1}^{n} X_{3,i} \\[10pt]
\sum_{i=1}^{n} X_{1b,i} & \sum_{i=1}^{n} X_{1b,i}^2 & \sum_{i=1}^{n} X_{1b,i} X_{1c,i} & \sum_{i=1}^{n} X_{1b,i} X_{2e,i} & \sum_{i=1}^{n} X_{1b,i} X_{3,i} \\[10pt]
\sum_{i=1}^{n} X_{1c,i} & \sum_{i=1}^{n} X_{1c,i} X_{1b,i} & \sum_{i=1}^{n} X_{1c,i}^2 & \sum_{i=1}^{n} X_{1c,i} X_{2e,i} & \sum_{i=1}^{n} X_{1c,i} X_{3,i} \\[10pt]
\sum_{i=1}^{n} X_{2e,i} & \sum_{i=1}^{n} X_{2e,i} X_{1b,i} & \sum_{i=1}^{n} X_{2e,i} X_{1c,i} & \sum_{i=1}^{n} X_{2e,i}^2 & \sum_{i=1}^{n} X_{2e,i} X_{3,i} \\[10pt]
\sum_{i=1}^{n} X_{3,i} & \sum_{i=1}^{n} X_{3,i} X_{1b,i} & \sum_{i=1}^{n} X_{3,i} X_{1c,i} & \sum_{i=1}^{n} X_{3,i} X_{2e,i} & \sum_{i=1}^{n} X_{3,i}^2 \\
\end{pmatrix}
\\
\end{equation}
\begin{equation}
\mathbf{X}^\top \mathbf{y} =
\begin{pmatrix}
\sum_{i=1}^{n} y_i \\
\sum_{i=1}^{n} X_{1b,i} y_i \\
\sum_{i=1}^{n} X_{1c,i} y_i \\
\sum_{i=1}^{n} X_{2e,i} y_i \\
\sum_{i=1}^{n} X_{3,i} y_i \\
\end{pmatrix}
\end{equation}
Solving now for each single $\beta_i$ results in:
\begin{equation}
% For β0
\beta_0 = \frac{ \left( \sum_{i=1}^{n} y_i \right) - \left( \beta_{1b} \sum_{i=1}^{n} X_{1b,i} + \beta_{1c} \sum_{i=1}^{n} X_{1c,i} + \beta_{2e} \sum_{i=1}^{n} X_{2e,i} + \beta_3 \sum_{i=1}^{n} X_{3,i} \right) }{n}
\end{equation}
\begin{equation}
% For β1b
\beta_{1b} = \frac{ \left( \sum_{i=1}^{n} X_{1b,i} y_i \right) - \left( \beta_0 \sum_{i=1}^{n} X_{1b,i} + \beta_{1c} \sum_{i=1}^{n} X_{1b,i} X_{1c,i} + \beta_{2e} \sum_{i=1}^{n} X_{1b,i} X_{2e,i} + \beta_3 \sum_{i=1}^{n} X_{1b,i} X_{3,i} \right) }{\sum_{i=1}^{n} X_{1b,i}^2}
\end{equation}
\begin{equation}
% For β1c
\beta_{1c} = \frac{ \left( \sum_{i=1}^{n} X_{1c,i} y_i \right) - \left( \beta_0 \sum_{i=1}^{n} X_{1c,i} + \beta_{1b} \sum_{i=1}^{n} X_{1c,i} X_{1b,i} + \beta_{2e} \sum_{i=1}^{n} X_{1c,i} X_{2e,i} + \beta_3 \sum_{i=1}^{n} X_{1c,i} X_{3,i} \right) }{\sum_{i=1}^{n} X_{1c,i}^2}
\end{equation}
\begin{equation}
% For β2e
\beta_{2e} = \frac{ \left( \sum_{i=1}^{n} X_{2e,i} y_i \right) - \left( \beta_0 \sum_{i=1}^{n} X_{2e,i} + \beta_{1b} \sum_{i=1}^{n} X_{2e,i} X_{1b,i} + \beta_{1c} \sum_{i=1}^{n} X_{2e,i} X_{1c,i} + \beta_3 \sum_{i=1}^{n} X_{2e,i} X_{3,i} \right) }{\sum_{i=1}^{n} X_{2e,i}^2}
\end{equation}
\begin{equation}
% For β3
\beta_3 = \frac{ \left( \sum_{i=1}^{n} X_{3,i} y_i \right) - \left( \beta_0 \sum_{i=1}^{n} X_{3,i} + \beta_{1b} \sum_{i=1}^{n} X_{3,i} X_{1b,i} + \beta_{1c} \sum_{i=1}^{n} X_{3,i} X_{1c,i} + \beta_{2e} \sum_{i=1}^{n} X_{3,i} X_{2e,i} \right) }{\sum_{i=1}^{n} X_{3,i}^2}
\end{equation}
Here is the R code:
### preparation
X1.b <- X[,2]
X1.c <- X[,3]
X2.e <- X[,4]
X3. <- X[,5]
n <- nrow(X)
b0 <- coef(model)[1]
b1b <- coef(model)[2]
b1c <- coef(model)[3]
b2e <- coef(model)[4]
b3 <- coef(model)[5]
sum.X1b <- sum(X1.b)
sum.X1c <- sum(X1.c)
sum.X2e <- sum(X2.e)
sum.X3 <- sum(X3.)
sum.Y <- sum(Y)
sum.X1b.sqare <- sum(X1.b^2)
sum.X1c.square <- sum(X1.c^2)
sum.X2e.sqare <- sum(X2.e^2)
sum.X3.sqare <- sum(X3.^2)
sum.X1b.Y <- sum(X1.b*Y)
sum.X1c.Y <- sum(X1.c*Y)
sum.X2e.Y <- sum(X2.e*Y)
sum.X3.Y <- sum(X3*Y)
sum.X1b.X1c <- sum(X1.b*X1.c)
sum.X1b.X2e <- sum(X1.b*X2.e)
sum.X1b.X3 <- sum(X1.b*X3.)
sum.X1c.X2e <- sum(X1.c*X2.e)
sum.X1c.X3 <- sum(X1.c*X3.)
sum.X2e.X3 <- sum(X2.e*X3.)
### Solving equations
beta.0 <- unname(sum.Y - (b1b*sum.X1b + b1c*sum.X1c + b2e*sum.X2e + b3*sum.X3))/n
beta.0
beta.1b <- unname((sum.X1b.Y - (b0*sum.X1b + b1c*sum.X1b.X1c + b2e*sum.X1b.X2e + b3*sum.X1b.X3))/sum.X1b.sqare)
beta.1b
beta.1c <- unname((sum.X1c.Y - (b0*sum.X1c + b1b*sum.X1b.X1c + b2e*sum.X1c.X2e + b3*sum.X1c.X3))/sum.X1c.square)
beta.1c
beta.2e <- unname((sum.X2e.Y - (b0*sum.X2e + b1b*sum.X1b.X2e + b1c*sum.X1c.X2e + b3*sum.X2e.X3))/sum.X2e.sqare)
beta.2e
beta.3 <- unname((sum.X3.Y - (b0*sum.X3 + b1b*sum.X1b.X3 + b1c*sum.X1c.X3 + b2e*sum.X2e.X3))/sum.X3.sqare)
beta.3
2 Add-Ons to the great previous post from Thomas:
(1) The reference category can also be changed, it doesn't have to be necessarily the first level of your variables.
For example setting level "b" in X1 as reference category - new $\beta$ estimate for X1a instead of X1b:
X1 <- relevel(X1,ref=2)
model2 <- lm(Y ~ X1 + X2 + X3)
summary(model2)
More options can be found in this post:
https://stackoverflow.com/questions/3872070/how-to-force-r-to-use-a-specified-factor-level-as-reference-in-a-regression
(2) The following code provides a quick overview how the categorial variables are coded. Just pick all columns of categorial variables. It shows the "unique" rows of the model matrix:
unique(model.matrix(model)[,1:4])
model.matrix(model). By default, factors are dummy encoded, dropping the "first" level (as being linearly dependent). $\endgroup$