Recall that the least squares estimates minimize \[SSE(\beta_0, \beta_1) = \sum_{i = 1}^n \left( Y_i - \beta_0 - \beta_1 x_i \right)^2\]
Express this as the inner product \[SSE(\beta) = (Y - X\beta)^T (Y - X\beta)\]
Two approaches to finding the minimizer:
calculus
geometry
A dash of vector calculus
The gradient is the vector equivalent of a derivative.
If \(x \in \mathbb{R}^p\) and \(f(x)\) is a real-valued function the gradient is \[\nabla f(x) = \left[\begin{array}{cccc}
\frac{\partial f}{\partial x_1}
&\frac{\partial f}{\partial x_2}
&\cdots
&\frac{\partial f}{\partial x_p}
\end{array}\right]^T
\quad\text{where}\quad
\frac{\partial f}{\partial x_j} = \frac{d}{dx_j}f(x)\]
Your turn: If \(f(x) = 4x_1 + 5x_2\) then \(\nabla f(x) = \left[\hspace{1in}\right]^T\)
The calculus approach
Let \(A \in \mathbb{R}^{p\times p}\) be a symmetric matrix and \(a, x \in \mathbb{R}^p\) be vectors. Two rules: \[\nabla a^T x = \nabla x^T a = a
\quad\text{and}\quad
\nabla x^T A x = 2Ax\]
# variances of simulated estimatesapply(B, 1, var)
intercept slope
11.66160516 0.08338812
# theoretical variance64*solve(t(X) %*% X)
x
11.2202514 -0.82478585
x -0.8247859 0.07855103
Questions:
Which numbers suggest estimates are unbiased?
What error variance was used in the simulations?
What is -0.8248? How would you approximate it?
Wait, the estimators are correlated?
The plot shows simulated pairs \(\left(\hat\beta_0, \hat\beta_1\right)\).
Each point is one estimate.
cor(B['slope', ], B['intercept', ])
[1] -0.883547
Why would this happen?
Centering the predictor
If we instead fit
\[Y_i = \beta_0 + \beta_1 (x_i - \bar{x}) + \epsilon_i\] then \(\beta_0\) represents the mean of \(Y\) at \(\bar{x}\) and the estimators become uncorrelated.
and it is straightforward to show (HW2) that \[X^T X = \left[\begin{array}{cc}
n &0 \\
0 &(n - 1)s^2_x
\end{array}\right]
\quad\text{and thus}\quad
(X^T X)^{-1} = \left[\begin{array}{cc}
\frac{1}{n} & 0 \\
0 & \frac{1}{(n - 1)s^2_x}
\end{array}\right]\]
Centering the predictor
For example, centering age in the RFFT model:
# fit with centered predictorfit <-lm(rfft ~I(age -mean(age)), data = prevend)# inspectfit
So far we know that \(\mathbb{E}[ \hat{\beta} ] = \beta\) and \(\text{Var}[\hat{\beta}] = \sigma^2 (X^T X)^{-1}\). Could we do better?
We can’t improve on the bias, so “better” here would mean smaller variance.
Writing \(\hat{\beta} = AY\) with \(A = (X^TX)^{-1}X^T\), we see it is a linear function of \(Y\).
Gauss-Markov theorem. For any unbiased estimator of the form \(AY\), we have that \(\text{var}[AY] \succeq \text{var}[\hat{\beta}]\) under the model \(Y = X\beta + \epsilon\) with \(\mathbb{E}[\epsilon] = 0\) and \(\text{Var}(\epsilon) = \sigma^2 I\).
For a proof, see Rencher & Schaalje, Linear Models in Statistics.
So the least squares estimator is the best linear unbiased estimator, or BLUE for short.
Fitted values
We call \(\hat{Y} = X\hat\beta\) the “fitted values”.
Observe that \[\hat{Y} = \underbrace{X (X^TX)^{-1}X^T}_{H}\; Y\] where \(H\) is an orthogonal projection.
The hat matrix
We call \(H\) the “hat matrix” because it puts a hat on \(Y\): \[\hat{Y} = HY\]
\(H\) is an orthogonal projection in the formal sense, meaning it is
symmetric: \(H^T = H\)
idempotent: \(HH = H\)
The residuals can be written compactly using the hat matrix: \[e = Y - HY = (I - H)Y\]
and \(I - H\) is also a projection (see handout), specifically onto the orthogonal complement of \(\text{col}(X)\) in \(\mathbb{R}^n\), which we denote by \(\text{col}(X)^\perp\).
Estimating \(\sigma^2\)
Variance is always a little trickier.
Now since \((I - H)X = 0\), the residual vector is a projection of \(\epsilon\) onto \(\text{col}(X)^\perp\):
\[e = (I - H)Y = (I - H)(X\beta + \epsilon) = (I - H) \epsilon\]
which we can re-express as \[(I - H)\,\epsilon = \sum_{j} (u_j^T\, \epsilon)\, u_j\] where \(u_1, \dots, u_{n-2}\) is an orthonormal basis.
Estimating \(\sigma^2\)
Now the residual variance is proportional to \(\|e\|^2\). As a matter of algebra \[\|e\|^2 = [(I - H)\,\epsilon]^T [(I - H)\, \epsilon] = \sum_{j} (u_j^T\,\epsilon)^2\]
Now since \(\mathbb{E}[u_j^T\epsilon] = 0\) we have that \(\text{var}[u_j^T\epsilon] = \mathbb{E}[(u_j^T \epsilon)^2]\), and then \[\mathbb{E}\left[\|e\|^2\right]
= \sum_j \text{var}\left(u_j^T\,\epsilon\right)
= \sum_j u_j^T\, (\sigma^2 I)\, u_j
= \sum_j \sigma^2 u_j^T u_j
= (n - 2)\, \sigma^2\]
And finally, we obtain an unbiased estimate of the residual variance \[\hat{\sigma}^2 = \frac{1}{n - 2}\, \|e\|^2\]
Uncertainty quantification
All told, we now have
least squares estimates \(\hat{\beta} = (X^T X)^{-1} X^T Y\)
Put these together and we get point estimates \(\hat{\beta}\) with standard errors given by \[\widehat{\text{var}}\left(\hat{\beta}\right) = \hat{\sigma}^2 (X^T X)^{-1}\]
So far we have that under the model \[Y = X\beta + \epsilon\,,
\quad
\mathbb{E}[\epsilon] = 0\,,
\quad
\text{Var}[\epsilon] = \sigma^2 I\] the least squares estimates are BLUE with \(\text{var}(\hat\beta) = \sigma^2 (X^T X)^{-1}\).
If in addition we assume the errors are normal, we get that \[\epsilon \sim N(0, \sigma^2 I)
\quad\Longrightarrow\quad
\hat\beta \sim N(\beta\,, \sigma^2 (X^T X)^{-1})\]
This provides a basis for statistical inference on \(\beta\)
