A More Careful Development of the Linear Regression Model 1 Estimating the Coefficients in the Multiple Regression Model

The multiple regression model that we have been working with may be written as

y = β_{0} + β_{1}x_{1} + β_{2}x_{2} + ⋯ + β_{k}x_{k} + ϵ

where ϵ ∼ N(0, σ^{2}). As usual, this stands for

y_{1} = β_{0} + β_{1}x_{11} + β_{2}x_{12} + ⋯ + β_{k}x_{1k} + ϵ_{1}
y_{2} = β_{0} + β_{1}x_{21} + β_{2}x_{22} + ⋯ + β_{k}x_{2k} + ϵ_{2}
⋮
y_{n} = β_{0} + β_{1}x_{n1} + β_{2}x_{n2} + ⋯ + β_{k}x_{nk} + ϵ_{n}

or, in matrix notation,

\underline{y} = X\underline{β} + \underline{ϵ}, \underline{ϵ} ∼ N(\underline{0}, σ^{2}I_{n×n})

i.e.

[ y_{1} ; y_{2} ; ⋮ ; y_{n} ] = [ \begin{matrix}1 & x_{11} & x_{12} & ⋯ & x_{1k} \\ 1 & x_{21} & x_{22} & ⋯ & x_{2k} \\ & ⋮ & ⋮ & ... & ⋮ \\ 1 & x_{n1} & x_{n2} & ⋯ & x_{nk} \end{matrix} ] [ β_{0} ; β_{1} ; β_{2} ; ⋮ ; β_{k} ] + [ ϵ_{1} ; ϵ_{2} ; ⋮ ; ϵ_{n} ]

We know that

Mean: E[y_{i}] = E[β_{0} + β_{1}x_{i1} + … β_{k}x_{ik} + ϵ_{i}]
= β_{0} + β_{1}x_{i1} + … β_{k}x_{ik} + E[ϵ_{i}]
= β_{0} + β_{1}x_{i1} + … β_{k}x_{ik} + 0
= β_{0} + β_{1}x_{i1} + … β_{k}x_{ik}, and
Variance: Var(y_{i}) = Var(β_{0} + β_{1}x_{i1} + … β_{k}x_{ik} + ϵ_{i})
= Var(ϵ_{i})
= σ^{2}.

which we write compactly in matrix notation as

E[\underline{y}] = X\underline{β}
Var(\underline{y}) = σ^{2}I_{n×n}, and so
\underline{Y} ∼ N(X\underline{β}, σ^{2}I_{n×n})

If we have some guesses

ˆ\underline{β} = (β̂_{0}, β̂_{1}, …, β̂_{k})^{T}

for \underline{β}, then the fitted values are ˆ\underline{y} = X ˆ\underline{β} and the residuals are

\underline{E} = \underline{y} − ˆ \underline{y}

To find good estimates of \underline{β}, we may minimize the sum of squares of residuals,

Q(β̂_{0}, β̂_{1}, …, β̂_{k}) = ∑_{i=1}^{n} E_{i}^{2} = ∑_{i=1}^{n}(y_{i} − [β̂_{0} + β̂_{1}x_{i1} + … β̂_{k}x_{ik}])^{2}

= ∑_{i=1}^{n}(y_{i} − ŷ_{i})^{2}
= (\underline{y} − ˆ \underline{y})^{T} (\underline{y} − ˆ \underline{y})
= (\underline{y} − X ˆ\underline{β})^{T} (\underline{y} − X ˆ\underline{β})

One way to proceed is to take derivatives (\frac{∂}{∂\underline{x}}A\underline{x} = \underline{x}, and \frac{∂}{∂\underline{x}} \underline{x}^{T} A\underline{x} = 2A\underline{x}) and equate to zero:

\frac{∂}{∂ˆ\underline{β}} [(\underline{y} − X ˆ\underline{β})^{T} (\underline{y} − X ˆ\underline{β})]
= \frac{∂}{∂ˆ\underline{β}} (\underline{y} − X ˆ\underline{β})^{T} · 2(\underline{y} − X ˆ\underline{β})
= −2X^{T} (\underline{y} − X ˆ\underline{β})

Setting this last expression equal to zero produces the normal equations,

X^{T} (\underline{y} − X ˆ\underline{β}) = X^{T} E = \underline{0}

because this asserts that the columns of X are orthogonal to, or in older mathematical language “normal to”, the residuals E = \underline{y} − X ˆ\underline{β}. Distributing X^{T} across the (\underline{y} − X ˆ\underline{β}) term we obtain

X^{T} \underline{y} = X^{T} X ˆ\underline{β}

and—assuming that X^{T} X is an invertible matrix—we can premultiply this equation by (X^{T} X)^{−1} to obtain

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

Note that it is just as easy to obtain the fitted values ˆ \underline{y}:

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

The matrix H = X(X^{T} X)^{−1}X^{T} is often called the hat matrix for linear regression, because the effect of multiplying y by H is to “put the hat on”:

ˆ\underline{y} = Hy

Remarks:

• The fact that the fitted values ˆ \underline{y} are a \underline{linear transformation of the original responses} \underline{y} is \underline{really} what makes linear regression “linear”. Nonparametric regression methods such as

– Kernel smoothing – Spline smoothing

– \underline{but not loess}

can be shown to be linear smoothers in exactly this sense. Linearity makes obtaining distribution theory, large sample theory, etc., particularly easy.

• Note that

\frac{∂}{∂\underline{y}} ˆ \underline{y} = H

and in particular

\frac{∂}{∂y_{i}} ŷ_{i} = h_{ii},

the i^{th} diagonal element of the hat matrix H. Thus, the larger h_{ii} is, the more sensitive ŷ_{i} is to changes in y_{i}; this further explains why h_{ii} is a good place to start for influence diagnostics.

• We can also obtain the residuals using H:

\underline{E} = \underline{y} − ˆ \underline{y} = \underline{y} − H\underline{y} = (I − H)\underline{y}

• H has many remarkable properties. Let

P = H
Q = 1 − H

Then it is not hard to show, using the definition of H and/or the normal equations, that

P^{T} = P, Q^{T} = Q : symmetric
P^{2} = P, Q^{2} = Q : idempotent
P Q = 0 :  orthogonal
P + Q = I :  complete

and indeed P and Q are “projection matrices” in n-dimensional real space: every vector x in ℜ^{n} can be written uniquely as an orthogonal decomposition

x = x_{P} + x_{Q}

where x_{P} = P x ∈ ColSpace(P), x_{Q} = Qx ∈ ColSpace(Q), and x_{P} ⊥ x_{Q} (they are orthogonal).

• Calculating the projection and properties of the projection matrix are important, but the algebraic formula H = X(X^{T} X)^{−1}X^{T} involving matrix inversion need not be the best computational formula. See Thisted’s Statistical Computing text for linear regression estimation formulae that are actually useful computationally.

2 Some Geometry

Consider a dependent variable \underline{y} = (y_{1}, …, y_{n})^{T} and k = 2 independent variables,

\underline{X_{1}} = ( x_{11} ; x_{21} ; ⋮ ; x_{n1} ) \underline{X_{2}} = ( x_{12} ; x_{22} ; ⋮ ; x_{n2} )

We wish to consider the model y_{i} = β_{0} + β_{1}x_{i1} + β_{2}x_{i2} + e_{i}. Letting

X = ( \begin{matrix}1 & x_{11} & x_{12} \\ 1 & x_{21} & x_{22} \\ ⋮ & ⋮ & ⋮ \\ 1 & x_{n1} & x_{n2} \end{matrix} ) = (\underline{1}, \underline{X_{1}}, \underline{X_{2}}),

\underline{ϵ} = (ϵ_{1}, …, ϵ_{n})^{T}, and \underline{β} = (β_{1}, β_{2}, β_{3})^{T}, we may again write this in matrix form as

\underline{y} = X\underline{β} + \underline{ϵ}

The least-squares fitted values ˆ \underline{y} minimize the squared distance ∥\underline{y} − ˆ \underline{y}∥^{2}, and so, looking at the picture and thinking about elementary geometry, it’s clear that the vector \underline{E} = \underline{y} − ˆ \underline{y} should be perpendicuar to M = Span(\underline{1}, \underline{X_{1}}, \underline{X_{2}}) = ColSpace(X), the vector space spanned by the columns of X.

On the other hand, if X^{T} X is of full rank, then the algebra we did above asserts that the least squares estimates satisfy

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

and the fitted or predicted values for the y’s satisfy

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

This equation shows immediately that ˆ \underline{y} ∈ ColSpace(X), and since \underline{E} = (I − H)\underline{y} we also know that, for any \underline{β} (least-squares estimates or not)

(X\underline{β})^{T} \underline{E} = \underline{β}^{T} X^{T} (I − H)\underline{y} = \underline{β}^{T} (X^{T} − X^{T} X(X^{T} X)^{−1}X^{T}) = \underline{β}^{T} (X^{T} − X^{T}) = \underline{0}

so that the residuals \underline{E} = \underline{y} − ˆ \underline{y} really are orthogonal to every possible linear combination of the columns of X.

3 Assumptions and Standard Error Calculations

All of the preceding was just algebra; lets do some statistics. The model is

\underline{y} = X\underline{β} + \underline{ϵ}

and we assume that 1. X is of full rank (this is easily fixed)

2. E(\underline{ϵ}) = \underline{0}

3. Var(\underline{ϵ}) = σ^{2}\underline{I}_{n×n} which incorporates the equal variance and zero covariance assumptions

4. \underline{ϵ} ∼ N(\underline{0}, σ^{2}\underline{I}) a multivariate normal distribution.

3.1 Digression — Multivariate Normal Distribution

Definition

A random vector

\underline{x} = ( x_{1} ; x_{2} ; ⋮ ; x_{n} )

is said to have a multivariate normal distribution N(\underline{μ}, Σ) if

\underline{μ} = ( μ_{1} ; μ_{2} ; ⋮ ; μ_{n} )

Σ = ( \begin{matrix}σ_{11} & σ_{12} & … & σ_{1n} \\ σ_{21} & σ_{22} & … & σ_{2n} \\ ⋮ & ... & ... & ⋮ \\ σ_{n1} & σ_{n2} & … & σ_{nn} \end{matrix} ), and

E(\underline{x}) = \underline{μ}, and
Cov(x_{i}, x_{j}) = σ_{i,j}

AND the distribution of

\underline{a}^{T} \underline{x} is N(\underline{a}^{T} \underline{μ}, \underline{a}^{T} Σ\underline{a})

for all vectors \underline{a}. In other words it is multivariate normal iff every possible linear combination is univariate normal.

Transformation Properties

If \underline{x} ∼ N_{n}(\underline{μ}, Σ) and W_{p×n} is a matrix of constants then

W \underline{x} ∼ N_{p}(W \underline{μ}, W ΣW^{T})

Now let Σ_{n×n} be decomposed as V_{n×n}V_{n×n}^{T}, where each V is of full rank (called a square root of W. This can always be done (in many ways). Then

V^{−1}\underline{x} ∼ N_{n}(V^{−1}\underline{μ}, \underbrace{V^{−1}(V V^{T})V^{−T}}_{I_{n×n}})

Thus, if \underline{x} is a random vector with mean μ and variance Σ and Σ = V V^{T}, then V^{−1}\underline{x} has mean V^{−1}\underline{μ} and variance I.

3.2 Sampling Distributions of ˆ \underline{β} and ˆ \underline{y}

We are now in the position to say something about ˆ \underline{β} = (X^{T} X)^{−1}X^{T} \underline{y} and ˆ \underline{y} = H\underline{y}.

For starters, note that

E(\underline{y}) = X\underline{β}, and
Var(\underline{y}) = σ^{2} · I.

Thus

\begin{aligned}E(ˆ\underline{β}) & = & (X^{T} X)^{−1}X^{T} E(\underline{y}) \\ & = & (X^{T} X)^{−1}X^{T} X\underline{β} & = & \underline{β} \end{aligned}

and

Var(ˆ\underline{β}) = (X^{−T} X)^{−1}X^{T} Var(\underline{y})X(X^{−T} X)^{−1}
= (X^{−T} X)^{−1}X^{T} (σ^{2} · I)X(X^{−T} X)^{−1}
= σ^{2}(X^{T} X)^{−1}

The fitted values ˆ \underline{y} and residuals \underline{E} are even easier

E[ˆ\underline{y}] = E[H\underline{y}] = HX\underline{β} = X\underline{β} = E[\underline{y}]
E[ ˆ\underline{E}] = E[(I − H)\underline{y}] = E[\underline{y}] − E[H\underline{y}] = E[\underline{y}] − E[\underline{y}] = \underline{0}

and

Var(ˆ\underline{y}) = Var(H\underline{y}) = HVar(\underline{y})H^{T} = Hσ^{2}IH^{T} = σ^{2}HH^{T} = σ^{2}H
Var(\underline{E}) = Var((I − H)\underline{y}) = ⋯ = σ^{2}(I − H)

Together with the assumption that \underline{ϵ} ∼ N(\underline{0}, σ^{2}I), or \underline{y} ∼ N(X\underline{β}, σ^{2}I), these facts show that

ˆ\underline{β} ∼ N(\underline{β}, σ^{2}(X^{T} X)^{−1}) ˆ \underline{y} ∼ N(X\underline{β}, σ^{2}H) \underline{E} ∼ N(\underline{0}, σ^{2}(I − H))

How do we estimate σ^{2}? From our earlier work we expect that we can get an unbiased estimate as a multiple of RSS. To calculate E[RSS] we need a few preliminaries.

First we need a fact from linear algebra, namely that if both AB and BA are legal matrix multiplications, then

trace(AB) = trace(BA)

where trace (A) = ∑_{i} a_{ii} is the sum of the diagonal elements of A. We immediately put this to good use, to see that

trace(H) = trace(X(X^{T} X)^{−1}X^{T}) = trace(X^{T} X(X^{T} X)^{−1}) = trace(I_{(k+1)×(k+1)}) = k + 1

Next we need a result about the expected value of a quadratic form

E[\underline{y}^{T} A\underline{y}] = trace(AVar(\underline{y})) + E[\underline{y}]^{T} AE[\underline{y}]

(this is analogous to the univariate fact that E[cy^{2}] = c^{2}σ^{2} + c^{2}E[y]^{2}). Finally we note that

RSS = (\underline{y} − ˆ \underline{y})^{T} (\underline{y} − ˆ \underline{y})
= \underline{y}^{T} [I − H]^{T} [I − H]\underline{y}
= \underline{y}^{T} [I − H]\underline{y}

With these three bits of information we calculate

E(RSS) = E[\underline{y}^{T} [I − H]\underline{y}]
= trace([I − H]σ^{2}I) + (Xβ)^{T} (I − H)(Xβ)
= σ^{2}trace(I − H) + \underline{0}
= σ^{2}[n − (k + 1)]

Consequently,

E [\frac{RSS}{n − (k + 1)}] = σ^{2}

so that

σ̂^{2} = \frac{RSS}{n − (k + 1)}

is an unbiased estimate of σ^{2}, as we expected.

3.3 Variances of Fitted Values and Predicted Values

We have done part of this already. Note that

ŷ_{i} = \underline{X_{i}}^{T} · ˆ \underline{β}
= \underline{X_{i}}^{T} (X^{T} X)^{−1}X^{T} \underline{y}

so that

Var(ŷ_{i}) = σ^{2}\underline{X_{i}}^{T} (X^{T} X)^{−1}\underline{X_{i}}

Notice that the variance of the fitted value depends on the vector X where y_{i} was observed.

At a new value \underline{x}^{∗} that was not in the original data set, the fitted value is still

ŷ^{∗} = \underline{x^{∗}}^{T} ˆ \underline{β}

with variance

Var(ŷ^{∗}) = σ^{2}\underline{x^{∗}}^{T} (X^{T} X)^{−1}\underline{x^{∗}}

Note that \underline{x^{∗}} does not have to be part of the original X matrix. The variance calculated here is for estimating E[y^{∗}] in the equation

y^{∗} = \underline{x^{∗}}^{T} β + ϵ^{∗}

Sometimes it is desirable to try to predict the actual value y^{∗} itself, and not just the expected value E[y^{∗}]. Assuming as usual that ϵ^{∗} is independent of the other ϵ’s in the problem, and comes from the same N(0, σ^{2}) distribution, we have simply

Var(ŷ_{pred}^{∗}) = Var(ŷ^{∗}) + Var(ϵ^{∗})
= σ^{2}\underline{x^{∗}}^{T} (X^{T} X)^{−1}\underline{x^{∗}} + σ^{2}
= σ^{2} [1 + \underline{x^{∗}}^{T} (X^{T} X)^{−1}\underline{x^{∗}}]

4 Some final remarks

Aside from being a clean way to derive and generalize all of the standard formulas from simple regression, there are several reasons for bothering with the linear algebra formulation:

• The geometry suggested above is sometimes quite helpful in figuring out what vague notions in regression mean. For example, “collinearity” of \underline{X_{1}} and \underline{X_{2}} simply means that the angle between \underline{X_{1}} and \underline{X_{2}} is so small that the position of the hyperplane Span (\underline{1}, \underline{X_{1}}, \underline{X_{2}}) is very sensitive to small changes in \underline{X_{1}} or \underline{X_{2}},

• Standard error calculations for the β’s are clearly straightforward (though messy calculations are hidden in the matrix multiplications!). Since we also have expressions for the covariances, we can easily obtain standard errors for any linear combination (contrast!) of the parameters.

• The geometry is also useful for suggesting what is the “same” model—any parametrization that generates the same H matrix is the “same” model as far as linear regression is concerned.

• Finally the geometry, mixed with a little distribution theory, will be very helpful in creating and understanding standard errors that are corrected for multiple inference, such as in Scheffé’s method.

You may also be interested in knowing that, just as in the simple regression case, the least-squares estimates are the same as the maximum likelihood estimates of the β’s; and as before the MLE of σ^{2} is the biased estimate RSS/n.