Chapter 3 Topics in Regression

3.1 Regression Through the Origin

Applications occasionally arise in which it is of benefit to force a regression line to pass through the origin. To illustrate such applications, return to Examples 1.1 and 1.3 in which Bob and Cheryl each had the number of sales per week as an independent variable X. In both of the examples, [latex]X = 0[/latex] sales per week corresponds to [latex]Y = 0[/latex] commissions (for Bob) and [latex]Y = 0[/latex] revenue per week (from Cheryl’s sales). In these settings it is sensible to force the regression line to pass through the origin; estimating a population intercept does not make sense. The resulting regression model does not contain the β₀ parameter. The simple linear regression model forced through the origin, sometimes abbreviated RTO for regression through the origin, is defined next.

The regression parameter β₁ can be estimated using least squares from the data pairs [latex](X_i, \, Y_i)[/latex] for [latex]i = 1, \, 2, \, \ldots, \, n[/latex].

The next example conducts a hypothesis test to determine whether it is appropriate to drop the intercept term from the simple linear regression model based on the data pairs, and then proceeds to fit the reduced model.

The next example revisits the regression modeling of the stopping distance as a function of the speed of a car in the built-in cars data frame.

This ends the discussion of forcing the regression line through the origin. Occasions arise in regression modeling in which it is more appropriate to fit a statistical model with fewer parameters. Some of the results from the full simple linear regression model generalize to simple linear regression forced through the origin. The point estimate for β₁, for example, is unbiased. Three examples of results that do not generalize are (a) the residuals do not necessarily sum to zero, (b) the regression line does not necessarily pass through the point [latex]\big( \bar X, \, \bar Y \big)[/latex], and (c) it is possible that SSE can exceed the total sum of squares SST, which can result in a negative value of R².

3.2 Diagnostics

Diagnostic procedures are applied to fitted regression models to assess their conformity to the assumptions (for example, constant variance of the error terms) implicit in the simple linear regression model. We have already considered one such diagnostic procedure from the previous chapter, which is the examination of the residuals to assess their independence, constant variance, and normality. Two other diagnostic procedures will be examined here, which are the identification of data pairs known as leverage points and the identification of data pairs known as influential points. The subsequent section considers remedial procedures, which can be applied to a regression model that fails to satisfy one or more of the assumptions implicit in a regression model.

3.3 Remedial Procedures

The diagnostic procedures presented in the previous section are designed to identify assumptions associated with the simple linear regression model that are not satisfied for a particular set of n data pairs. But these diagnostic procedures do not suggest remedies when model assumptions are not satisfied. This section considers remedial procedures.

Reasons that simple linear regression model with normal error terms can fail to satisfy the assumptions given in Definition 2.1 include

• the regression function is not linear,
• the regression model has not included an important independent variable,
• the error terms have a variance that varies with X,
• the error terms are not independent,
• the error terms are not normally distributed,
• the scope of the regression model is too wide,
• the scope of the regression model is too narrow, and
• an influential point has an unusually strong effect on the regression line.

Two common approaches to handling a regression model which violates one or more of the assumptions are (a) formulate and fit a regression model with nonlinear terms, and (b) transform the X-values or the Y-values (or both) in a fashion so that the simple linear regression assumptions are satisfied. Regression models with nonlinear terms will be considered in a subsequent section in this chapter; transformations will be considered here. Transformations will be illustrated in a single (long) example.

The previous example took a trial-and-error approach to determining an appropriate transformation to apply to the raw data pairs in order to satisfy the assumptions implicit in a simple linear regression model with normal error terms. There are templates that can give a more systematic approach to determining these transformations.

There is a nice synergy between matrix algebra and regression, which will be presented in the next section.

3.4 Matrix Approach to Simple Linear Regression

So far, a purely algebraic approach has been taken to simple linear regression modeling. This section considers a matrix-based approach. There are (at least) four reasons to take this approach. First, the mathematical expressions are in many cases much more compact; summations from the algebraic approach are often equivalent to matrix multiplications. Second, matrix algebra can easily be implemented on a computer. Third, the matrix approach generalizes very easily to the multiple regression case in which there are several independent variables. Fourth, the matrix approach generalizes very easily to weighted least squares, which will be introduced in the next section.

We begin the matrix approach by defining certain critical matrices, which will be set in boldface. Let [latex]{\bf X}[/latex] be an [latex]n \times 2[/latex] matrix whose first column is all ones and whose second column contains the observed values of the independent variable, [latex]{\bf Y}[/latex] be an [latex]n \times 1[/latex] vector which holds the observed values of the dependent variable, [latex]\pmb{\beta}[/latex] be a [latex]2 \times 1[/latex] vector which holds the population intercept and slope, and [latex]\pmb{\epsilon}[/latex] be an [latex]n \times 1[/latex] vector which holds the error terms:

$$\begin{array}{l}\mathbf{X}=\left[\begin{array}{cc}1& X_1\\ 1& X_2\\ ⋮& ⋮\\ 1& X_n\end{array}\right],\mathbf{Y}=\left[\begin{array}{c}{Y}_1\\ Y_2\\ ⋮\\ Y_n\end{array}\right],\mathit{\beta }=\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\end{array}\right],\text{and}\mathit{ϵ}=\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ ⋮\\ {ϵ}_n\end{array}\right].\end{array}$$

The [latex]\bf X[/latex] matrix is known as the design matrix.

As before, the values of the independent variable (the second column of [latex]\bf X[/latex]) are assumed to be fixed constants observed without error with at least two distinct values, the values of the dependent variable contained in [latex]\bf Y[/latex] are assumed to be continuous random responses, and the elements of the vector [latex]\pmb\epsilon[/latex] are assumed to be mutually independent random variables, each with population mean 0 and finite positive population variance σ². Stated another way, the expected value of [latex]\pmb\epsilon[/latex] is the zero vector and the variance–covariance matrix of [latex]\pmb\epsilon[/latex] is

$$\begin{array}{l}\left[\begin{array}{cccc}{\sigma }^2& 0& \cdots & 0\\ 0& {\sigma }^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\sigma }^2\end{array}\right].\end{array}$$

The simple linear regression model

$$\begin{array}{l}{Y}_i={\beta }_0+{\beta }_1{X}_i+{ϵ}_i\end{array}$$

for [latex]i = 1, \, 2, \, \ldots, \, n[/latex], can be written more explicitly in terms of each observed data pair as

$$\begin{array}{ll}{Y}_1& ={\beta }_0+{\beta }_1{X}_1+{ϵ}_1\\ Y_2& ={\beta }_0+{\beta }_1{X}_2+{ϵ}_2\\ & ⋮\\ Y_n& ={\beta }_0+{\beta }_1{X}_n+{ϵ}_n\end{array}$$

which, in matrix form, is

$$\begin{array}{l}\left[\begin{array}{c}{Y}_1\\ Y_2\\ ⋮\\ Y_n\end{array}\right]=\left[\begin{array}{cc}1& X_1\\ 1& X_2\\ ⋮& ⋮\\ 1& X_n\end{array}\right]\cdot \left[\begin{array}{c}{\beta }_0\\ {\beta }_1\end{array}\right]+\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ ⋮\\ {ϵ}_n\end{array}\right]\end{array}$$

or simply

$$\begin{array}{l}\mathbf{Y}=\mathbf{X}\mathit{\beta }+\mathit{ϵ}.\end{array}$$

This explains why the artificial column of ones appears as the first column of the [latex]{\bf X}[/latex] matrix; it is to account for the intercept term. To force a regression line through the origin, simply omit the column of ones in the [latex]{\bf X}[/latex] matrix. Taking the expected value of both sides of this equation results in

$$\begin{array}{l}E[\mathbf{Y}]=\mathbf{X}\mathit{\beta }\end{array}$$

because [latex]E[ \epsilon_i ] = 0[/latex], for [latex]i = 1, \, 2, \, \ldots, \, n[/latex], (that is, [latex]E[ \pmb{\epsilon} ] = {\bf 0}[/latex]). The left-hand side of this equation, [latex]E[ {\bf Y} ][/latex], is an n-element column vector with elements [latex]E[Y_1], \, E[Y_2], \ldots , \, E[Y_n][/latex]. The sum of squares which is to be minimized to find the least squares estimators is

$$\begin{array}{l}S={(\mathbf{Y}-\mathbf{X}\mathit{\beta })}^{\mathrm{\prime <!--\; \prime \; -->}}(\mathbf{Y}-\mathbf{X}\mathit{\beta }).\end{array}$$

With this notation established, the algebraic results concerning the simple linear regression model can be restated more compactly in terms of these matrices. The results have already been proved, so there is no need to prove them again when stated in matrix form. The [latex]\, ^ \prime[/latex] superscript denotes transpose. It is a good exercise to perform the algebra necessary to see that the algebraic and matrix versions of these definitions and theorems match. The dimensions of the matrices should be checked for conformity.

• Definition 1.1. The simple linear regression model is

  $$\begin{array}{l}\mathbf{Y}=\mathbf{X}\mathit{\beta }+\mathit{ϵ},\end{array}$$

  where [latex]E[ \pmb{\epsilon} ] = {\bf 0}[/latex], [latex]V[ \pmb{\epsilon} ] = \sigma ^ {\, 2} {\bf I}[/latex], and [latex]{\bf I}[/latex] is the [latex]n \times n[/latex] identity matrix.

• Theorem 1.1. The least squares estimators of [latex]\pmb{\beta}[/latex], denoted by [latex]\hat{\pmb{\beta}} = \big( \hat \beta_0 , \, \hat \beta_1 \big) ^ \prime[/latex], solve the normal equations

  $$\begin{array}{l}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}\hat{\mathit{\beta }}={\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y}.\end{array}$$

  The [latex]{\bf X}[/latex] matrix has rank 2 because there are at least two distinct X_i values. So [latex]{\bf X} ^ \prime {\bf X}[/latex] is invertible and the normal equations have the unique solution

  $$\begin{array}{l}\hat{\mathit{\beta }}=({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y},\end{array}$$

  by premultiplying both sides of the normal equations by [latex]\big( {\bf X} ^ \prime {\bf X} \big) ^ {-1}[/latex].

• Theorem 1.2. The least squares estimator of [latex]\pmb{\beta}[/latex] in a simple linear regression model is an unbiased estimator of [latex]\pmb{\beta}[/latex] because

  $$\begin{array}{l}E[\hat{\mathit{\beta }}]=\mathit{\beta }.\end{array}$$

• Theorem 1.3. The least squares estimators of [latex]\pmb{\beta}[/latex] in the simple linear regression model can be written as linear combinations of the dependent variables:

  $$\begin{array}{l}\hat{\mathit{\beta }}=({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y},\end{array}$$

  where the coefficients in the linear combinations are given by [latex]\big( {\bf X} ^ \prime {\bf X} \big) ^ {-1} {\bf X} ^ \prime[/latex].

• Theorem 1.4. The variance–covariance matrix of the least squares estimators of [latex]\pmb{\beta}[/latex] is

  $$\begin{array}{l}{\sigma }^2({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}.\end{array}$$

• Theorem 1.5 (Gauss–Markov theorem). The least squares estimators of [latex]\pmb{\beta}[/latex] in a simple linear regression model, [latex]\hat{\pmb{\beta}} = ( {\bf X} ^ \prime {\bf X} ) ^ {-1} {\bf X} ^ \prime {\bf Y}[/latex], have the smallest population variance amongst all linear unbiased estimators of [latex]\pmb{\beta}[/latex].
• Definition 1.2. The vector of fitted values in a simple linear regression model is the [latex]n \times 1[/latex] column vector

  $$\begin{array}{l}\widehat{\mathbf{\text{Y}}}=\mathbf{X}\hat{\mathit{\beta }}=\mathbf{X}({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y},\end{array}$$

  which is a linear combination of the dependent variables. The vector of residuals is the [latex]n \times 1[/latex] column vector

  $$\begin{array}{ll}\mathbf{e}& =\mathbf{Y}-\widehat{\mathbf{\text{Y}}}\\ & =\mathbf{Y}-\mathbf{X}\hat{\mathit{\beta }}\\ & =\mathbf{Y}-\mathbf{X}({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y}\\ & =(\mathbf{I}-\mathbf{X}({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}})\mathbf{Y},\end{array}$$

  which is also a linear combination of the dependent variables. The matrix [latex]{\bf I}[/latex] is the [latex]n \times n[/latex] identity matrix.

• Theorem 1.6. For the simple linear regression model with fitted values [latex]\bf\hat{Y}[/latex] and residuals [latex]{\bf e}[/latex],
  • [latex]{\bf e} ^ \prime {\bf 1} = 0[/latex],
  • [latex]{\bf Y} ^ \prime {\bf 1} = \bf\hat{Y}^ {\kern 0.19em \prime} {\bf 1}[/latex]
  • [latex]\bf\hat{Y}^ {\kern 0.19em \prime} {\bf e} = 0[/latex],

  where [latex]{\bf 1}[/latex] is an n-element column vector of ones.

• Theorem 1.7. An unbiased estimator of σ² in a simple linear regression model is

  $$\begin{array}{l}{\hat{\sigma }}^2=MSE=\frac{{\mathbf{e}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{e}}{n-2}.\end{array}$$

• Theorem 1.8. The sums of squares can be partitioned in a simple linear regression model as [latex]SST = SSR + SSE[/latex] or

  $$\begin{array}{l}(\mathbf{Y}-\overline{\mathbf{Y}}{)}^{\mathrm{\prime <!--\; \prime \; -->}}(\mathbf{Y}-\overline{\mathbf{Y}})=(\widehat{\mathbf{\text{Y}}}-\overline{\mathbf{Y}}{)}^{\mathrm{\prime <!--\; \prime \; -->}}(\widehat{\mathbf{\text{Y}}}-\overline{\mathbf{Y}})+(\mathbf{Y}-\widehat{\mathbf{\text{Y}}}{)}^{\mathrm{\prime <!--\; \prime \; -->}}(\mathbf{Y}-\widehat{\mathbf{\text{Y}}}),\end{array}$$

  where [latex]\bar {\bf Y}[/latex] is an n-element column vector with identical elements which are each the sample mean of the values of the dependent variable.

• Definition 1.3. The coefficient of determination in a simple linear regression model is

  $$\begin{array}{l}{R}^2=\frac{SSR}{SST}=\frac{(\widehat{\mathbf{\text{Y}}}-\overline{\mathbf{Y}}{)}^{\mathrm{\prime <!--\; \prime \; -->}}(\widehat{\mathbf{\text{Y}}}-\overline{\mathbf{Y}})}{(\mathbf{Y}-\overline{\mathbf{Y}}{)}^{\mathrm{\prime <!--\; \prime \; -->}}(\mathbf{Y}-\overline{\mathbf{Y}})},\end{array}$$

  when [latex]\big( {\bf Y} - \bar {\bf Y} \big) ^ \prime \big( {\bf Y} - \bar {\bf Y} \big) \ne 0[/latex]. The coefficient of correlation is

  $$\begin{array}{l}r=\pm \sqrt{R^2},\end{array}$$

  where the sign associated with r is positive when [latex]\hat \beta_1 \ge 0[/latex] and negative when [latex]\hat \beta_1 < 0[/latex].

• Definition 2.1. The simple linear regression model with normal error terms is

  $$\begin{array}{l}\mathbf{Y}=\mathbf{X}\mathit{\beta }+\mathit{ϵ},\end{array}$$

  where [latex]\pmb{\epsilon} \sim N\left( {\bf 0}, \, \sigma ^ {\, 2} {\bf I} \right)[/latex].

• Theorem 2.1. For the simple linear regression model with normal error terms, the maximum likelihood estimators of [latex]\pmb{\beta}[/latex] are

  $$\begin{array}{l}\hat{\mathit{\beta }}=({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y}\end{array}$$

  and the maximum likelihood estimator of σ² is

  $$\begin{array}{l}{\hat{\sigma }}^2=\frac{1}{n}(\mathbf{Y}-\mathbf{X}\hat{\mathit{\beta }}{)}^{\mathrm{\prime <!--\; \prime \; -->}}(\mathbf{Y}-\mathbf{X}\hat{\mathit{\beta }}).\end{array}$$

  Since the vector of error terms [latex]\pmb{\epsilon}[/latex] consists of independent and identically distributed normal random variables, [latex]{\bf Y} = {\bf X} \pmb{\beta} + \pmb{\epsilon}[/latex] is a vector of independent and identically distributed normal random variables, and the linear transformation [latex]\hat{\pmb{\beta}} = ( {\bf X} ^ \prime {\bf X} ) ^ {-1} {\bf X} ^ \prime {\bf Y}[/latex] has normally distributed elements.

• Theorem 2.2. For the simple linear regression model with normal error terms,

  $$\begin{array}{l}\frac{{\mathbf{e}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{e}}{{\sigma }^2}\sim {\chi }^2(n-2),\end{array}$$

  and is independent of [latex]\hat{\pmb{\beta}}[/latex].

• Theorem 2.3. For the simple linear regression model with normal error terms, an exact two-sided [latex]{100(1 - \alpha)}\%[/latex] confidence interval for σ² is

  $$\begin{array}{l}\frac{{\mathbf{e}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{e}}{{\chi }_{n-2,\alpha /2}^2}<{\sigma }^2<\frac{{\mathbf{e}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{e}}{{\chi }_{n-2,1-\alpha /2}^2}.\end{array}$$

• Theorems 2.4 and 2.7. For the simple linear regression model with normal error terms,

  $$\begin{array}{l}\hat{\mathit{\beta }}\sim N(\mathit{\beta },{\sigma }^2({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}).\end{array}$$

• Theorem 2.12. For the simple linear regression model with normal error terms, an exact two-sided [latex]{100(1 - \alpha)}\%[/latex] confidence interval for [latex]E[ Y_h ][/latex] for a given value of the independent variable X_h is

  $$\begin{array}{l}{\mathbf{X}}_h^{\mathrm{\prime <!--\; \prime \; -->}}\hat{\mathit{\beta }}-t_{n-2,\alpha /2}\sqrt{{\hat{\sigma }}^2{\mathbf{X}}_h^{\mathrm{\prime <!--\; \prime \; -->}}({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}_h}<E[Y_h]<{\mathbf{X}}_h^{\mathrm{\prime <!--\; \prime \; -->}}\hat{\mathit{\beta }}+t_{n-2,\alpha /2}\sqrt{{\hat{\sigma }}^2{\mathbf{X}}_h^{\mathrm{\prime <!--\; \prime \; -->}}({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}_h},\end{array}$$

  where [latex]{\bf X}_h = (1, \, X_h) ^ \prime[/latex] and [latex]\hat{\sigma} ^ 2 = MSE[/latex].

• Theorem 2.15. For the simple linear regression model with normal error terms, an exact two-sided [latex]{100(1 - \alpha)}\%[/latex] prediction interval for [latex]Y_h^\star[/latex] for a given value of the independent variable X_h is

  $$\begin{array}{l}{\mathbf{X}}_h^{\mathrm{\prime <!--\; \prime \; -->}}\hat{\mathit{\beta }}-t_{n-2,\alpha /2}\sqrt{{\hat{\sigma }}^2(1+{\mathbf{X}}_h^{\mathrm{\prime <!--\; \prime \; -->}}({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}_h)}<Y_h^{\star }<{\mathbf{X}}_h^{\mathrm{\prime <!--\; \prime \; -->}}\hat{\mathit{\beta }}+t_{n-2,\alpha /2}\sqrt{{\hat{\sigma }}^2(1+{\mathbf{X}}_h^{\mathrm{\prime <!--\; \prime \; -->}}({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}_h)},\end{array}$$

  where [latex]{\bf X}_h = (1, \, X_h) ^ \prime[/latex] and [latex]\hat{\sigma} ^ 2 = MSE[/latex].

• Theorem 2.16. Under the simple linear regression model with normal error terms and parameters estimated from the data pairs [latex]\left( X_1, \, Y_1 \right), \, \left( X_2, \, Y_2 \right), \, \ldots , \, \left( X_n, \, Y_n \right)[/latex],

  $$\begin{array}{l}\frac{(\hat{\mathit{\beta }}-\mathit{\beta }{)}^{\mathrm{\prime <!--\; \prime \; -->}}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}(\hat{\mathit{\beta }}-\mathit{\beta })}{2\cdot MSE}\sim F(2,n-2).\end{array}$$

• Theorem 2.17. Under the simple linear regression model with normal error terms and parameters estimated from the data pairs [latex]\left( X_1, \, Y_1 \right), \, \left( X_2, \, Y_2 \right), \, \ldots , \, \left( X_n, \, Y_n \right)[/latex], the values of β₀ and β₁ satisfying

  $$\begin{array}{l}\frac{(\hat{\mathit{\beta }}-\mathit{\beta }{)}^{\mathrm{\prime <!--\; \prime \; -->}}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}(\hat{\mathit{\beta }}-\mathit{\beta })}{2\cdot MSE}\le F_{2,n-2,\alpha }\end{array}$$

  form an exact joint [latex]{100(1 - \alpha)}\%[/latex] confidence region for β₀ and β₁.

• Definition 3.2. Under the simple linear regression model, the hat matrix is

  $$\begin{array}{l}\mathbf{H}=\mathbf{X}{\left({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}.\end{array}$$

  The diagonal elements of the hat matrix are the leverages. The matrix equation

  $$\begin{array}{l}\widehat{\mathbf{\text{Y}}}=\mathbf{H}\mathbf{Y}\end{array}$$

  indicates that [latex]{\bf H}[/latex] transforms [latex]{\bf Y}[/latex] to [latex]\bf\hat{Y}[/latex]. The hat matrix is symmetric (that is, [latex]{\bf H} = {\bf H} ^ \prime[/latex]) and idempotent (that is, [latex]{\bf H} {\bf H} = {\bf H}[/latex]).

The matrix approach applied to a simple linear regression model is illustrated for a small sample size next.

Theorem 2.2 stated that under the simple linear regression model with normal errors,

$$\begin{array}{l}\frac{SSE}{{\sigma }^2}\sim {\chi }^2(n-2).\end{array}$$

An outline of the proof of Theorem 2.2 was given in Chapter 2 in purely algebraic terms. An outline of the proof to the result using the matrix approach to simple linear regression is given here to contrast the difference between the two approaches.

The matrix approach gives an alternative way of computing measures of interest in a simple linear regression. Using matrices also allows the following two helpful extensions to simple linear regression.

• Removing the first column of the [latex]{\bf X}[/latex] matrix that consists entirely of ones corresponds to forcing a regression line through the origin.
• Adding additional columns to the [latex]{\bf X}[/latex] matrix corresponds to including additional independent variables to the regression model, which is known as multiple linear regression. This is the topic of the next section.

3.5 Multiple Linear Regression

Multiple linear regression can often be applied when there are several independent variables (or predictors) [latex]X_1, \, X_2, \, \ldots, \, X_p[/latex] which can be used to explain a continuous dependent (or response) variable Y. Three examples are listed below.

• The asking price of a home Y is a function of
  • the number of square feet in the home,
  • the number of bedrooms, and
  • acreage of the land associated with the home.
• The annual amount of money a person donates to charity Y is a function of
  • the nationality of the person,
  • the annual income of the person,
  • the net worth of the person,
  • the religious affiliation of the person,
  • the age of the person, and
  • the gender of the person.
• The stopping distance of a car Y is a function of
  • the speed of the car,
  • the weight of the car, and
  • the type of brakes installed on the car.

One way to formulate a multiple linear regression model is to treat the left-hand side of the model as an expected value:

$$\begin{array}{l}E[Y]={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_p{X}_p.\end{array}$$

Since [latex]E[Y][/latex] denotes a conditional expectation of Y given the values of the p independent variables [latex]X_1, \, X_2, \, \ldots, \, X_n[/latex], a more careful way to write this model is

$$\begin{array}{l}E[Y|X_1,X_2,\dots ,X_n]={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_p{X}_p.\end{array}$$

So far, there has been no consideration of the probability distribution of the error terms, and that is addressed in the formal definition of a multiple linear regression model given next.

To estimate the parameters in a multiple linear regression model, we collect n observations which each consist of the p independent variables and the associated dependent variable. In most applications, [latex]p > n[/latex]. Occasions arise (often in biostatistical applications) in which [latex]p > n[/latex]. The formulation of the simple linear regression model with notation included for the n observations is

$$\begin{array}{l}{Y}_i={\beta }_0+{\beta }_1{X}_{i1}+{\beta }_2{X}_{i2}+\cdots +{\beta }_p{X}_{ip}+{ϵ}_i\end{array}$$

for [latex]i = 1, \, 2, \, \ldots, \, n[/latex]. So [latex]X_{ij}[/latex] denotes the value of the jth independent variable collected on the ith observational unit. In the real estate example given at the beginning of this section, X₈₃ is the value of the third independent variable (acreage) collected on the 8th home collected by the analyst. The associated asking price of the 8th home is Y₈.

Figure 3.15 shows a portion of the population regression plane [latex]E[Y] = \beta_0 + \beta_1 X_1 + \beta_2 X_2[/latex] for a multiple linear regression model with [latex]p = 2[/latex] independent variables X₁ and X₂. The plane extends outward from the portion shown in Figure 3.15. The regression parameters β₀, β₁, and β₂ are fixed constants. The intercept β₀ is positive in Figure 3.15 because the plane strikes the Y-axis above the origin. Based on the inclination of the population regression plane relative to the X₁– and X₂-axes it is clear that [latex]\beta_1 < 0[/latex] and [latex]\beta_2 > 0[/latex]. To avoid clutter and highlight the geometry and notation, only the ith data triple [latex](X_{i1}, \, X_{i2}, \, Y_i)[/latex] and the associated error term ϵ_i are shown in the figure.

Figure 3.16 shows a portion of the estimated regression plane [latex]Y = \hat \beta_0 + \hat \beta_1 X_1 + \hat \beta_2 X_2[/latex] for a multiple linear regression model with [latex]p = 2[/latex] independent variables X₁ and X₂. The estimated regression parameters [latex]\hat \beta_0[/latex], [latex]\hat \beta_1[/latex], and [latex]\hat \beta_2[/latex] are random variables which are estimated from n data triples [latex]\left( X_{11}, \, X_{12} , \, Y_1 \right), \, \left( X_{21}, \, X_{22}, \, Y_2 \right), \, \ldots , \, \left( X_{n1}, \, X_{n2}, \, Y_n \right)[/latex]. The estimated regression parameters are random variables because the dependent variable values [latex]Y_1, \, Y_2, \, \ldots, \, Y_n[/latex] are random variables. The estimated intercept [latex]\hat \beta_0[/latex] is positive in Figure 3.16 because the plane strikes the Y-axis above the origin. Based on the inclination of the estimated regression plane relative to the X₁– and X₂-axes it is clear that [latex]\hat \beta_1 < 0[/latex] and [latex]\hat \beta_2 > 0[/latex].  To avoid clutter and highlight the geometry and notation, just the ith data triple [latex](X_{i1}, \, X_{i2}, \, Y_i)[/latex], the associated fitted value [latex](X_{i1}, \, X_{i2}, \, \hat{Y}_i)[/latex], and the associated residual e_i are shown in the figure.

When there are [latex]p > 2[/latex] independent variables, the estimated regression model is a hyperplane in [latex]{\cal R} ^ {p + 1}[/latex]. Residual i is the distance [latex]e_i = Y_i - \hat Y_i[/latex], for [latex]i = 1, \, 2, \, \ldots, \, n[/latex].

When the error terms are assumed to be normally distributed, this model is known as the multiple linear regression model with normal error terms. This additional assumption allows for statistical inference concerning parameters and predicted values in a similar manner to that described in Chapter 2.

The multiple linear regression model can also be expressed in terms of matrices. Relative to the simple linear regression model, additional columns are appended to the [latex]{\bf X}[/latex] matrix, and the [latex]\pmb{\beta}[/latex] vector is expanded to include the parameters associated with the additional parameters:

$$\begin{array}{l}\mathbf{X}=\left[\begin{array}{ccccc}1& X_{11}& X_{12}& \cdots & X_{1p}\\ 1& X_{21}& X_{22}& \cdots & X_{2p}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& X_{n1}& X_{n2}& \cdots & X_{np}\end{array}\right],\mathbf{Y}=\left[\begin{array}{c}{Y}_1\\ Y_2\\ ⋮\\ Y_n\end{array}\right],\mathit{\beta }=\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_p\end{array}\right],\text{and}\mathit{ϵ}=\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ ⋮\\ {ϵ}_n\end{array}\right].\end{array}$$

The vectors [latex]{\bf Y}[/latex] and [latex]\pmb{\epsilon}[/latex] remain unchanged from the simple linear regression formulation. The first row of [latex]{\bf X}[/latex] corresponds to the values of the independent variables collected on the first observational unit, the second row of [latex]{\bf X}[/latex] corresponds to the values of the independent variables collected on the second observational unit, etc. As was the case in simple linear regression, [latex]{\bf X}[/latex] is known as the design matrix.

The good news about the matrix approach to multiple linear regression is that the definitions and results from simple linear regression only require some minor tweaking in order to generalize to multiple regression. Several of these definitions and results are given below. In many cases, it is just a matter of replacing the word “simple” with the word “multiple” or updating the degrees of freedom to account for the p independent variables. It is assumed that the [latex]{\bf X}[/latex] matrix has rank [latex]p + 1[/latex] (that is, a full rank matrix), which means that the columns of [latex]{\bf X}[/latex] are linearly independent.

• The multiple linear regression model is

  $$\begin{array}{l}\mathbf{Y}=\mathbf{X}\mathit{\beta }+\mathit{ϵ},\end{array}$$

  where [latex]E[ \pmb{\epsilon} ] = {\bf 0}[/latex], [latex]V[ \pmb{\epsilon} ] = \sigma ^ {\, 2} {\bf I}[/latex], and [latex]{\bf I}[/latex] is the [latex]n \times n[/latex] identity matrix.

• The least squares estimators of [latex]\pmb{\beta}[/latex], denoted by [latex]\hat{\pmb{\beta}} = \big( \hat \beta_0 , \, \hat \beta_1 , \, \ldots , \, \hat \beta_p \big) ^ \prime[/latex], solve the normal equations

  $$\begin{array}{l}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}\hat{\mathit{\beta }}={\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y}.\end{array}$$

  Since [latex]{\bf X}[/latex] has full rank, [latex]{\bf X} ^ \prime {\bf X}[/latex] is invertible and the normal equations have the unique solution

  $$\begin{array}{l}\hat{\mathit{\beta }}=({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y},\end{array}$$

  by premultiplying both sides of the normal equations by [latex]\big( {\bf X} ^ \prime {\bf X} \big) ^ {-1}[/latex].

• The least squares estimator of [latex]\pmb{\beta}[/latex] in a multiple linear regression model is an unbiased estimator of [latex]\pmb{\beta}[/latex] because

  $$\begin{array}{l}E[\hat{\mathit{\beta }}]=\mathit{\beta }.\end{array}$$

• The least squares estimators of [latex]\pmb{\beta}[/latex] in the multiple linear regression model can be written as linear combinations of the dependent variables:

  $$\begin{array}{l}\hat{\mathit{\beta }}=({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y},\end{array}$$

  where the coefficients in the linear combinations are given by [latex]\big( {\bf X} ^ \prime {\bf X} \big) ^ {-1} {\bf X} ^ \prime[/latex].

• The variance–covariance matrix of the least squares estimators of [latex]\pmb{\beta}[/latex] is

  $$\begin{array}{l}{\sigma }^2({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}.\end{array}$$

• (Gauss–Markov theorem) The least squares estimators of [latex]\pmb{\beta}[/latex] in a multiple linear regression model, [latex]\hat{\pmb{\beta}} = ( {\bf X} ^ \prime {\bf X} ) ^ {-1}{\bf X} ^ \prime {\bf Y}[/latex], have the smallest population variance amongst all linear unbiased estimators of [latex]\pmb{\beta}[/latex].
• The vector of fitted values in a multiple linear regression model is the [latex]n \times 1[/latex] column vector

  $$\begin{array}{l}\widehat{\mathbf{\text{Y}}}=\mathbf{X}\hat{\mathit{\beta }}=\mathbf{X}({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y},\end{array}$$

  which is a linear combination of the dependent variables. The vector of residuals is the [latex]n \times 1[/latex] column vector

  $$\begin{array}{ll}\mathbf{e}& =\mathbf{Y}-\widehat{\mathbf{\text{Y}}}\\ & =\mathbf{Y}-\mathbf{X}\hat{\mathit{\beta }}\\ & =\mathbf{Y}-\mathbf{X}({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y}\\ & =(\mathbf{I}-\mathbf{X}({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}})\mathbf{Y},\end{array}$$

  which is also a linear combination of the dependent variables. The matrix [latex]{\bf I}[/latex] is the [latex]n \times n[/latex] identity matrix.

• The multiple linear regression model with normal error terms is

  $$\begin{array}{l}\mathbf{Y}=\mathbf{X}\mathit{\beta }+\mathit{ϵ},\end{array}$$

  where [latex]\pmb{\epsilon} \sim N\left( {\bf 0}, \, \sigma ^ {\, 2} {\bf I} \right)[/latex].

• For the multiple linear regression model with normal error terms, the maximum likelihood estimators of [latex]\pmb{\beta}[/latex] are

  $$\begin{array}{l}\hat{\mathit{\beta }}=({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y}\end{array}$$

  and the maximum likelihood estimator of σ² is

  $$\begin{array}{l}{\hat{\sigma }}^2=\frac{1}{n}(\mathbf{Y}-\mathbf{X}\hat{\mathit{\beta }}{)}^{\mathrm{\prime <!--\; \prime \; -->}}(\mathbf{Y}-\mathbf{X}\hat{\mathit{\beta }}).\end{array}$$

  Since the vector of error terms [latex]\pmb{\epsilon}[/latex] consists of independent and identically distributed normal random variables, [latex]{\bf Y} = {\bf X} \pmb{\beta} + \pmb{\epsilon}[/latex] is a vector of independent and identically distributed normal random variables. Since [latex]\hat{\pmb{\beta}}[/latex] is a linear transformation of Y, [latex]\hat{\pmb{\beta}} \sim N \left( {\pmb{\beta}}, \, \sigma ^ {\, 2} \big( {\bf X} ^ \prime {\bf X} \big) ^ {-1} \right)[/latex].

• Under the multiple linear regression model, the [latex]n \times n[/latex] hat matrix is

  $$\begin{array}{l}\mathbf{H}=\mathbf{X}{\left({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}.\end{array}$$

  The diagonal elements of the hat matrix are the leverages. The matrix equation

  $$\begin{array}{l}\widehat{\mathbf{\text{Y}}}=\mathbf{H}\mathbf{Y}\end{array}$$

  indicates that [latex]{\bf H}[/latex] transforms [latex]{\bf Y}[/latex] to [latex]\bf\hat{Y}[/latex]. The hat matrix is symmetric (that is, [latex]{\bf H} = {\bf H} ^ \prime[/latex]) and idempotent (that is, [latex]{\bf H} {\bf H} = {\bf H}[/latex]). The trace of the hat matrix is [latex]\sum_{i\,=\,1}^n h_{ii} = p + 1[/latex].

The example of multiple linear regression that follows considers [latex]p = 2[/latex] predictors of the sales price of a home.

A multiple linear regression model can easily be adapted to include nonlinear terms. A multiple regression model with two independent variables X₁ and X₂, for example, with a linear relationship between X₁ and Y and a quadratic relationship between X₂ and Y which includes an intercept term is

$$\begin{array}{l}Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_2^2+ϵ.\end{array}$$

Using the R lm function to estimate the coefficients will be illustrated in Section 3.7.

Multiple linear regression has many more modeling issues that arise than simple linear regression. The subsections that follow consider the following topics within multiple regression: (a) handling categorical independent variables which fall in categories rather than quantitative values, (b) handling the case in which independent variables have interactive effects, (c) extending the ANOVA table to multiple independent variables, (d) calculation of the coefficient of determination for multiple linear regression, and an adjustment that can be made to reduce its bias, (e) the effect of multicollinearity among the independent variables, and (f) algorithms for model selection.

3.5.2 Interaction Terms

The multiple linear regression model

$$\begin{array}{l}Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_p{X}_p+ϵ\end{array}$$

assumes a linear relationship between each independent variable and Y and the slope associated with an independent variable is identical at all values of the other independent variables within the scope of the multiple linear regression model. This relationship is illustrated for some selected data points of smaller homes from the Ames, Iowa housing data set from Examples 2.9 and 3.10. In this case, X₁ is the interior square footage, X₂ is an indicator variable reflecting the lot size,

$$\begin{array}{l}{X}_2=\{\begin{array}{lll}0& & \text{lot size is less than or equal to 10,000 square feet}\\ 1& & \text{lot size is greater than 10,000 square feet,}\end{array}\end{array}$$

and Y is the sales price. The multiple linear regression model with the [latex]p = 2[/latex] independent variables is

$$\begin{array}{l}Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+ϵ.\end{array}$$

Figure 3.17 shows a scatterplot of the interior square footage and sales price of homes on smaller lots ([latex]X_2 = 0[/latex] as open points) and larger lots ([latex]X_2 = 1[/latex] as solid points). The values of [latex]\hat \beta_0[/latex], [latex]\hat \beta_1[/latex], and [latex]\hat \beta_2[/latex] are indicated on the graph. The estimated intercept [latex]\hat \beta_0 = \text{21,473}[/latex], although slightly outside of the scope of the model, gives the estimated sales price of a small lot containing no dwelling as $21,473. The estimated regression coefficient [latex]\hat \beta_1 = 31.33[/latex] indicates that the sales price of a home increases by an estimated $31.33 for each additional interior square foot, adjusted for lot size. The estimated regression coefficient [latex]\hat \beta_2 = \text{35,693}[/latex] indicates that homes on larger lots cost $35,693 more, on average, than homes on smaller lots, adjusted for interior square feet. Notice that this formulation of the multiple linear regression model forces the slopes of the two lines in Figure 3.17 to be identical, regardless of the value of X₂.

Long Description for Figure 3.17

The horizontal axis labeled X 1 and ranges between 0 to 2500 in increments of 500. The vertical axis labeled Y ranges from 0 to 120,000 in increments of 40,000. The data pairs representing sales price of smaller lots, X equals 0, are indicated in open circle sand those for larger lots, X 2 equals 1 are represented as solid points. Two parallel regression lines, with positive slopes originate from points (0, 20,000) and (0, 50,000). The data pair for smaller lots are clustered around the regression line originating at (0, 20,000) and those for larger lots are clustered around the regression line originating at (0, 50,000). The slope of the regression lines are indicated as beta cap 1. The distance between horizontal axis and regression line originating at (0, 20,000) is indicated as beta cap 0. The distance between the two regression lines is indicated as beta cap 2.

But is the assumption of equal slopes of the two lines in Figure 3.17 justified? Separate simple linear regression models are fitted to the homes built on smaller and larger lots, and the results are plotted in Figure 3.18. The lines do not appear to be parallel in this case, indicating that a more complex regression model is warranted. There appears, in this case, to be an interaction effect between X₁ and X₂. This means that the effect of one independent variable (X₁, for example, the interior size) on Y is altered based on the value of another independent variable (X₂, the lot size indicator).

Long Description for Figure 3.18

The horizontal axis labeled X 1 ranges between 0 to 2500 in increments of 500. The vertical axis labeled Y ranges from 0 to 120,000 in increments of 40,000. The data pairs representing the sales price of smaller lots, X equals 0, are indicated in open circle sand those for larger lots, X 2 equals 1 are represented as solid points. Two regression lines, with positive slopes originate from points (0, 30,000) and (0, 40,000). Most of the data pair fall between 500 and 1500 of horizontal axis and 30,000 and 60,000 of vatical axis values. The data pair for smaller lots are clustered around the less steeper regression line originating at (0, 30,000). The data pairs for larger lots are clustered around the regression line, which is a steeper line, originating at (0, 40,000). Most of the data pairs lie between 1000 and 2000 on the horizontal axis and between 70,000 and 120,000 on the vertical axis.

Regression analysts account for this interaction by including cross-product terms in the regression model. In this Ames housing data set example, the regression model with an interaction term is

$$\begin{array}{l}Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_1{X}_2+ϵ.\end{array}$$

If the regression parameter [latex]\hat \beta_3[/latex] differs statistically from 0, then the inclusion of the interaction term is warranted. Notice that when [latex]X_2 = 0[/latex] (smaller lots), the model reduces to

$$\begin{array}{l}Y={\beta }_0+{\beta }_1{X}_1+ϵ,\end{array}$$

which is a simple linear regression model with intercept parameter β₀ and slope parameter β₁. On the other hand, when [latex]X_2 = 1[/latex] (larger lots), the model reduces to

$$\begin{array}{l}Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2+{\beta }_3{X}_1+ϵ\end{array}$$

or

$$\begin{array}{l}Y={\beta }_0+{\beta }_2+({\beta }_1+{\beta }_3)X_1+ϵ\end{array}$$

which is a simple linear regression model with intercept parameter [latex]\beta_0 + \beta_2[/latex] and slope parameter [latex]\beta_1 + \beta_3[/latex]. It is in this fashion that the two non-parallel lines depicted in Figure 3.18 can be estimated in a single regression model. Not surprisingly, it requires four parameters, β₀, β₁, β₂, and β₃, to do so. The multiple linear regression model with an interaction term can be fitted using the lm function in R by simply replacing the usual + in the formula with *. All four parameters are statistically significant at the 0.05 level in this case, so the inclusion of an interaction term is warranted.

3.5.5 Multicollinearity

In many settings, the values of the independent variables are correlated. In the housing data set from Example 3.10, for example, the independent variables X₁ (interior square footage) and X₂ (lot size) are probably positively correlated. Intuition suggests that larger homes are built on larger lots, on average. In the extreme case, what if homes in Ames were required by some bizarre municipal code to all be single story homes with the square footage of the lot always exactly four times the square footage of the interior of the home? In this case, [latex]X_2 = 4 X_1[/latex], so knowing the value of either X₁ or X₂ allows you to know the value of the other. Intuitively, one of the two independent variables is superfluous. When this is the case, the design matrix [latex]{\bf X}[/latex] has two columns which are multiples of one another, so these columns are linearly dependent and the matrix does not have full rank. This implies that the matrix [latex]{\bf X} ^ \prime {\bf X}[/latex] (which is used in computing the estimates of the regression coefficients) is singular, so it does not have an inverse. In this case, the usual formula for the regression coefficients,

$$\begin{array}{l}\hat{\mathit{\beta }}={\left({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{Y},\end{array}$$

is undefined because the matrix [latex]{\bf X} ^ \prime {\bf X}[/latex] does not have an inverse. In the case in which [latex]X_2 = 4 X_1[/latex], all pairs of the independent variables fall on a line, so it is impossible to know the proper tilt of the fitted regression plane in [latex]{\cal R} ^ 3[/latex]. There are many planes that minimize the sum of squared errors.

Multicollinearity is the condition associated with independent variables that are highly correlated among themselves in a multiple regression model. More specifically, multicollinearity occurs when two or more of the independent variables have a high correlation. This can appear as an approximately linear relationship between two of the independent variables. Multicollinearity is a condition associated with the design matrix [latex]{\bf X}[/latex] rather than the values of the dependent variable [latex]{\bf Y}[/latex] or the model [latex]{\bf Y} = {\bf X} {\pmb{\beta}} + {\pmb{\epsilon}}[/latex]. In cases in which multicollinearity exists, the matrix [latex]{\bf X} ^ \prime {\bf X}[/latex] has an inverse, but it is ill-conditioned and subject to slight variations in the data or is unstable because of large differences in the magnitudes of the various values of the independent variables. One of the key practical issues when multicollinearity is present is that an estimated regression coefficient for a particular independent variable depends on whether the other independent variables are included or left out of the model.

So multicollinearity has been loosely defined as high correlation among the independent variables. There is redundancy to the information contained in the independent variables. The next paragraphs describe how to detect multicollinearity, its consequences, and some remedies.

Although the hypothetical perfect correlation between the interior space and the lot size of a home from Ames, Iowa described previously occurs seldom in practice, highly correlated independent variables can result in some unusual behavior of regression coefficients as a regression model is constructed. Some signs that multicollinearity might be present in a multiple linear regression model include the following.

• Large values of the estimated standard deviations of the regression coefficients.
• Including or not including an independent variable in the model results in large changes to the estimated regression coefficients.
• An estimated regression coefficient that is statistically significant when the associated independent variable is considered alone, but becomes insignificant when one or more other independent variables are added to the model.
• An estimated regression coefficient with a sign that is inconsistent with expected sign or inconsistent with previous similar data sets.
• The pairwise sample correlation among the independent variables is high. The cor function in R can be used to assess the correlation among independent variables. The R statement
  

  

  for example, calculates the correlation matrix for the columns of the built-in data frame named swiss. The off-diagonal elements of this matrix range from [latex]-0.69[/latex] to 0.70, indicating that multicollinearity is present.

All of the criteria listed above are informal. A more formal way to determine whether multicollinearity is present is to introduce a statistic which reflects multicollinearity. The estimate of the variance of [latex]\hat \beta _ j[/latex] can be written as

$$\begin{array}{l}\hat{V}[{\hat{\beta }}_j]=\frac{1}{1-R_j^2}\left[\frac{MSE}{\sum _{i=1}^n(X_{ij}-{\overline{X}}_j{)}^2}\right],\end{array}$$

where [latex]\bar X_j = \sum_{i\,=\,1}^n X_{ij}[/latex], [latex]MSE = SSE / (n - p - 1)[/latex] for the full multiple regression model, and [latex]R_j^2[/latex] is the coefficient of determination obtained by conducting a multiple linear regression with X_j as the dependent variable and the other [latex]p - 1[/latex] X-values as the independent variables, for [latex]j = 1, \, 2, \, \ldots, \, p[/latex]. The coefficient on the right-hand side of this equation,

$$\begin{array}{l}VI{F}_j=\frac{1}{1-R_j^2},\end{array}$$

is known as a variance inflation factor for independent variable j, for [latex]j = 1, \, 2, \, \ldots, \, p[/latex]. In the extreme case when [latex]R_j^2 = 0[/latex], the associated variance inflation factor is [latex]VIF_j = 1[/latex]. This corresponds to the case in which X_j is not linearly related to the other independent variables. As [latex]R_j^2[/latex] increases, [latex]VIF_j[/latex] also increases, corresponding to increased correlation between the independent variables. When the largest of the [latex]VIF_j[/latex] values exceeds the threshold value of 10, one can conclude that the multicollinearity is present among the independent variables.

The R code below calculates the variance inflation factors for the data values in the swiss data frame, where the independent variables

• X₁, the percentage of males involved in agriculture as an occupation,
• X₂, the percentage of draftees receiving the highest make on an army examination,
• X₃, the percentage of draftees with education beyond the primary school,
• X₄, the percentage of Catholics, and
• X₅, the percentage of live births who live less than one year,

are used to predict Y, a common standardized fertility measure, from the [latex]n = 47[/latex] French-speaking provinces of Switzerland in about the year 1888. The R code below computes the variance inflation factors for the [latex]p = 5[/latex] independent variables.

The variance inflation factors for the [latex]p = 5[/latex] independent variables are

$$\begin{array}{l}VI{F}_1=2.28,VI{F}_2=3.68,VI{F}_3=2.77,VI{F}_4=1.94,VI{F}_5=1.11.\end{array}$$

Since none of these five values exceeds 10, we can conclude that the multicollinearity that exists in the independent variables is not strong enough to cause concern. (Some regression analysts use 5 as a threshold rather than 10.) Some keystrokes can be saved by using the vif function from the car package on a multiple linear regression model fitted by the lm function.

One popular remedy for multicollinearity is known as ridge regression, which is a parameter estimation technique that abandons the requirement of unbiased parameter estimates. The approach taken with ridge regression is to choose estimates for the regression parameters that are biased, but have a smaller variance than the ordinary least squares estimates. The goal is to generate parameter estimates with tolerable bias but smaller variance. The typical approach used in statistics to overcome this bias/variability trade-off is to use the estimates that minimize the mean square errors. Assuming that the X and Y values have been centered, we can dispense with the need for an intercept term in the multiple regression model. Rather than minimizing the usual sum of squared errors

$$\begin{array}{l}S=\sum _{i=1}^n{(Y_i-{\beta }_1{X}_{i1}-{\beta }_2{X}_{i2}-\cdots -{\beta }_p{X}_{ip})}^2,\end{array}$$

ridge regression minimizes

$$\begin{array}{l}{S}_R=\sum _{i=1}^n{(Y_i-{\beta }_1{X}_{i1}-{\beta }_2{X}_{i2}-\cdots -{\beta }_p{X}_{ip})}^2+\lambda \sum _{j=1}^p{\beta }_j^2.\end{array}$$

There are now two terms in the modified sum of squares. The second term in S_R is known as the penalty term. The new parameter λ is known as the penalty parameter. When [latex]\lambda = 0[/latex], [latex]S_R[/latex], reduces to the ordinary least squares case and achieves a value SSE at the ordinary least squares estimators. As λ increases, the estimators converge to [latex]\hat \beta_1 = \hat \beta_2 = \cdots = \beta_p = 0[/latex]. We desire a λ value that introduces some bias into the parameter estimates, but also have a reduced variance.

The geometry associated with ridge regression for [latex]p = 2[/latex] independent variables X₁ and X₂ in a multiple linear regression model is illustrated in Figure 3.19. The ellipses are level surfaces of the first term in [latex]S_R[/latex]. The center of the ellipses is the ordinary least squares estimators of [latex](\beta_1, \, \beta_2) = \big( \hat \beta_1, \, \hat \beta_2 \big)[/latex], which are the values that minimize the first term of [latex]S_R[/latex]. The circles centered at the origin are level surfaces of the second term in [latex]S_R[/latex]. The ridge regression estimators for β₁ and β₂ will occur at the intersection of one of elliptical and circular contours. In Figure 3.19 the two outermost level surfaces intersect at a point, which is a value of the ridge regression estimates of β₁ and β₂ which correspond to one particular value of the penalty parameter λ. The point at which this intersection occurs is a function of the penalty parameter λ. In higher dimensions, the circles become spheres and the ellipses become ellipsoids.

Long Description for Figure 3.19

The horizontal axis labeled beta 0 contains the points 0 and beta cap 0. The vertical axis labeled beta 1 contains the points 0 and beta cap 1. A point plotted at beta cap 0, beta cap 1 is surrounded by two concentric elliptical contours. Another point plotted at 0, 0 is surrounded by two concentric circles. The outermost levels of the elliptical and circle contours intersect at a point between 0 and beta cap 0 of the horizontal axis and 0 and beta cap 1 of the vertical axis.

Determining the value of the penalty parameter is critical in ridge regression, but its choice depends on the regression model and associated data set. A common technique for determining an optimal value for λ is known as k-fold cross-validation. There are several functions in R which can perform ridge regression: the lm.ridge function from the MASS package, the linearRidge function from the ridge package, and the glmnet function from the glmnet package. Ridge regression is related to the lasso (least absolute shrinkage and selection operator) estimator and elastic net regularization, two other popular parameter estimation techniques that are often applied for large values of p.

Is there a way to completely avoid multicollinearity? In some settings, the answer is yes. When the values of the independent variables are chosen so that they are uncorrelated, the regression coefficients associated with a simple linear regression model of each independent variable separately match the regression coefficients of any model involving more independent variables. This fact provides a strong argument for a designed experiment which can result in uncorrelated independent variables whenever the setting of the regression problem make this possible.

3.6 Weighted Least Squares

The three approaches to estimating the parameters in a simple linear regression model that we have encountered thus far,

• the algebraic approach,
• the matrix approach,
• using the R lm (linear model) function,

all have the same assumptions regarding the independent variable, the dependent variable, and the model [latex]Y = \beta_0 + \beta_1 X + \epsilon[/latex]. In all three approaches, the error terms are assumed to be mutually independent random variables, each with population mean 0 and population variance–covariance matrix [latex]V[\epsilon] = \sigma _ Z ^ {\, 2} I[/latex], where I is the [latex]n \times n[/latex] identity matrix. This means that [latex]V[ \epsilon_i ] = \sigma _ Z ^ {\, 2}[/latex], for [latex]i = 1, \, 2, \, \ldots, \, n[/latex]. There is also an implicit assumption that each of the data pairs [latex](X_i, \, Y_i)[/latex] are each given equal weight in the regression.

Settings occasionally arise in which some data values should be given different weights. There might be evidence that some of the Y_i values have more precision than others. Weights can be placed on each of the data pairs to account for this difference in precision. This leads to a weighted least squares approach to estimating the coefficients in a regression model.

In the standard simple linear regression model, the assumption

$$\begin{array}{l}V[{ϵ}_i]={\sigma }_Z^2,\end{array}$$

for [latex]i = 1, \, 2, \, \ldots, \, n[/latex], means that the variance of the dependent variable from the regression line is equal for all of the n data pairs, regardless of the value of the independent variable. In weighted least squares modeling, the positive weights [latex]w_1, \, w_2, \, \ldots , \, w_n[/latex] are determined so that

$$\begin{array}{l}V[{ϵ}_i]={\sigma }_Z^2/w_i\end{array}$$

for [latex]i = 1, \, 2, \, \ldots, \, n[/latex], which means that certain data pairs have more precision than other data pairs. The weights are fixed constants. There is no requirement that the weights sum to one. Data pairs with larger weights are assumed to have a lower variability to their error terms. This allows for a population variance that changes from one data pair to another.

As an illustration, the values of the dependent variable Y might be sample means at the various values of the independent variable X. Furthermore, if the sample sizes associated with the sample means are known and unequal, then we would like to assign higher weights to the data pairs associated with larger sample sizes. If n_i is the sample size for data pair i, for [latex]i = 1, \, 2, \, \ldots, \, n[/latex], then the appropriate weight for data pair i is [latex]w_i = n_i[/latex] so that

$$\begin{array}{l}V[{ϵ}_i]={\sigma }_Z^2/n_i\end{array}$$

for [latex]i = 1, \, 2, \, \ldots, \, n[/latex].

So rather than minimizing the sum of squares

$$\begin{array}{l}S=\sum _{i=1}^n(Y_i-{\beta }_0-{\beta }_1{X}_i{)}^2\end{array}$$

as was the case in the standard simple linear regression model, weighted least squares minimizes the weighted sum of squares

$$\begin{array}{l}S=\sum _{i=1}^n{w}_i(Y_i-{\beta }_0-{\beta }_1{X}_i{)}^2.\end{array}$$

Notice that this reduces to the ordinary sum of squares when [latex]w_1 = w_2 = \cdots = w_n = 1[/latex]. As before, calculus can be used to minimize S with respect to β₀ and β₁ to arrive at the least squares estimators [latex]\beta_0[/latex] and [latex]\beta_1[/latex]. The partial derivatives of S with respect to β₀ and β₁ are

$$\begin{array}{l}\frac{\partial S}{\partial {\beta }_0}=-2\sum _{i=1}^n{w}_i(Y_i-{\beta }_0-{\beta }_1{X}_i)=0\end{array}$$

and

$$\begin{array}{l}\frac{\partial S}{\partial {\beta }_1}=-2\sum _{i=1}^n{w}_i{X}_i(Y_i-{\beta }_0-{\beta }_1{X}_i)=0.\end{array}$$

These can be simplified to give the normal equations

$$\begin{array}{l}{\beta }_0\sum _{i=1}^n{w}_i+{\beta }_1\sum _{i=1}^n{w}_i{X}_i=\sum _{i=1}^n{w}_i{Y}_i\end{array}$$

and

$$\begin{array}{l}{\beta }_0\sum _{i=1}^n{w}_i{X}_i+{\beta }_1{w}_i{X}_i^2=\sum _{i=1}^n{w}_i{X}_i{Y}_i.\end{array}$$

The normal equations are a system of two linear equations in the two unknowns β₀ and β₁, given the data pairs [latex](X_1, \, Y_1), \, (X_2, \, Y_2), \, \ldots , \, (X_n, \, Y_n)[/latex] and the weights [latex]w_1, \, w_2, \, \ldots, \, w_n[/latex]. The normal equations can be solved to yield the weighted least squares estimators. This derivation constitutes a proof of the following theorem.

The matrix approach can also be applied to weighted least squares. Define the [latex]{\bf X}[/latex], [latex]{\bf Y}[/latex], [latex]\pmb{\beta}[/latex] and [latex]\pmb{\epsilon}[/latex] matrices as in Section 3.4:

$$\begin{array}{l}\mathbf{X}=\left[\begin{array}{cc}1& X_1\\ 1& X_2\\ ⋮& ⋮\\ 1& X_n\end{array}\right],\mathbf{Y}=\left[\begin{array}{c}{Y}_1\\ Y_2\\ ⋮\\ Y_n\end{array}\right],\mathit{\beta }=\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\end{array}\right],\text{and}\mathit{ϵ}=\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ ⋮\\ {ϵ}_n\end{array}\right].\end{array}$$

In addition, assume that the matrix [latex]{\bf W}[/latex] is a diagonal matrix with the weights [latex]w_1, \, w_2, \, \ldots, \, w_n[/latex] on the diagonal:

$$\begin{array}{l}\mathbf{W}=\left[\begin{array}{cccc}{w}_1& 0& \cdots & 0\\ 0& w_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & w_n\end{array}\right].\end{array}$$

In this case, the normal equations can be written in matrix form as

$$\begin{array}{l}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{W}\mathbf{X}\mathit{\beta }={\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{W}\mathbf{Y}.\end{array}$$

Pre-multiplying both sides of this equation by [latex]\left( {\bf X} ^ \prime {\bf W} {\bf X} \right) ^ {-1}[/latex] gives the least squares estimators for the regression parameters in matrix form as

$$\begin{array}{l}\hat{\mathit{\beta }}={\left({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{W}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{W}\mathbf{Y}.\end{array}$$

As before, the fitted values can also be written in matrix form as

$$\begin{array}{l}\widehat{\mathbf{\text{Y}}}=\mathbf{X}\hat{\mathit{\beta }}\end{array}$$

or

$$\begin{array}{l}\widehat{\mathbf{\text{Y}}}=\mathbf{X}{\left({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{W}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{W}\mathbf{Y}.\end{array}$$

The residuals [latex]e_i = Y_i - \hat{Y}_i[/latex] for [latex]i = 1, \, 2, \, \ldots, \, n[/latex], can also be written in matrix form as

$$\begin{array}{ll}\mathbf{e}& =\mathbf{Y}-\widehat{\mathbf{\text{Y}}}\\ & =\mathbf{Y}-\mathbf{X}\hat{\mathit{\beta }}\\ & =\mathbf{Y}-\mathbf{X}{\left({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{W}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{W}\mathbf{Y}\\ & =(\mathbf{I}-\mathbf{X}{\left({\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{W}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{\prime <!--\; \prime \; -->}}\mathbf{W})\mathbf{Y},\end{array}$$

where [latex]{\bf e}[/latex] is the column vector of residuals [latex]{\bf e} = (e_1, \, e_2, \, \ldots , \, e_n) ^ \prime[/latex]. These matrix results are summarized in the following theorem.

The algebraic approach, matrix approach, and R approach to weighted least squares problem will be illustrated in the next example. Establishing the weights [latex]w_1, \, w_2, \, \ldots, \, w_n[/latex] can be a nontrivial problem, and differs depending on the setting in which the weighted regression model is employed.

Using simple linear regression in the previous example, either weighted or unweighted, might not be the best approach. The dependent variable Y is the probability that an item of age X is functioning. This dependent variable must lie between 0 and 1, but the regression line could potentially fall outside of that range within the scope of the model. Two potential remedies are given in the next two sections: using a regression model with nonlinear terms such as X² or X³, or a survivor function of a lifetime model rather than a line, or a nonlinear model known as a logistic regression model, whose dependent variable necessarily lies between 0 and 1.

3.7 Regression Models with Nonlinear Terms

Regression models with nonlinear terms arise frequently in regression modeling. One simple example is polynomial regression. A quadratic regression model, for example, is

$$\begin{array}{l}Y={\beta }_0+{\beta }_{1}X+{\beta }_2{X}^2+ϵ,\end{array}$$

where β₀, β₁, and β₂ are the regression coefficients, and ϵ is a white noise term. This model is still linear in β₀, β₁, and β₂. One way to think about this model is to consider X and X² to be the [latex]p = 2[/latex] independent variables in a multiple regression model. The next example fits a quadratic model to the data pairs in which the independent variable X is the speed of an automobile and the dependent variable Y is its stopping distance.

Nonlinear regression modeling is not limited to just polynomial regression models. The next two examples fit the same data set concerning the national debt in the United States between 1970 and 2020 to a nonlinear regression model using two fundamentally different approaches. The first approach is to transform the nonlinear regression model to a linear regression model and then apply the standard techniques for parameter estimation to the transformed model. The second approach is to use numerical methods to minimize the sum of squares in the usual least squares fashion described previously.

There is a second approach to fitting an exponential regression model to the national debt data pairs that follows the standard approach to least squares estimation, which is given next.

One drawback that emerged from the survival function estimation example from the previous section (involving current status data) is that fitting a regression line results in a survival probability that can be negative or greater than one when extrapolated outside of the range of the independent variable in the data pairs. In addition, the estimated probability of survival at time zero for both the ordinary simple linear regression model and the weighted simple linear regression model seemed low. Typically, a brand-new item is not defective. A nonlinear regression function is an attractive alternative model in this particular setting. The next example combines a nonlinear regression model and weighted least squares estimators to provide an improved regression model.

3.8 Logistic Regression

Logistic regression is appropriate when the dependent variable Y can assume one of two values: zero and one. This is sometimes known as a binary or dichotomous response variable. For now, to keep the mathematics and interpretations simple, assume that there is a single predictor X. This is known as a simple logistic regression model, and is a special type of nonlinear regression model. Including multiple independent variables in a logistic regression model is a straightforward extension. For dichotomous data, instead of predicting 0 or 1, we predict the probability of getting a 1 [that is, [latex]P(Y = 1)[/latex]]. So we need a regression model that predicts values of the interval [latex][0, \, 1][/latex].

The following example will be used throughout this section to motivate the need for a special model to accommodate a binary dependent variable, and to illustrate the techniques for the estimation of the model parameters.

One of the initial considerations in developing a statistical model for the outcome of a field goal as a function of the length of the field goal attempt is to find a function that will only assume values between 0 and 1. A diagram that gives some guidance with regard to this function is to batch the data into 5-year increments. So the bins are all field goals that fall in the ranges [latex]20 \pm 2, \, 25 \pm 2, \, \ldots , \, 60 \pm 2[/latex]. This window is long enough so that the random sampling variability associated with nearby attempts is damped considerably, and yet short enough so that outcome patterns as a function of yardage are still apparent. The R code below batches the independent variable into the 5-yard increments and plots the estimated probability of success for attempts in each batch at its midpoint. This estimated probability is just the fraction of successful field goals within a particular range. Furthermore, the area of each point plotted is proportional to the number of attempts in that particular bin. For example, there were 79 attempts in the first bin (18–22 yards) and only 4 attempts in the last bin (58–62 yards). The R code below reads a data set off of the web that contains the results of [latex]n = 948[/latex] NFL field goal attempts during 2003. The data consists of columns that give the length of the field goal attempt and the outcome, failure ([latex]Y = 0[/latex]) or success ([latex]Y = 1[/latex]). The R code rounds each length to the nearest 5 yards, and plots the midpoint of the rounded field goal lengths versus the estimated probability of success.

While the performance of NFL field goal kickers varies from one kicker to the next, these points give us an idea of what we would like for a smooth regression function in this setting.

The results are shown in Figure 3.28. It is clear that the estimated probability of making a field goal decreases as the length of the field goal attempt increases, as one would expect. There is a strong relationship between the length of the field goal attempt and the probability of success. Our goal is to fit a nonlinear regression function to the raw data values that smooths the random sampling variability and can be used for the purpose of prediction.

When the dependent variable only takes on the values zero and one, the usual mean response function for the simple linear regression model

$$\begin{array}{l}Y={\beta }_0+{\beta }_{1}X+ϵ\end{array}$$

is

$$\begin{array}{l}E[Y]={\beta }_0+{\beta }_{1}X,\end{array}$$

where [latex]E[Y][/latex] denotes the conditional expected value of Y given a particular setting of the independent variable X. This mean response function does not limit the values of Y to just zero and one. With normally distributed error terms, this model would allow for Y values which could be less than 0 or greater than 1.

In logistic regression, this type of curve, regardless of whether it begins near one and ends near zero or it begins near zero and ends near one, is known as a sigmoidal response function. A natural choice for the sigmoidal response function is a cumulative distribution function associated with a random variable, or its complement (the survivor function). Three popular probability distributions whose cumulative distribution functions are used in logistic regression are the standard logistic distribution (also commonly called the logit model), the standard normal distribution (also commonly called the probit model), and the standard extreme value distribution (also commonly called the complimentary log-log model). These are described in the next paragraph.

The standard logistic distribution has probability density function

$$\begin{array}{l}f(x)=\frac{e^x}{{(1+e^x)}^2}-\mathrm{\infty }<x<\mathrm{\infty }\end{array}$$

and cumulative distribution function

$$\begin{array}{l}F(x)=\frac{e^x}{1+e^x}-\mathrm{\infty }<x<\mathrm{\infty }.\end{array}$$

The probability density function is symmetric about the population mean [latex]E[X] = 0[/latex] and has population variance [latex]V[X] = \pi ^ 2 / 3[/latex]. The standard normal distribution has probability density function

$$\begin{array}{l}f(x)=\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}-\mathrm{\infty }<x<\mathrm{\infty }\end{array}$$

and cumulative distribution function

$$\begin{array}{l}F(x)={\int }_{-\mathrm{\infty }}^{x}f(w)dw-\mathrm{\infty }<x<\mathrm{\infty }.\end{array}$$

The probability density function is also symmetric about the population mean [latex]E[X] = 0[/latex] and has population variance [latex]V[X] = 1[/latex]. The probability density function for the standard logistic distribution is similar in shape (that is, bell-shape) to that for the standard normal distribution, but has heavier tails. The symmetry of the probability density functions for the standard logistic distribution and the standard normal distribution limits the shape of the associated cumulative distribution function. A nonsymmetric distribution often provides a better fit. This leads to a search for a probability distribution with a nonsymmetric probability density function. One such probability distribution is the extreme value distribution. The standard extreme value distribution has probability density function

$$\begin{array}{l}f(x)=e^{x-e^x}-\mathrm{\infty }<x<\mathrm{\infty }\end{array}$$

and cumulative distribution function

$$\begin{array}{l}F(x)=1-e^{-e^x}-\mathrm{\infty }<x<\mathrm{\infty }.\end{array}$$

The population mean and the population variance are not mathematically tractable, but the numeric values, to ten digits, are

$$\begin{array}{l}E[X]=-0.5772156649\text{and}V[X]=1.644934067.\end{array}$$

The probability density function is not symmetric about the mean.

The R code below plots these three probability density functions on the same set of axes. The standard normal probability density function is taken directly from the formulas in the previous paragraph. The probability density functions for the standard logistic distribution and the standard extreme value distribution have been standardized (by subtracting their population mean and dividing by the population standard deviation) so that all three probability density functions can be viewed on an equal footing. The plot emphasizes the shape of the various probability density functions.

The results are displayed in Figure 3.29. All three probability distributions have support on the entire real number line [latex]-\infty < x < \infty[/latex], although the graph only includes the values within three standard deviation units from the population mean. As expected, the probability density functions for the standard normal distribution and the standardized version of the standard logistic distribution are symmetric and bell-shaped. The probability density function of the standardized version of the standard extreme value distribution is nonsymmetric. The R code below plots the cumulative distribution function associated with the standardized version of the standard logistic distribution.

The cumulative distribution function [latex]F(x) = P(X \le x)[/latex] is graphed in Figure 3.30. This cumulative distribution function is monotone increasing and satisfies [latex]\lim_{x \, \rightarrow \, - \infty} F(x) = 0[/latex] and [latex]\lim_{x \, \rightarrow \, \infty} F(x) = 1[/latex]. Notice that a plot of [latex]F(-x)[/latex] gives the complement of the cumulative distribution function. In other words, [latex]S(x) = 1 - F(x) = P(X \ge x)[/latex]. This function is monotone decreasing and satisfies [latex]\lim_{x \, \rightarrow \, - \infty} S(x) = 1[/latex] and [latex]\lim_{x \, \rightarrow \, \infty} S(x) = 0[/latex]. This function is known in survival analysis as the survivor function.

Now that cumulative distribution functions and their complements have been identified as a reasonable way to estimate the probability of success for the field goal data, we would like to establish a mechanism for incorporating the value of the predictor X into the probability model. The emphasis here will be on using the cumulative distribution function for the logistic distribution, since that seems to be the most commonly used in logistic regression.

The usual form of the mean response function for simple linear regression is

$$\begin{array}{l}E[Y]={\beta }_0+{\beta }_{1}X.\end{array}$$

But in the case of a binary outcome, the constraint

$$\begin{array}{l}0\le E[Y]\le 1\end{array}$$

must be imposed. This is done naturally using the cumulative distribution functions and their complements for the various probability distributions described earlier. Let [latex]\pi (X)[/latex] be the mean response function for a regression model with a binary response. Using the cumulative distribution function for the logistic distribution, the mean response function is

$$\begin{array}{l}\pi (X)=E[Y]=\frac{e^{{\beta }_0+{\beta }_{1}X}}{1+e^{{\beta }_0+{\beta }_{1}X}}.\end{array}$$

Since the random variable Y can only assume the values 0 and 1 for a particular value of X, it is a Bernoulli random variable with probability of success [latex]\pi(X)[/latex]. Since the expected value and the probability that a Bernoulli random variable assumes the value 1 are equal, the mean response function can also be expressed as

$$\begin{array}{l}\pi (X)=P(Y=1)=\frac{e^{{\beta }_0+{\beta }_{1}X}}{1+e^{{\beta }_0+{\beta }_{1}X}},\end{array}$$

where [latex]P(Y = 1)[/latex] is the probability that the dependent variable Y equals 1 for a particular fixed setting of the independent variable X. The parameters β₀ and β₁ assume the following roles.

• The sign of β₁ controls whether the mean response function is monotone increasing or decreasing. Table 3.9 shows the direction of the relationship associated with the sign of β₁. The statistical significance of the point estimator of β₁ depends on its magnitude.
• The magnitude of β₁ controls the steepness of the mean response function, with larger magnitudes corresponding to steeper mean response functions.
• The value of β₀ controls the location of the mean response function on the X-axis.

A graph that illustrates the effect of varying values of β₁ for the fixed value of [latex]\beta_0 = 0[/latex] on the mean response function [latex]\pi(X)[/latex] is given in Figure 3.31. As expected, the mean response function [latex]\pi(X)[/latex] is monotone decreasing for [latex]\beta_1 < 1[/latex] and monotone increasing for [latex]\beta_1 > 1[/latex]. The mean response function is steeper as the magnitude of β₁ increases.

A graph that illustrates the effect of varying values of β₀ for the fixed value of [latex]\beta_1 = 1[/latex] on the mean response function [latex]\pi(X)[/latex] is given in Figure 3.32. As expected, the mean response function [latex]\pi(X)[/latex] is monotone increasing in all cases because [latex]\beta_1 > 1[/latex]. The effect of varying β₀ is to shift the mean response functions horizontally. The rationale behind the horizontal shift can be seen by writing the mean response function with [latex]\beta_1 = 1[/latex] as

$$\begin{array}{l}\pi (X)=\frac{e^{{\beta }_0+X}}{1+e^{{\beta }_0+X}}.\end{array}$$

So the effect of increasing β₀ in this case is to shift the mean response function to the right (for [latex]\beta_1 < 1[/latex]) or to the left (for [latex]\beta_1 > 1[/latex]) relative to the [latex]\pi(X)[/latex] curve associated with [latex]\beta_0 = 0[/latex].

To summarize, the sign of β₁ controls the direction of the monotonicity of [latex]\pi(X)[/latex], the magnitude of β₁ controls the steepness of [latex]\pi(X)[/latex], and β₀ controls the location of [latex]\pi(X)[/latex] along the X-axis.

We now consider the estimation of the parameters β₀ and β₁ from a data set consisting of the n data pairs [latex](X_1, \, Y_1)[/latex], [latex](X_2, \, Y_2)[/latex], …, [latex](X_n, \, Y_n)[/latex]. The first components [latex]X_1, \, X_2, \, \ldots, \, X_n[/latex] are real numbers and the second components [latex]Y_1, \, Y_2, \, \ldots, \, Y_n[/latex] assume only the values 0 and 1. Since

$$\begin{array}{l}P(Y=1)=\pi (X)\text{and}P(Y=0)=1-\pi (X)\end{array}$$

the contribution to the likelihood function of the data pair [latex](X_i, \, Y_i)[/latex] is

$$\begin{array}{l}\pi (X_i{)}^{Y_i}{[1-\pi (X_i)]}^{1-Y_i}\end{array}$$

for [latex]i = 1, \, 2, \, \ldots, \, n[/latex]. When [latex]Y_i = 0[/latex], the contribution to the likelihood function is [latex]1 - \pi(X_i)[/latex], which is [latex]P(Y_i = 0)[/latex], where [latex]P(Y_i = 0)[/latex] is the probability that [latex]Y_i = 0[/latex] for the particular setting of the independent variable at X_i. When [latex]Y_i = 1[/latex], the contribution to the likelihood function is [latex]\pi(X_i)[/latex], which is [latex]P(Y_i = 1)[/latex]. Since X_i is assumed to be observed without error, Y_i is a random binary response, and the responses are assumed to be mutually independent random variables, the likelihood function is

$$\begin{array}{l}L({\beta }_0,{\beta }_1)=\prod _{i=1}^n\pi (X_i{)}^{Y_i}{[1-\pi (X_i)]}^{1-Y_i}.\end{array}$$

The log likelihood function is

$$\begin{array}{l}\ln L({\beta }_0,{\beta }_1)=\sum _{i=1}^n{Y}_i\ln [\pi (X_i)]+(1-Y_i)\ln [1-\pi (X_i)].\end{array}$$

This can be written in terms of β₀ and β₁ as

$$\begin{array}{l}\ln L({\beta }_0,{\beta }_1)=\sum _{i=1}^n{Y}_i[{\beta }_0+{\beta }_1{X}_i-\ln (1+e^{{\beta }_0+{\beta }_1{X}_i})]-(1-Y_i)\ln (1+e^{{\beta }_0+{\beta }_1{X}_i})\end{array}$$

or

$$\begin{array}{l}\ln L({\beta }_0,{\beta }_1)=\sum _{i=1}^n{Y}_i({\beta }_0+{\beta }_1{X}_i)-\ln (1+e^{{\beta }_0+{\beta }_1{X}_i}).\end{array}$$

The likelihood function and the log likelihood function are maximized at the same values of β₀ and β₁ because the natural logarithm is a monotonic transformation. The score vector is comprised of the partial derivatives of the log likelihood function with respect to β₀ and β₁:

$$\begin{array}{l}\frac{\partial \ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_0}=\sum _{i=1}^n(Y_i-\frac{e^{{\beta }_0+{\beta }_1{X}_i}}{1+e^{{\beta }_0+{\beta }_1{X}_i}})\end{array}$$

and

$$\begin{array}{l}\frac{\partial \ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_1}=\sum _{i=1}^n(X_i{Y}_i-\frac{X_i{e}^{{\beta }_0+{\beta }_1{X}_i}}{1+e^{{\beta }_0+{\beta }_1{X}_i}}).\end{array}$$

When these two equations are equated to zero, there is no closed form solution for [latex]\hat \beta_0[/latex] and [latex]\hat \beta_1[/latex], so numerical methods must be relied on to calculate these point estimates. The second derivatives of the log likelihood function after simplification are

$$\begin{array}{l}\frac{{\partial }^2\ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_0^2}=-\sum _{i=1}^n\frac{e^{{\beta }_0+{\beta }_1{X}_i}}{{(1+e^{{\beta }_0+{\beta }_1{X}_i})}^2},\end{array}$$

$$\begin{array}{l}\frac{{\partial }^2\ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_0\partial {\beta }_1}=-\sum _{i=1}^n\frac{X_i{e}^{{\beta }_0+{\beta }_1{X}_i}}{{(1+e^{{\beta }_0+{\beta }_1{X}_i})}^2},\end{array}$$

and

$$\begin{array}{l}\frac{{\partial }^2\ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_1^2}=-\sum _{i=1}^n\frac{X_i^2{e}^{{\beta }_0+{\beta }_1{X}_i}}{{(1+e^{{\beta }_0+{\beta }_1{X}_i})}^2}.\end{array}$$

The Fisher information matrix is the matrix of expected values of these partial derivatives:

$$\begin{array}{l}I({\beta }_0,{\beta }_1)=\left(\begin{array}{cc}E\left[\frac{-{\partial }^2\ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_0^2}\right]& E\left[\frac{-{\partial }^2\ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_0{\beta }_1}\right]\\ E\left[\frac{-{\partial }^2\ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_1{\beta }_0}\right]& E\left[\frac{-{\partial }^2\ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_1^2}\right]\end{array}\right).\end{array}$$

The expected values in this matrix can be determined because they do not contain any random variables. Their values cannot be calculated, however, because the values of the parameters β₀ and β₁ are unknown. The observed information matrix

$$\begin{array}{l}O({\hat{\beta }}_0,{\hat{\beta }}_1)={\left(\begin{array}{cc}\frac{-{\partial }^2\ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_0^2}& \frac{-{\partial }^2\ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_0{\beta }_1}\\ \frac{-{\partial }^2\ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_1{\beta }_0}& \frac{-{\partial }^2\ln L({\beta }_0,{\beta }_1)}{\partial {\beta }_1^2}\end{array}\right)}_{{\beta }_0={\hat{\beta }}_0,{\beta }_1={\hat{\beta }}_1}\end{array}$$

can be estimated from data values once the maximum likelihood estimators are computed. This matrix is the variance–covariance matrix of the score vector and its inverse is the asymptotic variance–covariance matrix of the maximum likelihood estimators. The square roots of the diagonal elements of this inverse matrix provide estimates of the standard errors of the maximum likelihood estimates.

The NFL field goal data set has a large sample size ([latex]n = 948[/latex]) and a strong statistical relationship between the length of the field goal attempt and the probability of success. The R code below again uses the optim function to calculate the parameter estimates. The first argument to optim are initial parameter estimates. The second argument to optim is the function to be minimized, so the negative of the log likelihood function is given as the second argument. Once the maximum likelihood estimates are calculated, the observed information matrix, standard errors, z-statistics, and associated p-values are calculated.

The results of the code are summarized in Table 3.10. The values of [latex]\hat \beta_0[/latex] and [latex]\hat \beta_1[/latex] are both statistically significant with p-values near zero. The observed information matrix for the NFL field goal data set