Model Diagnostics

Note

The following supplemental notes were created by Dr. Maria Tackett for STA 210. They are provided for students who want to dive deeper into the mathematics behind regression and reflect some of the material covered in STA 211: Mathematics of Regression. Additional supplemental notes will be added throughout the semester.

This document discusses some of the mathematical details of the model diagnostics - leverage, standardized residuals, and Cook’s distance. We assume the reader knowledge of the matrix form for multiple linear regression. Please see Matrix Form of Linear Regression for a review.

Introduction

Suppose we have nn observations. Let the ithi^{th} be (xi1,…,xip,yi)(x_{i1}, \ldots, x_{ip}, y_i), such that xi1,…,xipx_{i1}, \ldots, x_{ip} are the explanatory variables (predictors) and yiy_i is the response variable. We assume the data can be modeled using the least-squares regression model, such that the mean response for a given combination of explanatory variables follows the form in Equation 1.

y=β0+β1x1+…+βpxp(1) y = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p \qquad(1)

We can write the response for the ithi^{th} observation as shown in Equation 2.

yi=β0+β1xi1+…+βpxip+ϵi(2) y_i = \beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip} + \epsilon_i \qquad(2)

such that ϵi\epsilon_i is the amount yiy_i deviates from μ{y|xi1,…,xip}\mu\{y|x_{i1}, \ldots, x_{ip}\}, the mean response for a given combination of explanatory variables. We assume each ϵi∼N(0,σ2)\epsilon_i \sim N(0,\sigma^2), where σ2\sigma^2 is a constant variance for the distribution of the response yy for any combination of explanatory variables x1,…,xpx_1, \ldots, x_p.

Matrix Form for the Regression Model

We can represent the Equation 1 and Equation 2 using matrix notation. Let

𝐘=[y1y2⋮yn]𝐗=[x11x12…x1px21x22…x2p⋮⋮⋱⋮xn1xn2…xnp]𝛃=[β0β1⋮βp]𝛜=[ϵ1ϵ2⋮ϵn](3) \mathbf{Y} = \begin{bmatrix}y_1 \\ y_2 \\ \vdots \\y_n\end{bmatrix} \hspace{1cm} \mathbf{X} = \begin{bmatrix}x_{11} & x_{12} & \dots & x_{1p} \\ x_{21} & x_{22} & \dots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots\\ x_{n1} & x_{n2} & \dots & x_{np} \end{bmatrix} \hspace{1cm} \boldsymbol{\beta}= \begin{bmatrix}\beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{bmatrix} \hspace{1cm} \boldsymbol{\epsilon}= \begin{bmatrix}\epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{bmatrix} \qquad(3)

Thus,

𝐘=𝐗𝛃+𝛜\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{\epsilon}

Therefore the estimated response for a given combination of explanatory variables and the associated residuals can be written as

𝐘̂=𝐗𝛃̂𝐞=𝐘−𝐗𝛃̂(4) \hat{\mathbf{Y}} = \mathbf{X}\hat{\boldsymbol{\beta}} \hspace{1cm} \mathbf{e} = \mathbf{Y} - \mathbf{X}\hat{\boldsymbol{\beta}} \qquad(4)

Hat Matrix & Leverage

Recall from the notes Matrix Form of Linear Regression that 𝛃̂\hat{\boldsymbol{\beta}} can be written as the following:

𝛃̂=(𝐗T𝐗)−1𝐗T𝐘(5) \hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y} \qquad(5)

Combining Equation 4 and Equation 5, we can write 𝐘̂\hat{\mathbf{Y}} as the following:

𝐘̂=𝐗𝛃̂=𝐗(𝐗T𝐗)−1𝐗T𝐘(6) \begin{aligned} \hat{\mathbf{Y}} &= \mathbf{X}\hat{\boldsymbol{\beta}} \\[10pt] &= \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y}\\ \end{aligned} \qquad(6)

We define the hat matrix as an n×nn \times n matrix of the form 𝐇=𝐗(𝐗T𝐗)−1𝐗T\mathbf{H} = \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T. Thus Equation 6 becomes

𝐘̂=𝐇𝐘(7) \hat{\mathbf{Y}} = \mathbf{H}\mathbf{Y} \qquad(7)

The diagonal elements of the hat matrix are a measure of how far the predictor variables of each observation are from the means of the predictor variables. For example, hiih_{ii} is a measure of how far the values of the predictor variables for the ithi^{th} observation, xi1,xi2,…,xipx_{i1}, x_{i2}, \ldots, x_{ip}, are from the mean values of the predictor variables, x‾1,x‾2,…,x‾p\bar{x}_1, \bar{x}_2, \ldots, \bar{x}_p. In the case of simple linear regression, the ithi^{th} diagonal, hiih_{ii}, can be written as

hii=1n+(xi−x‾)2∑j=1n(xj−x‾)2 h_{ii} = \frac{1}{n} + \frac{(x_i - \bar{x})^2}{\sum_{j=1}^{n}(x_j-\bar{x})^2}

We call these diagonal elements, the leverage of each observation.

The diagonal elements of the hat matrix have the following properties:

  • 0≤hii≤10 \leq h_ii \leq 1
  • ∑i=1nhii=p+1\sum\limits_{i=1}^{n} h_{ii} = p+1, where pp is the number of predictor variables in the model.
  • The mean hat value is h‾=∑i=1nhiin=p+1n\bar{h} = \frac{\sum\limits_{i=1}^{n} h_{ii}}{n} = \frac{p+1}{n}.

Using these properties, we consider a point to have high leverage if it has a leverage value that is more than 2 times the average. In other words, observations with leverage greater than 2(p+1)n\frac{2(p+1)}{n} are considered to be high leverage points, i.e. outliers in the predictor variables. We are interested in flagging high leverage points, because they may have an influence on the regression coefficients.

When there are high leverage points in the data, the regression line will tend towards those points; therefore, one property of high leverage points is that they tend to have small residuals. We will show this by rewriting the residuals from Equation 4 using Equation 7.

𝐞=𝐘−𝐘̂=𝐘−𝐇𝐘=(1−𝐇)𝐘(8) \begin{aligned} \mathbf{e} &= \mathbf{Y} - \hat{\mathbf{Y}} \\[10pt] & = \mathbf{Y} - \mathbf{H}\mathbf{Y} \\[10pt] &= (1-\mathbf{H})\mathbf{Y} \end{aligned} \qquad(8)

Note that the identity matrix and hat matrix are idempotent, i.e. 𝐈𝐈=𝐈\mathbf{I}\mathbf{I} = \mathbf{I}, 𝐇𝐇=𝐇\mathbf{H}\mathbf{H} = \mathbf{H}. Thus, (𝐈−𝐇)(\mathbf{I} - \mathbf{H}) is also idempotent. These matrices are also symmetric. Using these properties and Equation 8, we have that the variance-covariance matrix of the residuals 𝐞\boldsymbol{e}, is

Var(𝐞)=𝐞𝐞T=(1−𝐇)Var(𝐘)T(1−𝐇)T=(1−𝐇)σ̂2(1−𝐇)T=σ̂2(1−𝐇)(1−𝐇)=σ̂2(1−𝐇)(9) \begin{aligned} Var(\mathbf{e}) &= \mathbf{e}\mathbf{e}^T \\[10pt] &= (1-\mathbf{H})Var(\mathbf{Y})^T(1-\mathbf{H})^T \\[10pt] &= (1-\mathbf{H})\hat{\sigma}^2(1-\mathbf{H})^T \\[10pt] &= \hat{\sigma}^2(1-\mathbf{H})(1-\mathbf{H}) \\[10pt] &= \hat{\sigma}^2(1-\mathbf{H}) \end{aligned} \qquad(9)

where σ̂2=∑i=1nei2n−p−1\hat{\sigma}^2 = \frac{\sum_{i=1}^{n}e_i^2}{n-p-1} is the estimated regression variance. Thus, the variance of the ithi^{th} residual is Var(ei)=σ̂2(1−hii)Var(e_i) = \hat{\sigma}^2(1-h_{ii}). Therefore, the higher the leverage, the smaller the variance of the residual. Because the expected value of the residuals is 0, we conclude that points with high leverage tend to have smaller residuals than points with lower leverage.

Standardized Residuals

In general, we standardize a value by shifting by the expected value and rescaling by the standard deviation (or standard error). Thus, the ithi^{th} standardized residual takes the form

std.resi=ei−E(ei)SE(ei) std.res_i = \frac{e_i - E(e_i)}{SE(e_i)}

The expected value of the residuals is 0, i.e. E(ei)=0E(e_i) = 0. From Equation 9), the standard error of the residual is SE(ei)=σ̂1−hiiSE(e_i) = \hat{\sigma}\sqrt{1-h_{ii}}. Therefore,

std.resi=eiσ̂1−hii(10) std.res_i = \frac{e_i}{\hat{\sigma}\sqrt{1-h_{ii}}} \qquad(10)

Cook’s Distance

Cook’s distance is a measure of how much each observation influences the model coefficients, and thus the predicted values. The Cook’s distance for the ithi^{th} observation can be written as

Di=(𝐘̂−𝐘̂(i))T(𝐘̂−𝐘̂(i))(p+1)σ̂(11) D_i = \frac{(\hat{\mathbf{Y}} -\hat{\mathbf{Y}}_{(i)})^T(\hat{\mathbf{Y}} -\hat{\mathbf{Y}}_{(i)})}{(p+1)\hat{\sigma}} \qquad(11)

where 𝐘̂(i)\hat{\mathbf{Y}}_{(i)} is the vector of predicted values from the model fitted when the ithi^{th} observation is deleted. Cook’s Distance can be calculated without deleting observations one at a time, since Equation 12 below is mathematically equivalent to Equation 11.

Di=1p+1std.resi2[hii(1−hii)]=ei2(p+1)σ̂2(1−hii)[hii(1−hii)](12) D_i = \frac{1}{p+1}std.res_i^2\Bigg[\frac{h_{ii}}{(1-h_{ii})}\Bigg] = \frac{e_i^2}{(p+1)\hat{\sigma}^2(1-h_{ii})}\Bigg[\frac{h_{ii}}{(1-h_{ii})}\Bigg] \qquad(12)