Odds and Ends from the Theory of Linear Regression

1 MLE Properties of Least-Squares

1.1 The Multivariate Normal Density

The univariate normal density (n = 1) we all remember is

φ(y) = \frac{1}{\sqrt{2π}σ} exp [−\frac{1}{2} (\frac{y − μ}{σ})^{2}]

for y ∈ R. We often write Y ∼ N(μ, σ^{2}) to convey that the random variable Y has this density.

The multivariate normal density, for \underline{y} ∈ R^{n}, is

φ(\underline{y}) = \frac{1}{(2π)^{n/2}|Σ̸|^{1/2}} exp[−\frac{1}{2}(\underline{y} − \underline{μ})^{T}Σ̸^{−1}(\underline{y} − \underline{μ})]  (1)

The mean vector \underline{μ} replaces the univariate mean μ and the variance-covariance matrix Σ̸ replaces the univariate variance σ^{2}. We often write \underline{Y} ∼ N_{n}(\underline{μ},Σ̸), or just \underline{Y} ∼ N(\underline{μ},Σ̸) if the dimension n is clear from context. Sometimes we also say that the components of \underline{Y} are “jointly normal” or “multivariate normal”.

One way to build up the multivariate normal is to start with a n-dimensional vector of independent identically distributed N(0, 1) variables \underline{Z} = (Z_{1}, …, Z_{n}). Since densities for independent random variables multiply, the density for \underline{Z} must be

φ(\underline{z}) = ∏_{1}^{n} φ(z_{i})
= \frac{1}{(2π)^{n/2}} exp[−\frac{1}{2} ∑_{1}^{n} z_{i}^{2}]
= \frac{1}{(2π)^{n/2}} exp[−\frac{1}{2}\underline{z}^{T} I^{−1}\underline{z}]

where I is the n × n identity matrix.

Then if \underline{Y} = A\underline{Z} + \underline{μ}, it is clear that E[\underline{Y}] = \underline{μ} and Var(\underline{Y}) = E[(\underline{Y} − \underline{μ})(\underline{Y} − \underline{μ})^{T}] = AA^{T} = Σ̸ (so for example it is immediately obvious that Σ̸ must by symmtric and positive definite).

The multivariate change-of-variables formula from calculus will ensure that the density for \underline{Y} looks like (1). In fact, it is generally true that if \underline{Y} = A\underline{Y} + \underline{b}, where A is m × n and \underline{b} ∈ R^{m}, and \underline{Y} ∼ N_{n}(\underline{μ},Σ̸), then \underline{Y} ∼ N_{m}(A\underline{μ} + \underline{b}, AΣ̸A^{T}).

1.2 Maximum Likelihood Estimation

Now let us derive maximum likelihood estimates for \underline{β} and σ_{ϵ}^{2} in the multiple regression model

\underline{y}_{n×1} = X_{n×k′} \underline{β}_{k′×1} + \underline{ϵ}_{n×1} \underline{ϵ} ∼ N_{n}(\underline{0}, σ_{ϵ}^{2}I_{n×n}).

Since this implies the y_{i}’s are independent

\underline{y} ∼ N_{n}(Xβ, σ_{ϵ}^{2}I_{n×n})

we can immediately write down the likelihood for \underline{β} and σ_{ϵ}^{2} by writing down the product of the individual normals^{1},

L(\underline{y}; \underline{β}, σ_{ϵ}^{2}) = ∏_{1}^{n} \frac{1}{\sqrt{2πσ_{ϵ}^{2}}} exp [−\frac{1}{2}(y_{i} − \underline{x}_{i}^{T} \underline{β})^{2}/σ_{ϵ}^{2}]

= \frac{1}{(2πσ_{ϵ}^{2})^{n/2}} exp [−\frac{1}{2} ∑_{1}^{n}(y_{i} − \underline{x}_{i}^{T} \underline{β})^{2}/σ_{ϵ}^{2}]
= \frac{1}{(2π)^{n/2}|σ^{2}I|^{1/2}} exp [−\frac{1}{2}(\underline{y} − X\underline{β})^{T} (σ_{ϵ}^{2}I)^{−1}(\underline{y} − X\underline{β})] ,

and hence the log-likelihood is

ℓ(\underline{β}, σ_{ϵ}^{2}) = log L((\underline{y}; \underline{β}, σ_{ϵ}^{2}) = −\frac{n}{2} log 2π − \frac{n}{2} log σ_{ϵ}^{2} − \frac{1}{2σ_{ϵ}^{2}} (\underline{y} − X\underline{β})^{T} (\underline{y} − X\underline{β})

Now we differentiate with respect to \underline{β} and σ_{ϵ}^{2} to obtain

\frac{∂ℓ}{∂\underline{β}} = −\frac{1}{2σ_{ϵ}^{2}} [−2X^{T} (y − X\underline{β})]
\frac{∂ℓ}{∂σ_{ϵ}^{2}} = −\frac{n}{2} \frac{1}{σ_{ϵ}^{2}} + \frac{1}{σ_{ϵ}^{4}} (\underline{y} − X\underline{β})^{T} (\underline{y} − X\underline{β})

Setting the first set of partial derivatives equal to zero produces the normal equations again, and hence

ˆ\underline{β} = (X^{T} X)^{−1}X^{T} \underline{y},

the usual least-squares regression estimates. Setting the second set of partial derivatives equal to zero produces

σ̂_{ϵ}^{2} = \frac{(\underline{y} − X ˆ\underline{β})^{T} (\underline{y} − X ˆ\underline{β})}{n} = \frac{RSS}{n}

which we now know is a biased estimate of σ_{ϵ}^{2} (the correct denominator should be n − k′).

1.3 Some Geometry; Eigensystems; Square Roots; Principal Components

A useful geometric property of the multivariate normal distribution is that it is “log-quadratic”; log φ(\underline{y}) is a simple linear function of the symmetric positive-definite quadratic form (\underline{y}−\underline{μ})^{T}Σ̸^{−1}(\underline{y}− \underline{μ}) in (\underline{y} − \underline{μ}). This means that the level surfaces \{\underline{y} : φ(\underline{y}) = K\}, or equivalently the level surfaces \{\underline{y} : log φ(\underline{y}) = c\} will be ellipsoids.

1Warning: In this handout, \underline{x}_{i}^{T} is the i^{th} \underline{row} of the regression matrix X.

Finding the principal axes of the ellipsoidal level sets (contours of φ(\underline{y})) is easy with Lagrange multipliers. To find the first principal axis for example, you need a point \underline{y} that is a maximum distance from \underline{μ} [maximize the squared distance (\underline{y}−\underline{μ})^{T} (\underline{y}−\underline{μ})], but is still on the level set [subject to the constraint (\underline{y} − \underline{μ})^{T}Σ̸^{−1}(\underline{y} − \underline{μ}) = c]. Differentiating the Lagrange-multiplier objective function

f(\underline{y}) = (\underline{y} − \underline{μ})^{T} (\underline{y} − \underline{μ}) − λ[(\underline{y} − \underline{μ})^{T}Σ̸^{−1}(\underline{y} − \underline{μ}) − c]  (2)

with respect to \underline{y} and setting these derivatives equal to zero leads to the eigenvalue/eigenvector problem

(Σ̸ − λI)(\underline{y} − \underline{μ}) = 0

So there must be multiple critical points for (2), one for one corresponding to each eigenvalue of σ. It now a matter of calculation to see that the eigenvalues will satisfy

λ_{1} ≥ λ_{2} ≥ ⋯ ≥ λ_{p} > 0

with corresponding eigenvectors \underline{e}_{i} = (\underline{y}_{i} − \underline{μ}) lying along the principal axes of the ellipsoid (\underline{y} − \underline{μ})^{T}Σ̸^{−1}(\underline{y} − \underline{μ}) = c; the half-lengths of the principal axes will be \sqrt{λc}.

Several interesting algebraic and geometric properties come out of this analysis. To every eigenvalue λ_{i} there corresponds an eigenvector \underline{e}_{i}, and it can be shown that all the \underline{e}_{i} are linearly independent. If two eigenvectors \underline{e}_{i} and \underline{e}_{j} correspond to different eigenvalues λ_{i} ≠ λ_{j} then, since \underline{e}_{j}^{T}Σ̸\underline{e}_{i} is a scalar (hence equal to its transpose),

λ_{i}\underline{e}_{j}^{T} \underline{e}_{i} = \underline{e}_{j}^{T}Σ̸\underline{e}_{i} = [\underline{e}_{j}^{T}Σ̸\underline{e}_{i}]^{T} = \underline{e}_{i}^{T}Σ̸\underline{e}_{j} = λ_{j}\underline{e}_{i}^{T} \underline{e}_{j}

so it must be that \underline{e}_{i}^{T} \underline{e}_{i} = 0: the eigenvectors from different eigenvalues are orthogonal. If λ_{i} = λ_{j}, there will still be two eigenvectors, and using the Gram-Schmidt procedure we can select two eigenvectors that are orthogonal to each other and to all the other eigenvectors in the problem. If the two eigenvalues are equal, then the cross-section of the ellipse in the plane containing the two linearly independent eigenvectors \underline{e}_{i} and \underline{y}_{j} will be circular (it would always be an ellipse, and now the major and minor principal axes will be of the same length 2\sqrt{λ_{i}c}); if three are equal the corresponding cross-section will be spherical, etc.

We can assume all the e_{i} are orthogonal; so let \underline{a}_{i} = \underline{e}_{i}/∥\underline{e}_{i}∥ and consider the matrix

A = [ \underline{a}_{1} \underline{a}_{2} ⋯ \underline{a}_{p} ]

You should verify that AA^{T} = I (the matrix is orthogonal—it has orthonormal columns); that Σ̸A = A diag (λ_{1}, …, λ_{p}) (the square p × p matrix with λ’s on the main diagonal and zeros everywhere else); and therefore that A^{T}Σ̸A = diag (λ_{1}, …, λ_{p}).

Now consider the random vector \underline{W} = A(\underline{Y} − \underline{μ}). This vector \underline{W} will be multivariate normal, and you should verify that it has mean \underline{0} (the zero vector) and variance-covariance matrix diag (λ_{1}, …, λ_{p}). Hence each W_{i} is uncorrelated with W_{j} (in fact, independent); and Var(W_{i}) = λ_{i}. The W_{i}’s are called the principal components of \underline{Y} and they play a large role in understanding the structure of \underline{Y}.

2 The Gauss-Markov Theorem

In the multiple regression model

\underline{y}_{n×1} = X_{n×k′} \underline{β}_{k′×1} + \underline{ϵ}_{n×1}

we will only assume that the error terms ϵ_{i} are independent, mean 0, constant variance σ_{ϵ}^{2}.

Let us define M = (X^{T} X)^{−1}X^{T}, so that ˆ \underline{β} = M\underline{y} are the least-squares estimators of \underline{β}, and let ˜\underline{β} = U\underline{y} be any unbiased linear estimate of β. Finally, let A = U − M.

Lemma 1 For any unbiased linear estimate U\underline{y} of \underline{β}, A = U − M is orthogonal to the rows of X; hence AM^{T} = MA^{T} = 0.

Proof: Since both ˜ \underline{β} and ˆ \underline{β} are unbiased,

\underline{β} = E[˜\underline{β}] = E[(M + A)\underline{y}] = E[M\underline{y}] + E[M\underline{y}] = E[ˆ\underline{β}] + AE[\underline{y}] = \underline{β} + AX\underline{β},

and since this must be true for all \underline{β} it follows that AX = 0. It immediately follows that AM^{T} = AX(X^{T} X)^{−1} = 0. □

Theorem 1 Gauss-Markov Theorem.

ˆ\underline{β} is the Best (minimum-variance) Linear Unbiased Estimate (BLUE) of \underline{β}.

Proof: We will do this by minimizing the variance of the arbitrary unbiased linear estimate ˜ \underline{β} = U\underline{y} = (M + A)\underline{y}. Since AM^{T} = MA^{T} = 0, it follows that

Var(˜\underline{β}) = Var[(M + A)\underline{y}]
= (M + A)Var(\underline{y})(M + A)^{T} = (M + A)σ_{ϵ}^{2}I(M + A)^{T}
= σ_{ϵ}^{2}(MM^{T} + MA^{T} + AM^{T} + AA^{T})
= σ_{ϵ}^{2}(MM^{T} + AA^{T})

and hence, reading off the diagonal elements of this variance-covariance matrix,

Var(β̃_{j}) = σ_{ϵ}^{2} (∑_{i=1}^{n} m_{ij}^{2} + ∑_{i=1}^{n} a_{ij}^{2}) .

Since the m_{ij}’s are fixed, we choose different unbiased linear estimates ˜ \underline{β} by choosing different a_{ij}’s; from the above expression it is clear that the only way to minimize the variance of all the β̃_{j}’s is to take a_{ij} ≡ 0; so that β̃_{j} = β̂_{j} for every j. □

Remark: Essentially the same technique can be used to prove the following result:

Corollary 1 Let φ = L\underline{β} be a linear function of the β’s. The BLUE of φ is Lˆ\underline{β} (i.e. plug in the least-squares estimates).