Deriving the regression parameters for simple linear regression: linear algebra

Andy Long
Spring, 2020

Clear["Global`*"] (* This removes any definitions that may be in place already; like if you've already been doing some Mathematica, and talking about other functions named f, x-values named a and b, etc. *)

To illustrate the derivation of the regression parameters that minimize the squares of the residuals, we're going to be using the Keeling data from Stewart, p. 25. 

We've already seen how we can find the regression line using Mathematica's  LinearModelFit, and multivariate calculus.

years={0,2,4,6,8,10,12,14,16,18,20,22,24,26,28}
co2={338.7,341.2,344.4,347.2,351.5,354.2,356.3,358.6,362.4,366.5,369.4,373.2,377.5,381.9,385.6}
n=Length[co2]

We want to find a linear function to fit the co2 values: we need to find the parameters a and b:

y(t)=a+bt

params={a ,b};
p=Length [params];
X=Table[{1,years[[i]]},{i,1,n}];
{MatrixForm[co2],"=", MatrixForm[X] ,MatrixForm[ params]}
{MatrixForm[co2],"=", MatrixForm[X . params]}

So we have 15 equations in 2 unknowns, which is called “over-determined” -- the chance of us finding an exact solution is essentially 0. So we relax, and try to turn our problem into a two-equation in two-unknown problem.

We will need the transpose of the matrix X as well:

Xt=Transpose[X];
MatrixForm[Xt]

The product XtX will be square, and 2x2:

{MatrixForm[Xt],MatrixForm[X]}
XtX=Xt.X;
MatrixForm[XtX]

The inverse of this matrix is easy (since it’s 2x2):

XtXinv=Inverse[XtX];
MatrixForm[XtXinv]

Now we find the values of a and b that give the least squares line, from
Xt co2 = Xt X beta
or
XtXinv  Xt  co2 = beta

params=XtXinv .( Xt. co2);
MatrixForm[params]

All your favorite diagnostics:

First we take a look at the fit and the residuals:

(* let's call co2 by our favorite name, y *)
y=co2;
p1=ListPlot[Transpose[{years, y}]];
p2=Plot[{1,x}.params,{x,0,28}];
Show[p1,p2,Frame->True]

(* these are the estimated values of y: *)
yhat=X.params; 
residuals=y-yhat;
ListPlot[Transpose[{years, residuals}]]
Histogram[residuals,4]

Then we use some matrix multiplication to compute some of these quantities:

ybar=Mean[y]; (* Null model *)
(* SSTO is the sum of squared errors of the null model, using the mean of y as the model (a horizontal line): *)
SSTO=(y-ybar).(y-ybar)
(* SSE is the minimal sum of squared errors, the minimum value possible of the sum of squared errors we sought to minimize: *)
SSE=(y-yhat).(y-yhat)
(* SSEreg is the variance of the residuals using the regression equation. *)
SSEreg=(yhat-ybar).(yhat-ybar)
(* MSE is the mean squared error -- the average error of a residual *)
MSE=SSE/(n-p)
(* this is what's called the unadjusted RSquared: how much of the variation in the data is explained by the regression model? *)
r2=SSEreg/SSTO

Finally we refer to the matrix inverse to compute t-statistics or confidence intervals for the parameters:

(* the square roots of the diagonals of the inverse matrix multiplied by the MSE give the standard errors of the parameters: *) 
ses= Sqrt[Diagonal[XtXinv*MSE]] 
(* the t-vals are parameters divided by these: *)
tvals=params/ses
(* and the associated p-values Mathematica gives are the two-sided probabilities of a value of t this large or larger: 
*)
2(1-CDF[StudentTDistribution[n-p],tvals])
(* to get the confidence intervals, we use the same standard errors ses, and add and subtract the appropriate width -- the margin of error -- corresponding to the confidence required (here 95%). 
I'm using a kludgy way to find the critical value of t -- there must be an easier way in Mathematica! *)
conf=.95;
width=Solve[CDF[StudentTDistribution[n-p],x]==1-(1-conf)/2,x,Reals][[1]]
width=x/.width
(* so here are the confidence intervals: *)
Transpose[{params-width*ses,params+width*ses}]