As a prequel to a question about linear-mixed models in R, and to share as a reference for beginner/intermediate statistics aficionados, I decided to post as an independent "Q&A-style" the steps involved in the "manual" computation of the coefficients and predicted values of a simple linear regression.

The example is with the R in-built dataset, mtcars,  and would be set up as miles per gallon consumed by a vehicle acting as the independent variable, regressed over the weight of the car (continuous variable), and the number of cylinders as a factor with three levels (4, 6 or 8) without interactions.

EDIT: If you are interested in this question, you will definitely find a detailed and satisfactory answer in this post by Matthew Drury outside CV.

Note: I've posted an expanded version of this answer on my website. 


  Would you kindly consider posting a similar answer with the actual R engine exposed?


Sure!  Down the rabbit hole we go.

The first layer is lm, the interface exposed to the R programmer.  You can look at the source for this by just typing lm at the R console.  The majority of it (like the majority of most production level code) is busy checking of inputs, setting of object attributes, and throwing of errors; but this line sticks out

lm.fit(x, y, offset = offset, singular.ok = singular.ok, 
                ...)


lm.fit is another R function, you can call it yourself.  While lm conveniently works with formulas and data frame, lm.fit wants matrices, so that's one level of abstraction removed.  Checking the source for lm.fit, more busywork, and the following really interesting line

z <- .Call(C_Cdqrls, x, y, tol, FALSE)


Now we are getting somewhere.  .Call is R's way of calling into C code.  There is a C function, C_Cdqrls in the R source somewhere, and we need to find it.  Here it is.

Looking at the C function, again, we find mostly bounds checking, error cleanup, and busy work.  But this line is different

F77_CALL(dqrls)(REAL(qr), &n, &p, REAL(y), &ny, &rtol,
        REAL(coefficients), REAL(residuals), REAL(effects),
        &rank, INTEGER(pivot), REAL(qraux), work);


So now we are on our third language, R has called C which is calling into fortran.  Here's the fortran code.

The first comment tells it all

c     dqrfit is a subroutine to compute least squares solutions
c     to the system
c
c     (1)               x * b = y


(interestingly, looks like the name of this routine was changed at some point, but someone forgot to update the comment).  So we're finally at the point where we can do some linear algebra, and actually solve the system of equations.  This is the sort of thing that fortran is really good at, which explains why we passed through so many layers to get here. 

The comment also explains what the code is going to do

c     on return
c
c        x      contains the output array from dqrdc2.
c               namely the qr decomposition of x stored in
c               compact form.


So fortran is going to solve the system by finding the $QR$ decomposition.

The first thing that happens, and by far the most important, is

call dqrdc2(x,n,n,p,tol,k,qraux,jpvt,work)


This calls the fortran function dqrdc2 on our input matrix x.  Whats this?  

 c     dqrfit uses the linpack routines dqrdc and dqrsl.


So we've finally made it to linpack.  Linpack is a fortran linear algebra library that has been around since the 70s.  Most serious linear algebra eventualy finds its way to linpack.  In our case, we are using the function dqrdc2

c     dqrdc2 uses householder transformations to compute the qr
c     factorization of an n by p matrix x.


This is where the actual work is done.  It would take a good full day for me to figure out what this code is doing, it is as low level as they come.  But generically, we have a matrix $X$ and we want to factor it into a product $X = QR$ where $Q$ is an orthogonal matrix and $R$ is an upper triangular matrix.  This is a smart thing to do, because once you have $Q$ and $R$ you can solve the linear equations for regression

$$ X^t X \beta = X^t Y $$

very easily.  Indeed

$$ X^t X = R^t Q^t Q R = R^t R $$

so the whole system becomes

$$ R^t R \beta = R^t Q^t y $$

but $R$ is upper triangular and has the same rank as $X^t X$, so as long as our problem is well posed, it is full rank, and we may as well just solve the reduced system

$$ R \beta = Q^t y $$

But here's the awesome thing.  $R$ is upper triangular, so the last linear equation here is just constant * beta_n = constant, so solving for $\beta_n$ is trivial.  You can then go up the rows, one by one, and substitute in the $\beta$s you already know, each time getting a simple one variable linear equation to solve.  So, once you have $Q$ and $R$, the whole thing collapses to what is called backwards substitution, which is easy.  You can read about this in more detail here, where an explicit small example is fully worked out.
