25.4. Multiple Regression#

Regression provides one way of predicting a numerical variable, called a response, based on other variables called predictor variables. The multiple regression model says in essence that

response = linear combination of predictor variables+ random noise 

You can think of the first term on the right hand side as a signal. The problem is that we don’t get to observe the signal. The observed response is the sum of the signal and the noise. The data scientist’s task is to use the observations to extract the signal as accurately as possible.

It is worth looking more closely at exactly what is linear in linear regression, now that we are allowing more than one predictor variable. For example, notice that you can fit a quadratic function of x by using the two predictor variables x1=x and x2=x2. Then the signal

β0+β1x1+β2x2 = β0+β1x+β2x2

is a quadratic function of x. But it is linear in the coefficients, and it is a linear combination of the two predictor variables x1 and x2.

25.4.1. The Model#

As in all of statistical inference, properties of estimates depend on the assumptions under which they are calculated. The multiple regression model is a commonly used set of assumptions that describes a particular kind of linear relation between a numerical response variable and a set of predictor variables. You should use it only if you believe that it makes sense for your data.

The model assumes that there are n individuals, on each of whom you have measured the response and the predictor variables. For 1in, the relation between the variables is assumed to be

Yi=β0+β1xi,1+β2xi,2++βp1xi,p1+ϵi

in the notation described below.

  • xi,1,xi,2,,xi,p1 are the observed constant values of p1 predictor variables for individual i. They are not random variables. If you prefer to think of the predictor variables as random, this model assumes that you have conditioned on them.

  • The intercept β0 and slopes β1,β2,,βp1 are unobservable constants and are parameters of the model. There are p of them, hence the notation p for “parameters”.

  • ϵi is an unobservable random error that has the normal (0,σ2) distribution for some unobservable σ2, and ϵ1,ϵ2,,ϵn are i.i.d.

  • Yi is the observable response of individual i. It is random because ϵi is one of its components.

We will assume that n>p, that is, we will assume we have more individuals than parameters. Indeed in this course it is fine for you to think of n as much larger than p.

Two special cases are already familiar.

p=1: Prediction by a Constant

When p=1 there is just one parameter: the intercept. There are no predictor variables at all. The model says that for each individual i, the response is Yi=β0+ϵi. This is a case of trying to estimate the response by a constant.

p=2: Simple Linear Regression

The two parameters are the intercept and a slope. The model says that for each individual i, the response is Yi=β0+β1xi,1+ϵi. That is, the response is the value on a hidden straight line, plus some normal noise. This is the simple regression model you used in Data 8.

25.4.2. Signal and Noise: Matrix Representation#

For any p, the model can be written compactly as

Y = Xβ+ϵ

in the matrix notation described below.

  • The design matrix X is an n×p matrix of real numbers, not random variables. Column 0 of X is a vector of 1’s and Column j for 1jp1 consists of the n observations on the jth predictor variable. For each i in the range 1 through n, Row i contains the values of all the predictor variables for individual i.

  • The parameter vector β=[β0  β1    βp1]T is a p×1 vector of the coefficients.

  • The error vector ϵ is an n×1 multivariate normal (0,σ2In) random vector. Its mean vector is an n×1 vector of 0’s and In is the n×n identity matrix.

  • The response vector Y is a random vector that is the sum of the linear signal Xβ and the normal noise ϵ.

25.4.3. Ordinary Least Squares#

Based on the observations of the predictor variables and the response, the goal is to find the best estimates of the intercept and slopes in the model.

These estimates can then be used to predict the response of a new individuals, assuming that the model holds for the new individual as well.

We must select a criterion by which we will decide whether one estimate is better than another. To develop one such criterion, start by noting that any linear function of the predictor variables can be written as Xγ where γ is some p×1 vector of coefficients. Think of Xγ as an estimate of Y. Then the error in the estimate is YXγ.

The goal of ordinary least squares (OLS) is to find the vector γ that minimises the mean squared error

MSE(γ) = 1ni=1n(Yi(Xγ)i)2

This is the same as the γ that minimizes the sum of squared errors

SSE(γ) = i=1n(Yi(Xγ)i)2

Again for compactness it will help to use matrix notation. For an n×1 vector w,

i=1nwi2 = wTw = ww = w2

which is sometimes called the squared norm of w.

In this notation, the goal of OLS is to find the p×1 vector β^ that minimizes YXγ2 over all vectors γ.

Typically you will also have to estimate the unknown error variance σ2. But we will not cover that in this class except in the case p=1.

See More

25.4.4. Guessing the Best Estimate of β#

Remember that we have assumed n>p. Assume also that X is of full column rank p, that is, none of the predictor variables is a linear combination of the others. By a theorem in linear algebra, it follows that the p×p matrix XTX has full rank p and is therefore invertible.

The claim is that OLS estimate of β is the vector β^ defined by

β^ = (XTX)1XTY

The claim is motivated by our earlier formula

b = ΣX1ΣXY

for the coefficients of the least squares linear predictor a random variable Y based on a random vector X. In fact the new formula is an application of the old one. But we will derive it afresh in our new setting.

The key idea is that of projection: the best β^ should be such that the error in the estimate is orthogonal to the linear space spanned by X.

The error in the best estimate is YY^=YXβ^. For this error to be orthogonal to linear transformations of X we must have

XT(YXβ^)=0 XTY = XTXβ^

We have assumed that X has full column rank, so XTX is invertible. So the natural guess for the best estimator β^ is

β^ = (XTX)1XTY

Before we go further, notice that β^ is a linear function of Y. This makes it straightforward to identify its distribution, which you will do in exercises.

Also note that the estimate of Y is

Y^ = Xβ^ = X(XTX)1XTY

which is also a linear function of Y.

25.4.5. Projection#

Define the ith residual as the prediction error ei=YiY^i. Then the n×1 vector of residuals is

e = YY^ = YXβ^

As we have seen repeatedly, the key to least squares is that the prediction error is orthogonal to the space of allowed functions. Our space of allowed functions is all linear functions of X. So we will show:

The residual vector is orthogonal to each column of X.

This is essentially true by construction. Formally, calculate the p×1 vector XTe. Each of its elements is the dot product of e and one column of X. We will show that each of the elements is 0.

XTe = XT(YY^) = XTYXTY^ = XTYXTX(XTX)1XTY = XTYXTY = 0

25.4.6. Least Squares#

Let γ be any p×1 vector. Then

SSE(γ) = YXγ2= (YXβ^)+(Xβ^Xγ)2= YXβ^2 + Xβ^Xγ2+2(Xβ^Xγ)T(YXβ^)= SSE(β^) + Xβ^Xγ2+2((X(β^γ))Te= SSE(β^) + Xβ^Xγ2+2(β^γ)TXTe= SSE(β^) + Xβ^Xγ2     by orthogonality SSE(β^)
See More

25.4.7. Signal and Noise, Revisited#

Our regression model is

Y = Xβ+ϵ

Here

  • Xβ is the unobservable but non-random true signal

  • ϵ is an unobservable random vector consisting of the deviations of Y from the true plane Xβ. Elements of ϵ are mutually independent.

Once we have carried out the regression, our estimate of the response vector Y is the vector Y^=Xβ^.

The residual vector is

e = YY^ = YXβ^

Therefore we have another expression for the response vector Y. This expression is our best attempt at separating the signal from the noise.

Y = Xβ^+e

It is important to note the distinction between this identity and the model.

  • Xβ^ is the observable random estimated signal.

  • e is the observable random vector consisting of the deviations of Y from the estimated plane Xβ^. Elements of e are not independent of each other, because they add up to 0.

In exercises you will show that β^ is an unbiased estimator of β and that both β^ and Y^ have normal (or multivariate normal) distributions. Both distributions have variance and covariance parameters that depend on the unknown error variance σ2.

25.4.8. Estimate of σ2#

It should come as no surprise that under the multiple regression model, there is an unbiased estimator of σ2 that has a chi-squared distribution. There is some work involved in establishing that this estimate is

S2 = 1npi=1n(YiY^i)2 = 1npe2

Some more work establishes that npσ2S2 has the chi-squared (np) distribution.

We’ll leave that work for another course. For now, just notice that if the number of data points n is large compared to the number of parameters p, then

S2 = 1npi=1n(YiY^i)2  1ni=1n(YiY^i)2 = σ^2

which is the natural mean squared error. If you have a lot of data, you don’t have to worry about fine points like dividing by np instead of n.

Special Case

As noted earlier, in the case p=1 you are trying to find the best constant by which to estimate Y.

You know that the least squares constant is Y¯, and you showed in exercises that

S2 = 1n1i=1n(YiY¯)2

is an unbiased estimate of σ2. You have shown in exercises that n1σ2S2 has the chi-squared n1 distribution under the assumption that the data are i.i.d. normal variables. This is the special case of the result stated above for general p.

25.4.9. Confidence Intervals#

The upshot of this discussion is that if n is large compared to p – that is, if you have a lot of observations compared to the number of predictors – then you can use ordinary normal theory to construct confidence intervals for the parameters β^.

For example, a 95% confidence interval for the parameter βi is β^[i]±2SD(β^[i]).

Here β^[i] is the ith element of the estimate vector β^.

The variance of β^[i] is the ith diagonal element of the covariance matrix of β^. You will see in exercises that this involves the error variance σ2. Typically, σ2 is unknown. But you can estimate it by the mean squared error σ^2 to get an approximate 95% confidence interval for βi.