Generalized F-test

Posted on October 14, 2024 by Dan Ma

This is a discussion of the F-test from a regression perspective and not from a analysis of variance perspective, even though there is also an ANOVA table listing out the sums of squares resulting from linear regression. From a regression perspective, when we think of F-test, we usually think of the F-test that is to test the significance of the entire regression model, which is an “all or nothing” test. If the F-statistic is sufficiently large, then we know the entire linear regression model is useful in explaining the variability in the response, i.e., some of the predictors are significant and have an impact on the response. Otherwise, we conclude that the predictors as a group do not improve the prediction of the response over that obtained by regressing the response on the intercept coefficient $\beta_0$ (this is the intercept-only model). For convenience, we call the “all or nothing” F-test as the whole model F-test. What if we are interested in knowing whether a subset of the predictors is significant? In other words, what if the comparison is not between the full model and the intercept-only model but is between the full model and a smaller model that is larger than the intercept-only model? We propose that the whole model F-test can be tweaked to test the significance of a reduced model versus the full model.

See here for an application of the generalized F-test.

The Whole Model F-Test

First, we start with the more familiar concept of the whole model F-test. As mentioned above, this test is “all or nothing.” The “all” is the full model with all the predictors and the “nothing” is the intercept-only model. In the intercept-only linear regression model, there are no slope coefficients and there is only the intercept coefficient. With the intercept-only model, the model fitting is easy since the fitted value is simply the sample mean of the observed response values. With “all or nothing,” it is about choosing between these two propositions: either we are doing no better than using the sample mean response as the prediction (going with the intercept-only model) or the linear regression model with all predictors is useful for explaining the variation in the response (going with the full model). These two models are described as follows:

(1)

……….

$Y=\beta_0+\beta_1 X_1+\beta_2 X_2 + \cdots+\beta_p X_p+ \epsilon$

….

(1)………. $Y=\beta_0+ \epsilon$ …………………………………………(Intercept-Only Model)

The hypotheses to be tested are the following:

(2)

……….

$H_0:\beta_1=\beta_2=\cdots=\beta_p=0$

….

(2)………. $H_1: \beta_j \ne 0$ for at least one $j$

If the null hypothesis $H_0$ holds, then we opt for the the intercept-only model, i.e., none of the predictors $X_1,\cdots,X_p$ makes any difference. If there is sufficient evidence to reject the null hypothesis, then we are in favor of the full model. As we will see below, the null hypothesis corresponds to the smaller model. The evidence for or against the null hypothesis is summarized by the following F-statistic, a quantity found in the ANOVA table,

(3)

……….

$\displaystyle F=\frac{\text{Reg SS} / p}{\text{RSS} / (n-p-1)}=\frac{(\text{TSS}-\text{RSS}) / p}{\text{RSS} / (n-p-1)}$

where $n$ is the sample size. See here for a more in-depth discussion of the ANOVA table and the F-statistic. Under the null hypothesis $H_0$ , the F-statistic follows an F distribution with degrees of freedom $p$ and $n-p-1$ with expected value close to 1. An observed value of the F-statistic that is substantially larger than 1 provides evidence to support the rejection of the null hypothesis.

In (3), RSS is the residual sum of squares, which is the sum of the squares of the differences “observed response minus fitted response.” On the other hand, TSS is the total sum of squares, which is the sum of squares of the differences “observed response minus sample mean response.” RSS is the amount of variability that is left unexplained after performing the regression, while TSS is the amount of variability inherent in the response before the regression is performed. As a result, TSS – RSS is the amount of variability in the response that is explained by performing the linear regression.

RSS and TSS can be regarded as as measurements of error because the difference in each sum of squares can be regarded as an error. The difference “observed response minus fitted response” measures the error associated with a fitted value in the full model. On the other hand, the difference “observed response minus sample mean response” measures the error associated with a fitted value in the intercept-only model. Thus RSS is a sum of squares that measures the error with respect to the full model (the smaller the RSS, the better the full model fits the data). In the terminology used here, both RSS and TSS are called errors sums of squares with respect to the appropriate models (for RSS, it’s the full model and for TSS, it’s the intercept-only model).

The model with more predictors will always be able to fit the data at least as well as the model with fewer predictors. Thus, the full model will give a better fit to the data than the intercept-only model, i.e., RSS $\le$ TSS. The hypothesis test described in (2) is to determine whether the full model will provide a significantly better fit to the data than the intercept-only model, i.e., whether RSS will be significantly lower than TSS. This is accomplished by the F-statistic given in (3).

Tweaking the Whole Model F-Test

The full model in (1) is still called the full model. The intercept-only model is now called the reduced model. It is called the reduced model because we drop the predictors $X_1,\cdots,X_p$ from the full model. The F-statistic in (3) is re-written as follows.

(4)

……….

$\displaystyle F=\frac{(\text{RSS}_0-\text{RSS}_1) / p}{\text{RSS}_1 / (n-p-1)}$

RSS with subscript 0 is error sum of squares of the reduced model, which is the TSS in (3). RSS with subscript 1 is the error sum of squares of the full model, which is the RSS in (3). In (4), $p$ is the number predictors dropped from the full model. The null hypothesis tested by the F-statistic in (4) is still the same as the one in (2), but it can be restated that the full model does not provide a significantly better fit than the reduced model. Under the null hypothesis, the F-statistic (4) will follow an F-distribution with degrees of freedom $p$ and $n-p-1$ .

The Generalized F-Test

With (4) as preparation, we are now ready to expand the F-test to handle “full versus reduced” and not just “full versus intercept-only.” First, we describe precisely the two models in question.

(5)

…

$Y=\beta_0+\beta_1 X_1+ \cdots + \beta_k X_k + \beta_{k+1} X_{k+1}+ \cdots+\beta_{k+r} X_{k+r}+ \epsilon$

(5)

…

(5)… $Y=\beta_0+\beta_1 X_1+ \cdots + \beta_k X_k+ \epsilon$
(5)…(Reduced Model)

We let $p=k+r$ so that the full model in (5) has the same number of predictors as the full model in (1). Here, we drop $r$ many predictors from the full model to form the reduced model. For clarity, we arrange the dropped predictors in the end of the model equation of the full model. The predictors in the reduced model are $X_1,X_2,\cdots,X_k$ . The predictors in the full model but not in the reduced model are $X_{k+1},X_{k+2},\cdots,X_{k+r}$ . The following are the hypotheses to be tested.

(6)

……….

$H_0:\beta_{k+1}=\beta_{k+2}=\cdots=\beta_{k+r}=0$

….

(6)………. $H_1: \beta_{k+j} \ne 0$ for at least one $j=1,2,\cdots,r$

The coefficients identified in the null hypothesis $H_0$ are the ones for the predictors dropped from the full model. Thus, $H_0$ corresponds to the reduced model. Thus, not rejecting the null hypothesis means that we choose the reduced model. On the other hand, rejecting the null hypothesis means that the predictors $X_{k+1},X_{k+2},\cdots,X_{k+r}$ as a whole should not be dropped from the full model, i.e., one of these predictors plays an important role in explaining the variability in the response. We use the following F-statistic to test the hypotheses in (6).

(7)

……….

$\displaystyle F=\frac{(\text{RSS}_0-\text{RSS}_1) / r}{\text{RSS}_1 / (n-p-1)}$

where $n$ is the sample size, $p$ is the number of predictors in the full model with $p=k+r$ and $r$ is the number of predictors dropped from the full model. Furthermore, RSS with subscript 0 is the RSS (the error sum of squares) of the reduced model and RSS with subscript 1 is the RSS (error sum of squares) of the full model. The following clarifies these error sums of squares.

$\displaystyle \text{RSS}_0$

$\displaystyle \text{RSS}_1$

The difference $\text{RSS}_0-\text{RSS}_1$ is the extra sum of squares. As noted above, having more predictors in the model lowers its RSS. As a result the difference $\text{RSS}_0-\text{RSS}_1 \ge 0$ . Here’s what we can say about the extra error sum of squares

$\displaystyle \text{RSS}_0-\text{RSS}_1$

A large extra sum of squares $\displaystyle \text{RSS}_0-\text{RSS}_1$ indicates that the predictors not found in the reduced model (or the predictors dropped from the full model) improve the prediction of the response over that obtained by regressing the response on the predictors in the reduced model. It follows from the preceding paragraph that the extra sum of squares $\text{RSS}_0-\text{RSS}_1$ can be used as a basis for testing the importance of the predictors $X_{k+1},X_{k+2},\cdots,X_{k+r}$ in the presence of $X_{1},X_{2},\cdots,X_{k}$ . In other words, $\text{RSS}_0-\text{RSS}_1$ can help us determine whether the additional predictors $X_{k+1},X_{k+2},\cdots,X_{k+r}$ can add significant predictive power over the predictors $X_{1},X_{2},\cdots,X_{k}$ that are already in the reduced model.

Thus, the test statistic for testing the null hypothesis in (6) should be a function of the extra sum of squares $\text{RSS}_0-\text{RSS}_1$ , with the direction that the larger the extra sum of squares, the larger the test statistic. The F-statistic in (7) is one such test statistics. It can be used to judge whether the extra sum of squares $\text{RSS}_0-\text{RSS}_1$ is large enough to warrant the rejection of the null hypothesis.

Under the null hypothesis $H_0$ in (6), it can be shown that the F-statistic in (7) follows an F-distribution with degrees of freedom $r$ and $n-p-1$ . How large does the extra sum of squares $\text{RSS}_0-\text{RSS}_1$ for us to be willing to reject the null hypothesis? The answer lies in evaluating the magnitude of the F-statistic, as described in the algorithm below.

The Algorithm

The denominator of the F-statistic in (7) is the mean squared error (MSE) of the full model. To calculate the F-statistic, we need to fit two regression models.

Run the full model as indicated by the model equation in (5). Obtain the error sum of squares $\text{RSS}_1$ and the MSE, i.e., $\text{RSS}_1 / (n-p-1)$ .
Run the reduced model as indicated by the model equation in (5). Obtain the error sum of squares $\text{RSS}_0$ .
Compute the F-statistic according to (7).
Calculate the p-value under the assumption that the null hypothesis $H_0$ in (6) is true, calculate the p-value under the assumption that the F-statistic in (7) follows an F-distribution with degrees of freedom $r$ and $n-p-1$ . Let $G(x)$ denote the cumulative distribution for this F-distribution. Recall that the p-value is the probability of obtaining an observed value of the F-statistic more extreme than the on observed from the data. It follows that the p-value is $1-G(F)$ where $F$ is the observed value of the F-statistic.
If the p-value is sufficiently small (e.g., 0.05 or 0.01), the F-statistic is large enough to warrant the rejection of the null hypothesis. We can then conclude that we should not drop the predictors $X_{k+1},X_{k+2},\cdots,X_{k+r}$ from the full model. The small p-value provides evidence that these predictors play an important role in explaining the variability in the response.
If the p-value is not below a pre-specified threshold (e.g., 0.05 or 0.01), the value of the F-statistic is not large enough to warrant the rejection of the null hypothesis. Adding the predictors $X_{k+1},X_{k+2},\cdots,X_{k+r}$ to the reduced model does not add sufficient predictive power to the model. In other words, we conclude that the additional predictor variables not in the reduced model do not improve the prediction of the response over that obtained by regressing the response on the predictors in the reduced model.

ANOVA Tables

As indicated in the preceding section, the generalized F-test requires fitting two separate linear regression models, one for the full model and one for the reduced model. We can use ANOVA tables to organize the information in such an F-test, one table for the full model and one for the reduced model.

ANOVA Table – Full Model

Source of Variation….	Sum of Squares…..	df……………..	Mean Square.
$\text{Regression}$	$\text{TSS}-\text{RSS}_1$	$p$	$(\text{TSS}-\text{RSS}_1) /p$
$\text{Error}$	$\text{RSS}_1$	$n-p-1$	$s^2=\text{RSS}_1/(n-p-1)$
$\text{Total}$	$\text{TSS}$	$n-1$

ANOVA Table – Reduced Model

Source of Variation….	Sum of Squares…..	df……………..	Mean Square.
$\text{Regression}$	$\text{TSS}-\text{RSS}_0$	$k$	$(\text{TSS}-\text{RSS}_0) /k$
$\text{Error}$	$\text{RSS}_0$	$n-k-1$	$s^2=\text{RSS}_0/(n-k-1)$
$\text{Total}$	$\text{TSS}$	$n-1$

In each of the tables, only the items enclosed with red borders are needed for the F-test. In the table for the full model, we need the error sum of squares (the RSS with subscript 1) and its corresponding mean square. In the table for the reduced model, we need the error sum of squares, the RSS with subscript 0. Then we can compute the F-statistic according to (7).

Example 1

We work some examples using the Boston dataset from the MASS package in R. The Boston dataset contains housing data in Boston MA and the surrounding area that were collected by the US Census. It is a relatively small dataset with only 506 records. Among the 14 variables in the dataset, the response variable is medv, the median value of owner-occupied homes in $1000’s. See here for more on the dataset and the R code. The following is the linear regression that is the full model (with 13 predictors).

Summary of MLR – All Predictors in Boston Dataset

All the predictors are statistically significant with the except of indus and age where indus is the proportion of non-retail business acres per town and age is the proportion of owner-occupied units built prior to 1940. Let’s drop these 2 predictors from the full model. We then fit a linear model with 11 predictors (the 13 predictors from the full model excluding indus and age). The ANOVA tables are shown below.

Example 1 – ANOVA Table – Full Model – Boston Dataset

Source of Variation….	Sum of Squares…..	df……………..	Mean Square.
Regression	31,637.51	13	2,433.655
Error	11,078.78	492	22.51785
Total	42,716.3	505

Example 1 – ANOVA Table – Reduced Model – Boston Dataset

Source of Variation….	Sum of Squares…..	df……………..	Mean Square.
Regression	31,634.93	11	2875.903
Error	11,081.36	494	22.43191
Total	42,716.3	505

Using the numbers enclosed by the red borders, the following is the calculation of the F-statistic.

(7)

……….

$\displaystyle F=\frac{(11,081.36-11,078.78) / 2}{22.51785}=0.057$

The value of the F-statistic is very small (close to 0)! This means that the extra sum of squares is too small. Dropping the two predictors indus and age does not deteriorate the model fit in any noticeable way. Observe that the two errors sums of squares are very close together (11,078.78 and 11,081.36). As a result, we will not reject the null hypothesis that the two dropped variables have predictive power. This implies that we should use the reduced model instead of the full model, i.e., it is OK to drop the two variables in question. In other words, we can conclude that the additional predictor variables not found in the reduced model do not improve the prediction of the response over that obtained by regressing the response on the predictors in the reduced model.

In this case, the value of the F-statistic is so small that there is no need to obtain the p-value. Just for the record, the p-value is 0.9819827, which obviously does not support the rejection of the null hypothesis.

Example 2

Example 1 shows that the predictors indus and age can be safely dropped without losing much predictive power. Now we use the reduced model in Example 1 as the full model. For the reduced model, we drop the predictor chas, which is the Charles River dummy variable (1 if tract bounds river; 0 otherwise). We use the F-test to determine whether chas is an important predictor of median house price. The two regression models are displayed in the following two ANOVA tables.

Example 2 – ANOVA Table – Full Model – Boston Dataset

Source of Variation….	Sum of Squares…..	df……………..	Mean Square.
Regression	31,634.93	11	2875.903
Error	11,081.36	494	22.43191
Total	42,716.3	505

Example 2 – ANOVA Table – Reduced Model – Boston Dataset

Source of Variation….	Sum of Squares…..	df……………..	Mean Square.
Regression	31,407.72	10	3140.772
Error	11308.58	495	22.84561
Total	42,716.3	505

Using the numbers enclosed by the red borders, the following is the calculation of the F-statistic.

(7)

……….

$\displaystyle F=\frac{(11308.58-11,081.36) / 1}{22.43191}=10.129$

The F-statistic appears to be substantially larger than 1. The p-value is 0.001551469, which is indeed small. The small p-value provides evidence to support the rejection of the null hypothesis. This implies that the predictor chas does indeed posses substantial predictive power even in the presence of all the predictors and should not be dropped from the full model. Thus, we can conclude that the additional predictor variable chas does improve the prediction of the response over that obtained by regressing the response on the predictors in the reduced model.

t-Test versus F-Test

The above discussion casts the F-test as a comparison of two nested models, i.e., the full model with all the predictors and a smaller model with a subset of the predictors. The t-test can also be viewed in such a framework. The t-test is to decide on these on these hypotheses (see here).

(8)

……….

$H_0:\beta_{j}=0$

….

(8)………. $H_1: \beta_{j} \ne 0$ ….where $j$ is fixed with $j=1,2,\cdots,p$

The null hypothesis in (8) is evaluated using the t-statistic as described below.

(9)

……….

$\displaystyle t=\frac{\hat{\beta_j}}{\text{Standard Error of } \hat{\beta_j}}$

Under the null hypothesis in (8), the t-statistic follows a t-distribution with degrees of freedom $n-p-1$ . The magnitude of the t-statistic is used to evaluate the significance of the coefficient $\beta_j$ (see here). Even though this test is to test the significance of one parameter, the test is performed in the presence of all the other parameters (or in the presence of all the other predictors). In other words, this test is to determine the significance of one predictor in the presence of all other predictors. This is a test of the full model versus the reduced model with one and only one predictor missing. Thus, the t-test as described in (8) and (9) is to choose one of the following models.

(10)

……….

$Y=\beta_0+\beta_1 X_1+\cdots+\beta_{j-1} X_{j-1}+\beta_j X_j+\beta_{j+1} X_{j+1} + \cdots+\beta_p X_p+ \epsilon$

…..

(10)………. $Y=\beta_0+\beta_1 X_1+\cdots+\beta_{j-1} X_{j-1}+\beta_{j+1} X_{j+1} + \cdots+\beta_p X_p+ \epsilon$ …………..(Reduced Model)

The t-test as described in (8) and (9) is equivalent to the F-test with the full model and reduced model described in (10). When calculating the F-statistic according to (7), $r=1$ since we are only dropping one predictor from the full model. In fact, the F-statistic in this case is the square of the t-statistic in (9). In Example 2, the reduced model is the result of dropping the variable chas. The value of the F-statistic is 10.129. As a result, the t-statistic is 3.183. The corresponding p-value is 0.001551, which is identical to the p-value if F-test is used.

To interpret the t-test, as in the F-test, we must do so in the presence of the other predictor. If the null hypothesis $H_0:\beta_j=0$ is rejected, we conclude that the inclusion of the predictor $X_j$ in the model will significantly improves the prediction of the response over that obtained by regressing on all the predictors excluding $X_j$ . If $H_0:\beta_j=0$ is not rejected, we conclude that including the predictor $X_j$ does not significantly improve the prediction of the response over that obtained by regressing the response on all the predictors excluding $X_j$ .

The F-Statistic in Another Form

The F-statistic in (7) is defined using the error sum of squares in two models. It can be restated using the $\text{R}^2$ statistics of the full model and reduced model.

(11)

……….

$\displaystyle F=\frac{(\text{R}_1^2-\text{R}_0^2) / r}{(1-\text{R}_1^2) / (n-p-1)}$

where $\text{R}_1^2$ is the $\text{R}^2$ statistic for the full model and $\text{R}_0^2$ is the $\text{R}^2$ statistic for the reduced model. (11) is derived by dividing the numerator and denominator of (7) by the TSS. Adding more predictors to a model always increases the $\text{R}^2$ statistic. Thus, $\text{R}_1^2 \ge \text{R}_0^2$ . The F-statistic (11) will help us determine whether the increase in $R^2$ is significant when adding predictors into the reduced model.

……….

This article was originally published here in a companion site.

$\text{ }$

Dan Ma Generalized F-Test
Daniel Ma Generalized F-Test

Dan Ma F-Test
Daniel Ma F-Test

Dan Ma full model versus reduced model
Daniel Ma full model versus reduced model

Dan Ma linear regression
Daniel Ma linear regression

Dan Ma simple linear regression
Daniel Ma simple linear regression

Dan Ma multiple linear regression
Daniel Ma multiple linear regression

$\copyright$ 2024 – Dan Ma

Posted: October 14, 2023

Ordinary least squares estimators are BLUE

Posted on October 13, 2024 by Dan Ma

The Gauss-Markov Theorem, which states that the least squares estimator in linear regression has the smallest variance among the unbiased estimators of the unknown coefficients that are linear combinations of the response $Y_i$ ‘s whenever some assumptions about the random error terms are met. As a result, the least squares estimator is BLUE (Best Linear Unbiased Estimator). We do not aim to prove the theorem. The goal is to unpack this theorem and put it in a proper context.

The Linear Regression Model

We first describe linear regression. In linear regression, the following model equation relates the response $Y$ to the predictors $X_1,X_2,\cdots,X_p$ ,

(1)

……….

$Y=\beta_0+\beta_1 X_1+\beta_2 X_2 + \cdots+\beta_p X_p+ \epsilon$

where $p$ is the number of predictors, $\beta_0$ is the intercept coefficient, $\beta_j$ is the model coefficient for the predictor $X_j$ with $j=1,2,\cdots,p$ , and $\epsilon$ is the random error term. The random error term $\epsilon$ is normally distributed with mean 0 and constant variance $\sigma^2$ .

In a training sample of $n$ observations, the $i$ th observation is of the form $(Y_i,X_{i1},X_{i2},\cdots,X_{ip})$ . As a result, there are $n$ model equations, one for each observation, that can be expressed as follows.

(2)

……….

$Y_i=\beta_0+\beta_1 X_{i1}+\beta_2 Y_{i2} + \cdots+\beta_p X_{ip}+ \epsilon_i$

(2)……….where $i=1,2,3,\cdots,n$

The above $n$ model equations can be better described in matrix notation.

(3)

……….

$\left[\begin{array}{c} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{array}\right] = \left[\begin{array}{ccccc} 1 & X_{11} & X_{12} & \cdots & X_{1p }\\ 1 & X_{21} & X_{22} & \cdots & X_{2p } \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & X_{n1} & X_{n2} & \cdots & X_{np } \end{array}\right] \left[\begin{array}{c} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{array}\right] + \left[\begin{array}{c} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{array}\right]$

(.)

……….

$\bold Y= \bold X \bold \beta + \bold \epsilon$

where

$\bold Y$ is the $n \times 1$ vector of response variables,
$\bold X$ is the matrix containing the values of the predictor variables, an $n \times (p+1)$ matrix (the matrix $\bold X$ is also known as the data matrix or design matrix),
$\beta$ is the $(p+1) \times 1$ vector of model coefficients,
$\epsilon$ is the $n \times 1$ vector of random error terms.

The least squares estimates are so named because we estimate the unknown coefficients indicated in (1) by choosing $\beta_0,\beta_1,\cdots,\beta_n$ such that the following sum of the squared differences is minimized.

(4)

……….

$\displaystyle \sum_{i=1}^n \biggl[Y_i-(\beta_0+\beta_1 X_{i1}+\beta_2 X_{i2}+ \cdots + \beta_p X_{ip}) \biggr]^2$

Minimizing (4) results in a system of $p+1$ simultaneous linear equations in $p+1$ unknowns, which are $\beta_0,\beta_1,\cdots,\beta_p$ . These equations are known as the normal equations. The normal equations can be compactly expressed in matrix form as follows.

(5)

……….

$\displaystyle ( \bold X^{'} \bold X) \hat{\beta}=\bold X^{'} \bold Y$

where the notation “‘” (read prime) refers to the transpose of a matrix. It can be shown that the solution to these equation will minimize (4). Provided that the inverse of the matrix $\bold X^{'} \bold X$ exists, the least squares estimates of the model coefficients $\beta_0,\beta_1,\cdots,\beta_p$ can be obtained by the following.

(6)

……….

$\displaystyle \hat{\beta}=( \bold X^{'} \bold X)^{-1} \bold X^{'} \bold Y$

The matrix calculation in (6) is a convenient way to produce the least squares estimates for any multiple linear regression model. Once the design matrix $\bold X$ is specified and once the vector $\bold Y$ is specified, the least squares estimates can be determined according to the matrix multiplications indicated in (6).

Note that (6) tells us that the least squares estimator $\hat{\beta}$ are linear combinations of the response variables $Y_i$ ‘s. The following makes it explicit.

(7)

……….

$\displaystyle \hat{\beta_k}=c_{k1} Y_1+c_{k2} Y_2+\cdots+c_{kn} Y_n$

……….

(7)

……….

…

$k=0,1,2,\cdots,p$

and that $(c_{k1},c_{k2},\cdots,c_{kn})$ is the $(k+1)$ st row of the matrix $( \bold X^{'} \bold X)^{-1} \bold X^{'}$ , which is a $(p+1) \times n$ matrix.

Furthermore, the estimator $\hat{\beta}$ is an unbiased estimator, which means that $E(\hat{\beta})=\beta$ (for the derivation, see here). To express this component-wise, we have $E(\hat{\beta_k})=\beta_k$ where $k=0,1,2,\cdots,p$ . In other words, the least squares estimator is correct on average. The expected value of the vector $\hat{\beta}$ of the least squares estimates is the vector $\hat{\beta}$ of the regression model coefficients. Thus, we have observed two properties of the least squares estimator, which are stated below.

(8)

……….

The least squares estimator as described in (6) is unbiased.
The least squares estimator as described in (6) are linear combinations of the response variables $Y_i$ ‘s.

Gauss-Markov Theorem

There are numerous ways to estimate the unknown coefficients in the linear model equation (1). The most common way is the least squares approach, which is to minimize the sum of squared errors described in (4). Another way to estimate the unknown model coefficients is to minimize the sum of 4th powers of errors (i.e., changing the exponent in (4) from 2 to 4). However, it is not clear if minimizing the sum of errors to the 4th power will result in an estimator that is unbiased and linear. With respect to the Gauss-Markov theorem, we are interested in all unbiased linear estimates of the unknown coefficients indicated in the model equation (1). According to (8), the least squares estimator is both unbiased and linear. The Gauss-Markov theorem says that out of all the unbiased linear estimators of the unknown model coefficients in (1), the least squares estimator has the lowest variance. Linear estimator means that the estimator is a linear combination of the response $Y_i$ ‘s. The following is the statement of the theorem.

Gauss-Markov Theorem
Provided that the random error terms $\epsilon_i$ in the model equations (2) have means zero and constant variance $\sigma^2$ and are uncorrelated, the least squares estimator described in (6) have the smallest variance among all unbiased linear estimators of the unknown model coefficients $\beta_0,\beta_1,\cdots,\beta_p$ . Here, a linear estimator is a linear combination of the response variables $Y_i$ ‘s.

Interestingly, the random error terms $\epsilon_i$ are not required to be normally distributed for the theorem to work. However, they are required to be normally distributed in order to perform other linear regression model processes such as inference and predictions (see here). Furthermore, if the normality assumption of the error terms holds, then it can be shown that the least squares estimator is also the maximum likelihood estimator, providing another justification of using least squares linear regression.

In multiple linear regression, the variance of an estimator is a variance covariance matrix. For example, the variance-covariance matrix of the least squares estimator $\hat{\beta}$ is given by the following.

(9).. $\displaystyle \text{Var}(\hat{\beta})=\left[\begin{array}{ccccc} \text{Var}(\hat{\beta_0}) & \text{Cov}(\hat{\beta_0},\hat{\beta_1}) & \text{Cov}(\hat{\beta_0},\hat{\beta_2}) & \cdots & \text{Cov}(\hat{\beta_0},\hat{\beta_p})\\ \text{Cov}(\hat{\beta_1},\hat{\beta_0}) & \text{Var}(\hat{\beta_1}) & \text{Cov}(\hat{\beta_1},\hat{\beta_2}) & \cdots & \text{Cov}(\hat{\beta_1},\hat{\beta_p}) \\ \text{Cov}(\hat{\beta_2},\hat{\beta_0}) & \text{Cov}(\hat{\beta_2},\hat{\beta_1}) & \text{Var}(\hat{\beta_2}) & \cdots & \text{Cov}(\hat{\beta_2},\hat{\beta_p}) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \text{Cov}(\hat{\beta_p},\hat{\beta_0}) & \text{Cov}(\hat{\beta_p},\hat{\beta_1}) & \text{Cov}(\hat{\beta_p},\hat{\beta_2}) & \cdots & \text{Var}(\hat{\beta_p}) \end{array}\right]$

which happens to be $\sigma^2 (\bold X^{'} \bold X)^{-1}$ . In the above variance-covariance matrix, each diagonal element is the variance of an individual estimator in the vector $\hat{\beta}$ and each off diagonal element is the covariance of a pair of individual estimators. Suppose $\tilde{\beta}$ is another linear unbiased estimator of the unknown model coefficients $\beta_0,\beta_1,\cdots,\beta_p$ . The same can be said about the variance-covariance matrix $\text{Var}(\tilde{\beta})$ , that is, each diagonal element is the variance of $\tilde{\beta_j}$ for some $j$ and each off diagonal element is the covariance of a pair of individual estimators $\tilde{\beta_i}$ and $\tilde{\beta_j}$ . According to Gauss-Markov Theorem, $\text{Var}(\hat{\beta})$ is no more than $\text{Var}(\tilde{\beta})$ , which means that $\text{Var}(\hat{\beta})-\text{Var}(\tilde{\beta})$ is a non-negative definite matrix. In a non-negative definite matrix, each diagonal element is non-negative. The fact that $\text{Var}(\hat{\beta})-\text{Var}(\tilde{\beta})$ is a non-negative matrix boils down to the fact that $\text{Var}(\hat{\beta_j}) \le \text{Var}(\tilde{\beta_j})$ for all $j=0,1,2,\cdots,p$ .

Remarks

The Gauss-Markov Theorem states that least squares estimator is the best among all the linear unbiased estimators of the coefficients in model equation (1). One conclusion is clear. If the plan is to use a linear unbiased estimator to estimate the model described in (1), we should look no further and use the least squares estimator. According to the theorem, least squares regression is the best. Why use anything else? Least squares estimator is the best only among the linear unbiased estimator. It may not be the best among the non-linear estimators of the coefficients in (1). On the other hand, if we are willing to sacrifice a little on the bias, we could potentially obtain an estimator with smaller variance (e.g., by using ridge regression).

……….

This article was originally published here in a companion site.

$\text{ }$

Dan Ma Gauss-Markov Theorem
Daniel Ma Gauss-Markov Theorem

Dan Ma BLUE
Daniel Ma BLUE

Dan Ma Best Linear Unbiased Estimator
Daniel Ma Best Linear Unbiased Estimator

Dan Ma linear regression
Daniel Ma linear regression

Dan Ma simple linear regression
Daniel Ma simple linear regression

Dan Ma multiple linear regression
Daniel Ma multiple linear regression

$\copyright$ 2024 – Dan Ma

Posted: October 13, 2024

Statistical significance in SLR but not in MLR

Posted on October 12, 2024 by Dan Ma

An effect may be statistically significant in simple linear regression (SLR) but is not necessarily significant in multiple linear regression (MLR). Simple linear regression models are obviously different from multiple linear regression models in that SLR is to describe how the response is related to only one particular predictor while MLR is to describe the relationship between the response and multiple predictors. The implication of this difference is definitely more profound than merely counting the number of variables in the model equations. For example, it is possible that an SLR model and an MLR model fitted on the same dataset can reach different conclusions on the same predictor. It is possible that a predictor is deemed statistically significant in the SLR setting while the same predictor is insignificant in the MLR setting. Though this may seem like a contradiction, such phenomenon is not a contradiction. We show how to interpret such seemingly contradictory results correctly, using the Boston dataset from the MASS package in R. The Boston dataset is used here to illustrate the notion of partial correlation, which can provide insight into the difference between SLR and MLR.

The Boston Dataset

The Boston dataset contains housing data in Boston MA and the surrounding area that were collected by the US Census. It is a relatively small dataset with only 506 records. Among the 14 variables in the dataset, the variable of interest to us is medv, the median value of owner-occupied homes in $1000’s. The following is a description of the 14 variables in the Boston dataset.

crim – per capita crime rate by town.
zn – proportion of residential land zoned for lots over 25,000 sq.ft.
indus – proportion of non-retail business acres per town.
chas – Charles River dummy variable (1 if tract bounds river; 0 otherwise).
nox – nitric oxides concentration (parts per 10 million).
rm – average number of rooms per dwelling.
age – proportion of owner-occupied units built prior to 1940.
dis – weighted distances to five Boston employment centres.
rad – index of accessibility to radial highways.
tax – full-value property-tax rate per $10,000.
ptratio – pupil-teacher ratio by town.
black – 1000(Bk – 0.63)^2 where Bk is the proportion of blacks by town.
lstat – % lower status of the population.
medv – Median value of owner-occupied homes in $1000’s.

We would like to predict medv using the information contained in the other 13 variables. All the model results and graphical displays presented here are produced in here.

The Fitted Regression Models

First, let’s describe the discrepancy between SLR and MLR. Four SLR models are fitted by regressing medv on each of these four predictors indus, age, dis and rad, The results are shown below.

Summary of SLRs – Boston Dataset

All four predictors are highly significant with very large t-statistics (in absolute value) and small p-values that are indistinguishable from zero. These results seem to provide evidence that each of these 4 variables is an important predictor explaining the median house price. Now consider the MLR setting. The following shows the results of the MLR model obtained by regressing medv on all other 13 variables as predictors.

Summary of MLR – All Predictors in Boston Dataset

Of the four variables, dis and rad are still significant in the MLR model with large t-statistics (in absolute value) and small p-values. However, indus and age are no longer significant. The t-statistic for indus is only 0.334 and the p-value is 0.74. The t-statistic for age is even smaller (0.052) with a p-value of 0.96.

The t-statistic in the MLR results shown above is to test the null hypothesis, which is that the unknown model coefficient of the predictor in question is zero. Let’s consider the predictor age. Let $\beta_{\text{age}}$ be the unknown regression coefficient for the predictor age. The t-statistic for age in an MLR model is to test the null hypothesis $H_0:\beta_{\text{age}}=0$ . The t-statistic is 0.052. Since it is small, the t-statistic does not provide sufficient evidence to reject the null hypothesis. Then we would conclude that including the predictor age does not significantly improve the prediction of median house price versus the model that utilizes all the predictors except age. Therefore, conclusion of the significance of one predictor is conditional on all the other predictors in the MLR model. The inference on the null hypothesis $H_0:\beta_{\text{age}}=0$ is conditional on all other coefficients in the model. The hypothesis $H_0:\beta_{\text{age}}=0$ provides information about whether age is related to the median house price, after adjusting for the other predictors. Thus, the t-statistic in an MLR model cannot be interpreted in isolation and must be interpreted in the presence of all other predictors.

It follows from the preceding paragraph that there is no contradiction between the SLR result and the MLR result. Let’s take age as an example again. In the SLR model with age as the predictor, age is an important predictor of median house price in the absence of the other predictors. In the MLR model displayed above, the predictor age offers very little information to explain median house price in the presence of the other predictors. The key fact is that the MLR model has accounted for the effects of the other predictors.

Looking at Correlation

We now look at the phenomenon from the perspective of correlation. Again, let’s use age as an example. In the SLR model, age is an important predictor of the median house price. However, the SLR model has not taken the other predictors into account. For example, there is a strong correlation between age and the variable lsat, which is 0.6023 according to the following correlation matrix.

Correlation Matrix – Boston Dataset

The correlation of 0.6023 is moderate but is large enough that it cannot be ignored. Thus, age and lstat moves together. The higher the value of lstat (the higher the percent of the lower status in the population), the higher the value of age (the higher the proportion of owner-occupied units built prior to 1940). The predictor lstat is a strongly significant predictor of median house price according to the MLR model. The MLR model shows that the higher the value of lstat, the smaller the median house price. Thus, the effect of age in the SLR model is a reflection of the effect of lstat and other predictors. Once the other predictors are accounted for in the MLR model, the predictor age no longer has impact on median house price.

Looking at Partial Correlation

Let’s continue to use age as an example. The partial correlation between the response medv and age is the correlation between medv and age after the effects of the other predictors have been filtered out. The correlation between medv and age is -0.3770. The partial correlation between medv and age is 0.002347 (see here on how partial correlation is calculated). Thus before accounting for the other predictors, medv and age are negatively correlated at a moderate level. After the effects of the other predictors have been filtered out, medv and age are virtually uncorrelated since 0.002347 is virtually zero. The near zero partial correlation quantifies the fact that age is no longer relevant once all the other predictors are accounted for in the MLR model. The following is a side-by-side comparison of the scatterplot displaying the correlation and the added variable plot displaying the parital correlation.

Scatterplot and Added Variable Plot

The plot on the left is the scatterplot of the response medv against age. It shows a negative association, a visual conformation of -0.3770. The plot on the right is the added variable plot, which is a scatterplot of the two sets of residuals that have the effects of the other predictors removed. The partial correlation is defined as the correlation between these two sets of residuals. Any pattern shown in this scatterplot would indicate the true nature of the relationship between medv against age. This added variable plot shows no discernible pattern. The trend line is essentially the horizontal axis, a visually conformation that the partial correlation is near zero.

The following is the side-by-side comparison of the scatterplot and the added variable plot for medv and idus. The added variable plot on the right is also a visual confirmation that the partial correlation is near zero after accounting for the other predictors.

Scatterplot and Added Variable Plot

The following gives two more side-by-side comparisons of scatterplots and added variable plots. These two predictors are significant in both SLR and MLR. However, there is a difference between SLR and MLR in that the association is reversed. In the comparison for the predictor dis, the response medv and dis are positively associated, but are negatively associated once the other predictors are accounted for. In the comparison for the predictor rad, the reversal is the opposite (from negative to positive).

Scatterplot and Added Variable Plot

Remarks

Examples of significance in SLR but not in MLR is quite common. When there are multiple predictors in the models, correlation among the predictors can mask the real influence on the response. MLR results for one predictor must be interpreted in the presence of all the other predictors. Feature selection (or variable selection) is also important for identifying predictors that have predictive power.

……….

This article was originally published here in a companion site.

$\text{ }$

Dan Ma added variable plot/font>
Daniel Ma added variable plot

Dan Ma partial correlation
Daniel Ma partial correlation

Dan Ma linear regression
Daniel Ma linear regression

Dan Ma simple linear regression
Daniel Ma simple linear regression

Dan Ma multiple linear regression
Daniel Ma multiple linear regression

$\copyright$ 2024 – Dan Ma

Posted: October 12, 2024

Multiple linear regression – partial correlation

Posted on October 11, 2024 by Dan Ma

We describe a statistical phenomenon where the association between two variables emerges, disappears or reverses when the effects of other predictors are accounted for. Suppose that we are investigating how a numeric response variable is related to several numeric independent variables (or predictors). Suppose we find that one of the predictor, say $X_1$ , is positively associated with the response $Y$ as indicated by their correlation. This tells us that $X_1$ and $Y$ move together, e.g., when one is large, the other is large. However, it is possible that the positive relationship between $X_1$ and $Y$ is due to one or more additional variables that cause both $X_1$ and $Y$ to be large. When these additional variables are added into the linear regression model, the positive association reverses and becomes a negative association. We use the Boston dataset from the MASS package in R to illustrate this phenomenon. Examining the notion of partial correlation will also shed light on this seemingly anomalous data phenomenon.

The Boston Dataset

crim – per capita crime rate by town.
zn – proportion of residential land zoned for lots over 25,000 sq.ft.
indus – proportion of non-retail business acres per town.
chas – Charles River dummy variable (1 if tract bounds river; 0 otherwise).
nox – nitric oxides concentration (parts per 10 million).
rm – average number of rooms per dwelling.
age – proportion of owner-occupied units built prior to 1940.
dis – weighted distances to five Boston employment centres.
rad – index of accessibility to radial highways.
tax – full-value property-tax rate per $10,000.
ptratio – pupil-teacher ratio by town.
black – 1000(Bk – 0.63)^2 where Bk is the proportion of blacks by town.
lstat – % lower status of the population.
medv – Median value of owner-occupied homes in $1000’s.

We would like to predict medv using the information contained in the other 13 variables. As a result, one of the first tasks is to examine the correlations between medv and each of the other 13 variables. We also fit a linear model by regressing medv on the 13 predictors. Table 1 in the next section shows the results.

The Calculation

Table 1 – How medv and the Predictors Relate

Predictor….	Coefficient Estimate.	t-Statistic.	Correlation with medv	Partial Correlation with medv
Intercept	36.36	7.1440
crim	-0.1080	-3.287	-0.3883	-0.1467
zn	0.04642	3.382	0.3604	0.1509
indus	0.02056	0.3340	-0.4837	0.01510
chas	2.6870	3.1180	0.1753	0.1393
nox	-17.7700	-4.6510	-0.4273	-0.2054
rm	3.8100	9.1160	0.6954	0.3805
age	0.0006922	0.0520	-0.3770	0.002347
dis	-1.4760	-7.3980	0.2499	-0.3167
rad	0.0306	4.6130	-0.3816	0.2038
tax	-0.01233	-3.2800	-0.4685	-0.1464
ptratio	-0.9527	-7.2830	-0.5078	-0.3122
black	0.009312	3.4670	0.3335	0.1546
lstat	-0.5248	-10.3470	-0.7377	-0.4231

The above table is based on a least squares linear model regressing the response medv on all 13 predictors and the correlation matrix (found below in the section for R code). The first column shows the names of the predictors. The next two columns show the coefficient estimates (orange header) and the t-statistics (blue header), respectively. The next column (green header) shows the correlations between the response medv and the predictors. The last column (lime header) shows the partial correlations between medv and the predictors, which is the focus of this discussion.

The coefficient estimate of any given predictor can be regarded as the impact of that predictor on the response holding all the other predictors fixed. Since all the predictors are numeric, the coefficient estimate of a predictor can be interpreted as expected change in the response medv (house price in $1000s) per unit increase in that predictor, given all other predictors are held fixed. Each t-statistic in the above table is way to test the significance of that predictor after controlling for the other predictors. Larger values of the t-statistic indicate significance of the predictor in question. The larger the t-statistic, the stronger the evidence for rejecting the null hypothesis that the coefficient of that predictor is zero. Because of the large sample size (506), a t-statistic (in absolute value) of around 2 corresponds to roughly a 5% p-value. Thus, all the t-statistics except for two are considered large.

The correlation is a measure of the relationship between each predictor and the response medv. This is where things get interesting. These are examples of associations that can disappear or reverse when other predictors are taken into consideration. Four of the correlations are highlighted in red (indus, age, dis and rad). The following shows the scatterplots for these 4 pairs of variables.

Selected Scatterplots – Boston Dataset

Take dis as an example with a correlation of 0.2499 with the response (the bottom left scatterplot). This indicates a positive association between dis and medv, the higher weighted distance to the employment centers, the higher the house price. Yet the coefficient estimate of dis is -1.4760, indicating a negative impact of dis on house prices. For each unit increase in weighted distance to the employment centers, the median house price is expected to decrease by $1,470, holding all other predictors fixed. Thus, the regression model is giving an opposite message, i.e., the higher the weighted distance, the lower the house prices. This seems like a contradiction.

It turns out that there is no contradiction. The correlation of 0.2499 between dis and medv is computed with only these two variables without filtering out the effects of other variables that might impact house prices. The positive correlation of 0.2499 may not be an accurate measure of the relationship between medv and dis because the positive correlation could be due to the influence of variables other than dis (these other influential variables are called confounding variables). The coefficient estimate of -1.4760 takes into the confounding variables into account. Once the these confounding variables are accounted for, the association between medv and dis reverses (in this case from positive to negative).

The other reversal of association indicated in Table 1 above is the correlation between medv and rad (-0,3816). The two examples of the disappearance of association are that between medv and indus (-0.4837) and the association between medv and age (-0.3770). The disappearance or reversal is due to the removal of the effects of the other predictors. To gain better insight, let’s focus on the column for partial correlation in the above table.

Partial Correlation

The partial correlation between the response and a predictor is the correlation between the two variables after the effects of the other predictors have been removed. Suppose we work with the following linear regression model.

(1)

……….

$Y=\beta_0+\beta_1 X_1+\beta_2 X_2 + \cdots+\beta_p X_p+ \epsilon$

The partial correlation between the response $Y$ and the predictor $X_j$ , $j=1,2,\cdots,p$ , is defined to be the correlation between $Y$ and $X_j$ after the effects of all other predictors are taken into account. How do we remove the effects of the predictors that are not $X_j$ ? There are two ways to calculate the partial correlation. One way is to use the t-statistics from regressing the response on all the available predictors.

Consider the linear model described in (1) where the response is $Y$ and the predictors are $X_1,X_2,\cdots,X_p$ . We fit one model by regressing $Y$ on all predictors $X_i$ . Let $t_i$ denote the t-statistic for the coefficient $\beta_i$ associated with the predictor $X_i$ . The following gives the partial correlation between the response $Y$ and the predictor $X_j$ after controlling for the effects of the other predictors.

(2)

……….

$\displaystyle \text{Partial Correlation}=\frac{t_j}{\sqrt{t_j^2+n-p-1}}$

where $n$ is the number of observations in the sample and $j=1,2,\cdots,p$ . The values of the partial correlation in Table 1 above are calculated according to (2) using the t-statistics in that table along with $n$ = 506 and $p$ = 13.

In Table 1 above, we see that the partial correlation between medv and dis is -0.3167, after controlling for the effects of the other predictors. In contrast, before accounting for the effects of the other 12 predictors in the Boston dataset, the correlation between medv and dis is 0.2499, which is due not to a causal relationship but rather to one or more additional predictors that influence both variables to be large. Once these additional predictors are accounted for, the coefficient estimate of dis is -1.476 (same sign as the partial correlation).

We would like to point out that the predictors indus and age, their partial correlations with medv are close to zero, 0.01510 and 0.002347, respectively. After filtering out the effects of the other predictors, these two predictors are virtually uncorrelated with the response medv. The two predictors are the disappearance examples, while dis and rad are the reversal examples.

Partial Correlation – A Better Definition

The partial correlation, as described in (2), is easy to calculate. We only need to fit one linear regression model with all the available predictors and make use of the t-statistics. However, the definition (2) does not provide much insight about the partial correlation. For example, how does the partial correlation reflect the removal of the effects of the confounding variables? We now present an equivalent definition that provides much more insight. This definition also provides a scatterplot called added variable plot, which provides a graphical representation of the relation between the response and the predictor. To calculate one value of partial correlation, the alternative definition requires fitting two linear regression models, as described below.

Consider the linear relationship described in (1).
Regress the response $Y$ on $X_1,\cdots,X_{j-1},X_{j+1},\cdots,X_p$ , i.e., regress $Y$ on all predictors $X_i$ where $i \ne j$ . The residuals from this regression are denoted by $d_1,\cdots,d_n$ where $d_i=y_i-\hat{y_i}$ .
Regress the predictor $X_j$ on $X_1,\cdots,X_{j-1},X_{j+1},\cdots,X_p$ , i.e., regress $X_j$ on all predictors $X_i$ where $i \ne j$ . The residuals from this regression are denoted by $e_1,\cdots,e_n$ where $e_i=x_{ij}-\hat{x_{ij}}$ , i.e., the $i$ th observed value of $X_j$ minus the $i$ th fitted value of $X_j$ .
The partial correlation between $Y$ and $X_j$ after controlling the effects of all other predictors is the sample correlation between the residuals $d_i$ and the residuals $e_i$ . The precise calculation of the partial correlation is given below.

(3)

……….

$\displaystyle \begin{aligned} \text{Partial Correlation}&=\text{Corr}(d,e) \\&=\frac{\sum_{i=1}^n (d_i-\overline{d}) (e_i-\overline{e})}{\sqrt{\sum_{i=1}^n (d_i-\overline{d})^2 \sum_{i=1}^n (e_i-\overline{e})^2}} \\&=\frac{\sum_{i=1}^n d_i \ e_i}{\sqrt{\sum_{i=1}^n d_i^2 \sum_{i=1}^n e_i^2}} \end{aligned}$

The plot of the residuals $d_i$ against the residuals $e_i$ is called the added variable plot, also called partial regression plot. The added variable plot provides a graphical representation of the relation between the response $Y$ and the predictor $X_j$ .

The definition using (3) is a lot more computationally involved. Each partial correlation calculation will require the fitting of two separate regression models. For the Boston dataset example, computing the partial correlations for the 13 predictors will require fitting 26 linear regression models. With the definition (2), we just need to fit one regression model housing all the predictors.

Though (2) is a simpler and more efficient way of calculating partial correlation, (2) by itself does not provide any intuition of partial correlation. Consider the two sets of residuals in (3). The first set of residuals results from regressing the response $Y$ on the other predictors and represents the part of $Y$ that is not explained by the other predictors. The second set of residuals results from regressing the predictor $X_j$ on the other predictors and represents the part of $X_j$ that is not explained by the other predictors. Thus the two sets of residuals contain no information about the predictors $X_i$ where $i \ne j$ . Thus, the information about all the other predictors are removed from these 2 sets of residuals. Therefore, the correlation between these 2 sets of residuals will give a more accurate picture of the relationship between the response and a given predictor.

Another advantage of the (3) as a definition of partial correlation is that the added variable plot will provide a visual confirmation of the partial correlation. The added variable plot not only provides a visualization of the relationship between the two sets of residuals. It also gives us a way to check for outliers and other data anomalies among the residuals. The following shows the added residual plots from the Boston dataset along with the scatterplots before adjusting for other predictors. Four added variable plots are shown below.

Selected Added Variable Plots – Boston Dataset

In each of the added variable plots, the trend line reflects the partial correlation. For example, in the upper two plots, the trend lines are horizontal, corresponding to the partial correlations of 0.01510 (between medv and indus) and 0.002347 (between medv and age). In the bottom two plots, the trend lines are clearly not horizontal, corresponding to the partial correlations of -0.3167 (between medv and dis) and 0.2038 (between medv and rad). Four side-by-side comparisons are shown below – the scatterplot prior to the removal of the effects of the other predictors and the corresponding added variable plot.

Scatterplot and Added Variable Plot

The above comparison shows a disappearance of the association between indus and medv. The trendline in the scatterplot on the left refers to a correlation of -0.4837. The added variable plot on the right refers to a partial correlation of 0.01510, which is virtually zero. After the effects of the other predictors are filtered out, indus and medv are no longer related.

Scatterplot and Added Variable Plot

The above comparison shows a disappearance of the association between age and medv. The trendline in the scatterplot on the left refers to a correlation of -0.3770. The added variable plot on the right refers to a partial correlation of 0.002347, which is virtually zero. After the effects of the other predictors are filtered out, age and medv are no longer related.

Scatterplot and Added Variable Plot

The above comparison shows a reversal of the association between dis and medv. The trendline in the scatterplot on the left refers to a correlation of 0.2499. The added variable plot on the right refers to a partial correlation of -0.3167. After the effects of the other predictors are filtered out, dis and medv are actually negatively associated.

Scatterplot and Added Variable Plot

The above comparison shows a reversal of the association between rad and medv. The trendline in the scatterplot on the left refers to a correlation of -0.3816. The added variable plot on the right refers to a partial correlation of 0.2038. After the effects of the other predictors are filtered out, rad and medv are actually positively associated.

Remarks

With the insight from the preceding section, we contrast the last two columns in Table 1. The predictors age and indus are negatively correlated with medv. Once all the other predictors that are negatively associated with median house price have been filtered out, age and indus are no longer correlated with median house price as evidenced by their partial correlations, 0.01510 and 0.002347, which are virtually zero. The small partial correlations also indicate that the t-statistics are small according to (2). Thus, age and indus are not related to median house prices after the effects of the other predictors are filtered out.

The definition (2) also provides valuable information since it is based on the t-statistics. A large value of the partial correlation indicates that the response $Y$ is strongly associated with the predictor $X_j$ and the association is “real” since the other predictors are no longer “muddying” the water. On the other hand, the t-statistic $t_j$ in (2) is for testing the significance of the coefficient $\beta_j$ . A large value of $t_j$ indicates that the coefficient $\beta_j$ is significant after adjusting for the effects of the other predictors. According to (2), a large value of the t-statistic corresponds to a large value of the partial correlation. Thus, the partial correlation and the t-statistic convey essentially the same information.

From the preceding paragraph, we know that small partial correlation means that the corresponding t-statistic is small and as a result, the variable in question is an insignificant predictor of the response. The variables age and indus have very small partial correlations with the response. Hence they are not statistically significant in the large model housing all the predictors. However, age and indus are significant if we regress medv on each of these two variables by itself. In fact, regressing medv on indus alone tells us that the coefficient estimate is -0.64849 and the t-statistic is -12.41, which is a very large t-statistic. Furthermore, regressing medv on age alone tells us that the coefficient estimate is -0.12316 and the t-statistic is -9.137. Here we have an interesting contrast between simple linear regression (SLR) and multiple linear regression (MLR). It is possible that a predictor is statistically significant in SLR but is insignificant in MLR. This discrepancy is intimately related to the notion of partial correlation. This is discussed in the next post. The summary of four SLR models are shown below.

Summary of SLRs – Boston Dataset

Another observation we would like to make is that even if there is no reversal in association, the partial correlation is usually smaller than the correlation without accounting for the effects of the other predictors. For example, the predictor lstat is the most significant predictor of medv with a t-statistic of -10.3470. The correlation with the response medv is -0.7377, indicating that the higher the percent of lower status of the population, the lower the median house price. This makes intuitive sense. This negative association still holds after filtering out the effects of the other predictors, resulting in a partial correlation of -0.4231, which is a modest level of correlation. Note that of all the predictors, lstat has the strongest correlation with medv both before and after adjusting for the other predictors.

The above discussion shows that there is a statistical phenomenon such that the association between two variables can emerge, disappear or reverse when the effects of other predictors are accounted for. With the Boston dataset, we give two examples of association that disappear and two examples of association that reverse. The disappearing examples are the association between indus and medv and the association between age and medv. In both of these examples, the partial correlations are near 0 (thus the disappearance of the association). The two examples of reversals are the association between dis and medv and the association between rad and medv. Such examples of disappearance or reversals seem like a contradiction. It turns out that there is no contradiction. Once we understand the roles played by the confounding variables, it is a matter of interpreting the correlations properly.

The Section for R Code

This section provides the R code for the graphical displays and the calculation used in the discussion above. See here for information on the fitting of multiple regression model.

First, we read the Boston dataset into a data frame and produce a summary of the dataset. This is performed by the following R code (Code Chunk 1).

library(MASS)
data(“Boston”)
summary(Boston)

The result of the summary of the Boston dataset is not shown. Since the topic is correlation, we would like to examine the correlation matrix. The following R code produces the correlation matrix along with scatterplots of medv against several interesting variables (Code Chunk 2).

cor(Boston)
……….
library(ggplot2)
library(gridExtra)
p1 <- ggplot(Boston, aes(x=indus,y=medv)) + geom_point() +
….. geom_smooth(method = “lm”, se = FALSE)
p2 <- ggplot(Boston, aes(x=age,y=medv)) + geom_point() +
….. geom_smooth(method = “lm”, se = FALSE)
p3 <- ggplot(Boston, aes(x=dis,y=medv)) + geom_point() +
….. geom_smooth(method = “lm”, se = FALSE)
p4 <- ggplot(Boston, aes(x=rad,y=medv)) + geom_point() +
….. geom_smooth(method = “lm”, se = FALSE)
……….
grid.arrange(p1, p2, p3, p4, ncol = 2)

….
Correlation Matrix – Boston Dataset

Selected Scatterplots – Boston Dataset

The calculation of the partial correlation involves regressing the response medv on 13 predictors. This is performed by the following R code (Code Chunk 3).

MLR.Boston <- lm(medv ~ ., Boston)
summary(MLR.Boston)

Summary of MLR – All Predictors in Boston Dataset

Table 1 above is based on information obtained from running Code Chunk 2 and Code Chunk 3. The partial correlations found in Table 1 can be calculated using formula (2) above based on the t-statistics. Though (2) is a computationally efficient way to calculate partial correlation, the calculation (3) using the residuals from two fitted models provide more insight. The following code chunks produce the added variable plots.

# Code Chunk 4
# Partial correlation for medv and indus
MLR1 <- lm(medv ~ .-indus, data = Boston)
MLR2 <- lm(indus ~ .-medv, data = Boston)
cor(MLR1$residuals, MLR2$residuals)
……….
Resid.Y <- MLR1$residuals
Resid.indus <- MLR2$residuals
ggplot(Boston, aes(x=Resid.Y,y=Resid.indus)) + geom_point() +
….geom_smooth(method = “lm”, se = FALSE)

# Code Chunk 5
# Partial correlation for medv and age
MLR1 <- lm(medv ~ .-indus, data = Boston)
MLR2 <- lm(age ~ .-medv, data = Boston)
cor(MLR1$residuals, MLR2$residuals)
……….
Resid.Y <- MLR1$residuals
Resid.age <- MLR2$residuals
ggplot(Boston, aes(x=Resid.Y,y=Resid.age)) + geom_point() +
….geom_smooth(method = “lm”, se = FALSE)

# Code Chunk 6
# Partial correlation for medv and dis
MLR1 <- lm(medv ~ .-dis, data = Boston)
MLR2 <- lm(dis ~ .-medv, data = Boston)
cor(MLR1$residuals, MLR2$residuals)
……….
Resid.Y <- MLR1$residuals
Resid.dis <- MLR2$residuals
ggplot(Boston, aes(x=Resid.Y,y=Resid.dis)) + geom_point() +
….geom_smooth(method = “lm”, se = FALSE)

# Code Chunk 7
# Partial correlation for medv and rad
MLR1 <- lm(medv ~ .-rad, data = Boston)
MLR2 <- lm(rad ~ .-medv, data = Boston)
cor(MLR1$residuals, MLR2$residuals)
……….
Resid.Y <- MLR1$residuals
Resid.rad <- MLR2$residuals
ggplot(Boston, aes(x=Resid.Y,y=Resid.rad)) + geom_point() +
….geom_smooth(method = “lm”, se = FALSE)

Selected Added Variable Plots – Boston Dataset

Some of the predictors are statistical significant in the SLR setting but not in MLR. The following R code produces 4 SLR models (Code Chunk 8).

# Code Chunk 8
SLR.indus <- lm(medv ~ indus, Boston)
summary(SLR.indus)
……….
SLR.age <- lm(medv ~ age, Boston)
summary(SLR.age)
……….
SLR.dis <- lm(medv ~ dis, Boston)
summary(SLR.dis)
……….
SLR.rad <- lm(medv ~ rad, Boston)
summary(SLR.rad)

Summary of SLRs – Boston Dataset

……….

This article was originally published here in a companion site.

$\text{ }$

Dan Ma added variable plot/font>
Daniel Ma added variable plot

Dan Ma partial correlation
Daniel Ma partial correlation

Dan Ma linear regression
Daniel Ma linear regression

Dan Ma simple linear regression
Daniel Ma simple linear regression

Dan Ma multiple linear regression
Daniel Ma multiple linear regression

$\copyright$ 2024 – Dan Ma

Posted: October 11, 2024

Multiple linear regression – adjusted R-squared

Posted on October 10, 2024 by Dan Ma

The $R^2$ statistic (discussed here), the coefficient of determination, is the ratio of the regression sum of squares (Reg SS) over the total sum of squares (TSS). Reg SS measures the amount of variability in the response $Y$ that is explained by performing the linear regression. As a result, the $\text{R}^2$ statistic is the measures of the proportion of the total variability in the response that is explained by performing the linear regression. The $R^2$ is between 0 and 1 and is a measure of the quality of the linear regression model. An $R^2$ that is near 1 indicates that a large proportion of the variability in the response has been explained by the regression model. An $R^2$ nears 0 indicates that linear regression did not explain much of the variability in the response. The higher the $R^2$ , the better the model quality. However, the $R^2$ statistic suffers from a serious shortcoming, which is that adding new predictors to the linear regression model always increases $R^2$ (at least not decreasing it). Thus, using $R^2$ as a criterion for selecting the best model would lead us to conclude that the best linear model is the one that houses all available predictors even if some of the added predictors are just noise variables, i.e., variables that have very little predictive power. We define adjust $R^2$ and show how it is a better measure of model quality. With the adjusted $\text{R}^2$ , a penalty is imposed on predictors that only add a minuscule amount of predictive power. As a result, the adjusted $\text{R}^2$ is a measure of model quality that is more capable of distinguishing predictors that have sufficient predictive power from those that do not.

Adjusted R-Squared

Recall that $R^2$ is defined as follows:

(1)

……….

$\displaystyle \text{R}^2=1-\frac{\text{RSS}}{\text{TSS}}$

……….
where $\text{TSS}=\sum_{i=1}^n (y_i-\overline{y})^2$ and $\text{RSS}=\sum_{i=1}^n (y_i-\hat{y_i})^2$ . TSS measures the total variability in the response $Y$ , i.e., the amount of variability inherent in the response before linear regression is performed. On the other hand, RSS measures the amount of variability in the response that is left unexplained after performing the regression. Regression sum of squares (RSS) is the difference between TSS and RSS (TSS – RSS), which is a measure of the variability in the response that is explained by the linear regression. Thus, by definition the $R^2$ statistic measures the proportion of the variability in the response that is explained by performing the linear regression.

As noted above, the $\text{R}^2$ statistic suffers from one serious shortcoming. Note that RSS always decreases when more predictors are added to the model. As a result, the $R^2$ statistic will always improve (increase) as more predictors are added to the model even if the added predictors do not add much predictive power. The model with all the available predictors will always have the lowest RSS and hence the largest $\text{R}^2$ . Thus, RSS and $\text{R}^2$ are not suitable as a model selection criterion. Among a collection of models with different number of predictors, the best model according to $\text{R}^2$ will be the one with all of the predictors. The adjusted $R^2$ is to remedy this shortcoming. It is defined as follows:

(2)

……….

$\displaystyle \text{Adjusted R}^2=1-\frac{\text{RSS} / (n-p-1)}{\text{TSS} / (n-1)}$

……….
The adjusted $R^2$ is derived by replacing the sums of squares in the definition of $R^2$ with the corresponding mean squares in the ANOVA table, i.e., replacing RSS by the error mean square and replacing TSS by the total mean square, which is the sample variance. Thus we can express adjusted $R^2$ as follows:

(3)

……….

$\displaystyle \text{Adjusted R}^2=1-\frac{\text{RSS} / (n-p-1)}{\text{TSS} / (n-1)}=1-\frac{s^2}{s_y^2}$

……….

where $s^2$ is the mean squared error (MSE) from the ANOVA table and $s_y^2$ is the sample variance of the response $Y$ . The adjusted $R^2$ can also be expressed using the ordinary $R^2$ as follows:

(4)

……….

$\displaystyle \text{Adjusted R}^2=1-\frac{n-1}{n-p-1} \ (1-\text{R}^2)$

……….

Penalty for Additional Predictors

The $\text{R}^2$ statistic always increases as more predictors are added to the model. The adjust $\text{R}^2$ does not necessarily increases as more predictors are added. The adjusted $\text{R}^2$ works by imposing a penalty on any additional predictor. Consider the dynamics in the fraction $\frac{\text{RSS}}{n-p-1}$ . The denominator decreases when predictors are added to the model (the $p$ increases). When the drop in RSS does not compensate for the decrease in the denominator, the ratio $\frac{\text{RSS}}{n-p-1}$ will increase and the adjust $\text{R}^2$ will decrease. Take the extreme example of adding a completely irrelevant predictor to the model. The RSS will not change but the denominator in $\frac{\text{RSS}}{n-p-1}$ is decreased by 1. Hence the adjusted $\text{R}^2$ will decrease, rather than increase. For the addition of a new predictor to increase the adjusted $\text{R}^2$ , the decrease in the RSS of the model must be large enough to offset the increase in $p$ . Thus, the adjusted $\text{R}^2$ is more capable of differentiating the predictors with predictive power from those that do not.

Similar to $\text{R}^2$ , a larger adjusted $\text{R}^2$ is preferred over a model with smaller adjusted $\text{R}^2$ . Given a collection of linear models fitted to the same dataset and with varying numbers of predictors, if using the adjusted $\text{R}^2$ as a model selection criterion, the model with the largest adjusted $\text{R}^2$ is the best model. In theory, we assume that the model with the largest adjusted $\text{R}^2$ will have only the correct predictors and will have no noise variables. This is the result of the adjusted $\text{R}^2$ imposing a penalty for the inclusion of unnecessary predictors into the model.

As predictors with no or little predictive power are added to the model, the resulting model tends to overfit the training data. In the parlance of statistical learning, the adjust $\text{R}^2$ and statistics such as AIC and BIC are to estimate the test error by making an adjustment to the training error (through altering the number of predictors) to account for the bias due to overfitting.

In (3), $s^2$ is the only quantity that depends on the model. Thus, choosing the model with the highest adjusted $\text{R}^2$ is equivalent to choosing the model with the smallest MSE. This means that we can achieve the same result by using MSE as a measure of model quality (the smaller the MSE, the better the model quality).

….

This article was published here in a companion site.

Dan Ma adjusted R-squared
Daniel Ma adjusted R-squared

Dan Ma linear regression
Daniel Ma linear regression

Dan Ma simple linear regression
Daniel Ma simple linear regression

Dan Ma multiple linear regression
Daniel Ma multiple linear regression

$\copyright$ 2024 – Dan Ma

Posted: October 10, 2024

Multiple linear regression – model fitting and basic topics

Posted on October 9, 2024 by Dan Ma

Simple linear regression (SLR) is a useful approach for predicting a response $Y$ on the basis of a single predictor variable $X$ . We often work with datasets with more than one predictor. In such datasets, the response variable is potentially related to multiple predictor variables. If the multiple factors influence the behavior of a variable of interest, we certainly would like to incorporate these factors into the analysis and potentially obtain better predictions. Multiple linear regression (MLR) is a natural extension of simple linear regression in that MLR can directly accommodate the multiple predictors. We discuss the MLR model fitting process and other basic topics such as inference, measures of model fit and predictions.

The discussion on MLR mirrors that on simple linear regression presented earlier. Here’s some contents on SLR: model fitting, inference, measures of model fit and predictions.

MLR Model Equation

With MLR, we would like to use the multiple predictors that are available at our disposal to obtain a better understanding of the behavior of the response variable. As the name suggests, the regression function in MLR is a linear function of the predictors. In MLR, the following model equation relates the response $Y$ to the predictors $X_1,X_2,\cdots,X_p$ .

(1)

……….

$Y=\beta_0+\beta_1 X_1+\beta_2 X_2 + \cdots+\beta_p X_p+ \epsilon$

where $p$ is the number of predictors, $\beta_0$ is the intercept coefficient, $\beta_j$ is the model coefficient for the $j$ th predictor $X_j$ for $j=1,2,\cdots,p$ and $\epsilon$ is the random error term. The random error term $\epsilon$ is normally distributed with zero mean and constant variance $\sigma^2$ .

In a sample of $n$ observations, the $i$ th observation is of the form $(Y_i,X_{i1},X_{i2},\cdots,X_{ip})$ . As a result, there are $n$ model equations, one for each observation, that can be expressed as follows.

(2)

……….

$Y_i=\beta_0+\beta_1 X_{i1}+\beta_2 Y_{i2} + \cdots+\beta_p X_{ip}+ \epsilon_i$

(2)……….where $i=1,2,3,\cdots,n$

Note that the random errors $\epsilon_i$ , in addition to being normally distributed with mean 0 and variance $\sigma^2$ , are independent. The above $n$ model equations can be better described in matrix notation.

(3)

……….

(.)

……….

$\bold Y= \bold X \bold \beta + \bold \epsilon$

where

$\bold Y$ is the $n \times 1$ vector of response variables,
$\bold X$ is the matrix containing the values of the predictor variables, an $n \times (p+1)$ matrix,
$\beta$ is the $(p+1) \times 1$ vector of model coefficients,
$\epsilon$ is the $n \times 1$ vector of random error terms.

Note that the matrix $\bold X$ is also known as the data matrix or design matrix.

Estimating MLR Coefficients using Least Squares Criterion

To estimate the unknown coefficients indicated in (1), we choose $\beta_0,\beta_1,\cdots,\beta_n$ such that the following sum of the squared differences is minimized.

(4)

……….

$\displaystyle \sum_{i=1}^n \biggl[Y_i-(\beta_0+\beta_1 X_{i1}+\beta_2 X_{i2}+ \cdots + \beta_p X_{ip}) \biggr]^2$

To minimize (4), take partial derivative of it with respect to each coefficient $\beta_j$ , resulting in $p+1$ expressions in $p+1$ unknowns. Setting these expressions equal to zero results in a system of $p+1$ simultaneous equations in $p+1$ unknowns. These equations are known as the normal equations. The normal equations can be compactly expressed in matrix form as follows.

(5)

……….

$\displaystyle ( \bold X^{'} \bold X) \hat{\beta}=\bold X^{'} \bold Y$

(6)

……….

$\displaystyle \hat{\beta}=( \bold X^{'} \bold X)^{-1} \bold X^{'} \bold Y$

When (6) is carried out for $p=1$ (simple linear regression), we obtain the following.

$\bold X^{'} \bold X=\left[\begin{array}{ccccc} 1 & 1 & 1 & \cdots & 1\\ X_1 & X_2 & X_3 & \cdots & X_n \end{array}\right] \left[\begin{array}{cc} 1 & X_1 \\ 1 & X_2 \\ \vdots & \vdots \\ 1 & X_n \end{array}\right]= \left[\begin{array}{cc} n & n \overline{X} \\ n \overline{X} & \sum_{i=1}^n X_i^2 \end{array}\right]$

……….

$(\bold X^{'} \bold X)^{-1}=\frac{1}{\displaystyle n \sum_{i=1}^n X_i^2- \biggl(n \overline{X} \biggr)^2} \left[\begin{array}{cc} \sum_{i=1}^n X_i^2 & -n \overline{X} \\ -n \overline{X} & n \end{array}\right]=\frac{1}{\displaystyle n \sum_{i=1}^n (X_i-\overline{X})^2} \left[\begin{array}{cc} \sum_{i=1}^n X_i^2 & -n \overline{X} \\ -n \overline{X} & n \end{array}\right]$

……….

$\bold X^{'} Y=\left[\begin{array}{c} n \overline{Y} \\ \sum_{i=1}^n X_i Y_i \end{array}\right]$

……….

$\displaystyle \begin{aligned} \hat{\beta}&=( \bold X^{'} \bold X)^{-1} \bold X^{'} \bold Y \\&=\frac{1}{\displaystyle n \sum_{i=1}^n (X_i-\overline{X})^2} \left[\begin{array}{cc} \sum_{i=1}^n X_i^2 & -n \overline{X} \\ -n \overline{X} & n \end{array}\right] \left[\begin{array}{c} n \overline{Y} \\ \sum_{i=1}^n X_i Y_i \end{array}\right] \\&=\left[\begin{array}{c} \overline{Y}-\hat{\beta_1} \overline{X} \\ \hat{\beta_1} \end{array}\right] \end{aligned}$

where $\displaystyle \hat{\beta_1}=\frac{\displaystyle \sum_{i=1}^n X_i Y_i-n \overline{X} \ \overline{Y}}{\displaystyle \sum_{i=1}^n X_i^2 - n \overline{X}}$ . Such simple forms for regression coefficient estimates are only possible for simple linear regression. When there are more than one predictor ( $p>1$ ), the forms for the coefficient estimates are complicated and cannot be easily stated as those for SLR. Typically, we use software solutions such R to obtain the least squares estimates.

Analysis of Variance

In simple linear regression, the total sum of squares (TSS) measures the total variance in the response $Y$ . It can be thought of as the amount of variability inherent in the response before the regression is performed. TSS can be partitioned into two parts. One of which is the residual sum of squares (RSS), which is the amount of variability that is left unexplained after performing the regression. The other part is the regression sum of squares (Reg SS), which is the amount of variability in the response $Y$ that is explained by performing the linear regression. The following relation holds true for simple linear regression as well as multiple linear regression (see here for the derivation).

(7)

……….

$\sum_{i=1}^n (y_i-\overline{y})^2=\sum_{i=1}^n (y_i-\hat{y_i})^2 + \sum_{i=1}^n (\hat{y_i}-\overline{y})^2$

(7)

……………….

TSS

………………….

RSS

………………..

Reg SS

The relation (7) can also be expressed in matrix form.

(8)

……….

$\bold Y^{'} \bold Y-n \overline{Y}^2=(\bold Y-\hat{\bold Y})^{'} (\bold Y-\hat{\bold Y}) + (\hat{\bold Y}^{'} \hat{\bold Y}-n \overline{Y}^2)$

(8)

……………..

TSS

…………………..

RSS

…………………..

Reg SS

As a result, the ANOVA table for MLR is similar in structure to the ANOVA table for SLR, except for the degrees of freedom due to the fact that there are more predictors in MLR.

ANOVA Table

Source of Variation…..	Sum of Squares…..	df………………….	Mean Square…..	F-statistic
$\text{Regression}$	$\text{Reg SS}$	$p$	$\text{Reg SS} /p$	$\displaystyle \frac{\text{Reg SS}/p}{s^2}$
$\text{Error}$	$\text{RSS}$	$n-p-1$	$s^2=\text{RSS}/(n-p-1)$
$\text{Total}$	$\text{TSS}$	$n-1$

The quantity $s^2$ is the mean square error (MSE), which is an unbiased estimator of the variance $\sigma^2$ of the random error term.

F-Test

The F-statistic in the ANOVA table can be used for testing the following hypotheses.

(9a)

……….

$H_0:\beta_1=\beta_2= \cdots =\beta_p=0$

…..

(Null Hypothesis)

versus

(9b)

……….

$H_1: \beta_j \ne 0$

$j$

………

(Alternative Hypothesis)

Under the null hypothesis $H_0$ , the F-statistic follows an F distribution with degrees of freedom $p$ and $n-p-1$ with expected value close to 1. For any given value of $n$ and $p$ , we can use statistical software package to compute the p-value associated with the F-statistic. Based on the computed p-value, we can determine whether or not to reject the null hypothesis $H_0$ .

When there is no relationship between the response and predictors, one would expect the F-statistic to take on a value close to 1. On the other hand, if $H_0$ is not true, we expect the F-statistic to be substantially greater than 1. Thus, an F-statistic that is sufficiently larger than 1 (as evidenced by a small p-value) provides evidence for rejecting the null hypothesis and in favor of the alternative hypothesis, which is that the linear regression model is useful for explaining the variation in the response $Y$ . Rejecting the null hypothesis $H_0$ means that we have strong evidence that at least one of the $p$ predictors is an important predictor for the response variable $Y$ . However, the F-test does not indicate which of the predictors is important.

A test of the null hypothesis $H_0$ is a test of the “significance of the regression”. This is because if $H_0$ is not rejected, we must conclude that the predictor variables as a group do not improve the prediction of the response $Y$ over that obtained by regressing $Y$ on the intercept coefficient $\beta_0$ (intercept-only model). Note that the least squares estimate would be $\hat{\beta_0}=\overline{Y}$ when using the intercept-only model. Thus, if the null hypothesis (9a) is shown to be not significant, the best estimate of the response would be its sample mean.

How large does the F-statistic have to be in order to reject $H_0$ and make the conclusion that there is a relationship? It turns out that the answer depends on the values of $n$ and $p$ . When $n$ is large, an F-statistic that is just a little larger than 1 might still provide evidence against $H_0$ . In contrast, a larger F-statistic is needed to reject $H_0$ if the sample size $n$ is small.

The Coefficient of Determination

The $R^2$ statistic, the coefficient of determination, is defined similarly as for simple linear regression. It is the ratio of the regression sum of squares (Reg SS) over the total sum of squares (TSS). Thus, by definition, $R^2$ measures the proportion of the variation in the response variable about its mean that can be explained by the least squares regression model.

(10)

……….

$\displaystyle \text{R}^2=\frac{\text{Reg SS}}{\text{TSS}}=1-\frac{\text{RSS}}{\text{TSS}}$

In simple linear regression, the $\text{R}^2$ statistic is the square of the sample correlation (between the response $y$ and single predictor $x$ ). In MLR, this is no longer the case since there are more than one predictor. However, the $\text{R}^2$ statistic is the square of the sample correlation between the observed $y$ and the fitted values $\hat{y}$ . That is $\text{R}^2=\text{Corr}(y,\hat{y})^2$ .

Sampling Distribution of the Least Squares Estimator

The least squares estimator $\hat{\beta}$ follows a multinormal distribution. We make use of the multinormal distribution to draw inference of the vector of least squares estimator $\hat{\beta}$ . From (6), we know that $\hat{\beta}=\bold A \bold Y$ where $\bold A=(\bold X^{'} \bold X)^{-1} \bold X^{'}$ . First, we examine the mean and variance of the least squares estimators $\hat{\beta}$ . The following derivation shows that $\hat{\beta}$ is an unbiased estimator of the unknown $\beta$ .

$\displaystyle \begin{aligned}E(\hat{\beta})&=E(\bold A \bold Y) \\&=\bold A E(\bold Y) \\&=\bold A \bold X \beta \\&=[(\bold X^{'} \bold X)^{-1} \bold X^{'}] \bold X \beta \\&=[(\bold X^{'} \bold X)^{-1}] [\bold X^{'} \bold X] \beta=\beta \end{aligned}$

Thus, the mean vector $E(\hat{\beta})$ is the vector $\beta$ of the unknown coefficients. With respect to the individual predictors, we have $E(\hat{\beta_j})=\beta_j$ for each $j=1,2,\cdots,p$ .

Note that the random errors $\epsilon_i$ described in (2) have means zero and constant variance $\sigma^2$ . They are also independent. As a result, the response $Y_i$ are independent with mean $E(Y_i \lvert X_1,\cdots,X_n)=\beta_0+\beta_1 X_1+\cdots+\beta_n X_n$ (a linear function of the predictors) and constant variance $\sigma^2$ . Thus the variance-covariance matrix of the vector $\bold Y=(Y_1,\cdots,Y_n)^{'}$ is $\text{Var}(\bold Y)=\sigma^2 \ \bold I_n$ where $\bold I_n$ is the $n \times n$ identity matrix. In this matrix, the diagonal elements are $\sigma^2$ and the off-diagonal elements are zero due to the fact that the variable $Y_i$ are independent. With this in mind, we now derive the variance of $\hat{\beta}$ .

$\displaystyle \begin{aligned} \text{Var}(\hat{\beta})&=\text{Var}(\bold A \bold Y) \\&=\bold A \text{Var}(\bold Y) \bold A^{'} \\&=\bold A \sigma^2 \ \bold I_n \bold A^{'} \\&=[(\bold X^{'} \bold X)^{-1} \bold X^{'}] \sigma^2 \bold I_n [(\bold X^{'} \bold X)^{-1} \bold X^{'}]^{'} * \\&=[(\bold X^{'} \bold X)^{-1} \bold X^{'}] \sigma^2 \bold I_n \bold X (\bold X^{'} \bold X)^{-1} \\&=\sigma^2 [(\bold X^{'} \bold X)^{-1} \bold X^{'} \bold X] (\bold X^{'} \bold X)^{-1} \\&=\sigma^2 (\bold X^{'} \bold X)^{-1}\end{aligned}$

For the step with *, we use these matrix facts: $(\bold M \bold N)^{'}=\bold N^{'} \bold M^{'}$ and $(\bold M^{'})^{-1}=(\bold M^{-1})^{'}$ . We also use that fact that $(\bold X^{'} \bold X)^{'}=\bold X^{'} \bold X$ . The last matrix $\sigma^2 (\bold X^{'} \bold X)^{-1}$ is the variance-covariance matrix of the estimator $\hat{\beta}$ .

(11).. $\displaystyle \text{Var}(\hat{\beta})=\left[\begin{array}{ccccc} \text{Var}(\hat{\beta_0}) & \text{Cov}(\hat{\beta_0},\hat{\beta_1}) & \text{Cov}(\hat{\beta_0},\hat{\beta_2}) & \cdots & \text{Cov}(\hat{\beta_0},\hat{\beta_p})\\ \text{Cov}(\hat{\beta_1},\hat{\beta_0}) & \text{Var}(\hat{\beta_1}) & \text{Cov}(\hat{\beta_1},\hat{\beta_2}) & \cdots & \text{Cov}(\hat{\beta_1},\hat{\beta_p}) \\ \text{Cov}(\hat{\beta_2},\hat{\beta_0}) & \text{Cov}(\hat{\beta_2},\hat{\beta_1}) & \text{Var}(\hat{\beta_2}) & \cdots & \text{Cov}(\hat{\beta_2},\hat{\beta_p}) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \text{Cov}(\hat{\beta_p},\hat{\beta_0}) & \text{Cov}(\hat{\beta_p},\hat{\beta_1}) & \text{Cov}(\hat{\beta_p},\hat{\beta_2}) & \cdots & \text{Var}(\hat{\beta_p}) \end{array}\right]=\sigma^2 (\bold X^{'} \bold X)^{-1}$

In the above variance-covariance matrix, the elements on the diagonal are the variances of the individual estimators $\hat{\beta_j}$ and the off-diagonal elements are covariances of pairs of predictors. Thus, the variance of $\hat{\beta_j}$ is the $(j+1)$ st diagonal element of the matrix $\sigma^2 (\bold X^{'} \bold X)^{-1}$ .

(12)

……….

$\displaystyle \text{Var}(\hat{\beta_j})=\sigma^2 \times (j+1)\text{st diagonal element of } (\bold X^{'} \bold X)^{-1}$

An unbiased estimator of $\sigma^2$ is given by $s^2=\frac{\text{RSS}}{n-p-1}$ where $s$ is called the residual standard error (RSE). Thus, replacing $\sigma^2$ by $s^2$ produces the estimated variance-covariance matrix of $\hat{\beta}$ . The following is the estimated variance of the estimator $\hat{\beta_j}$ for $j=1,2,\cdots,p$ .

(13)

……….

$\displaystyle \text{Est Var}(\hat{\beta_j})=s^2 \times (j+1)\text{st diagonal element of } (\bold X^{'} \bold X)^{-1}$

Because the random errors $\epsilon_i$ are uncorrelated and are normally distributed with zero mean and constant variance $\sigma^2$ , the responses $Y_i$ are uncorrelated, with means as a linear function of the predictors and constant variance $\sigma^2$ . As a result, each estimator $\hat{\beta_j}$ is normally distributed with mean $\beta_j$ and variance given in (13). Furthermore, the vector $\hat{\beta}$ follows a multinormal distribution with mean vector $\beta$ and variance-covariance matrix given in (11). This forms the basis for the inference discussed in the next section.

Inference

We are now ready to draw inference on the estimator $\hat{\beta_j}$ for $j=1,2,\cdots,p$ . As stated in the preceding section, each individual estimator $\hat{\beta_j}$ follows a normal distribution with mean $\beta_j$ and variance given in (12) and estimated variance given in (13). Then the following statistic follows a t-distribution with $n-p-1$ degrees of freedom.

(14)

……….

$\displaystyle t=\frac{\hat{\beta_j}-\beta_j}{\sqrt{s^2 \times (j+1) \text{st diagonal element of } (\bold X^{'} \bold X)^{-1}}}$

(14)……….where $j=1,2,\cdots,p$

The t-statistic forms the basis for constructing confidence intervals and performing test of hypotheses on the unknown coefficient $\beta_j$ .

First, consider test of the null hypothesis $H_0:\beta_j=b$ where $b$ is a value specified by the user performing the test. Then under $H_0$ , the t-statistic in (14) with $\beta_j=b$ follows a t-distribution with $n-p-1$ degrees of freedom. We can test $H_0:\beta_j=b$ against $H_1:\beta_j \ne b$ (two-sided test), against $H_1:\beta_j <b$ (one-sided test that is left sided), or against $H_1:\beta_j >b$ (one-sided test that is right-sided).

The null hypotheses of interest are $H_0:\beta_j=0$ where $j=1,2,\cdots,p$ . Each such test is to examine the significance of the coefficient $\beta_j$ and hence the significance of the term $\beta_j X_j$ . Since each of these hypotheses involves only a single parameter, we must be careful in interpreting such tests. The inference on the null hypothesis $H_0:\beta_j=0$ is conditional on all other coefficients in the model. Each hypothesis provides information about whether each individual predictor is related to the response, after adjusting for the other predictors.

For example, suppose $p=5$ and the regression function is as in (1) with $p=5$ . Suppose that $H_0:\beta_5=0$ is rejected. We then would conclude that the inclusion of $X_5$ variable in the regression function significantly improves the prediction of $Y$ versus the model that uses only the predictors $X_1,X_2,X_3,X_4$ . On the other hand, suppose $H_0:\beta_5=0$ is not rejected. Then we would conclude that including the predictor $X_5$ does not significantly improve the prediction of $Y$ versus the model that utilizes only the predictors $X_1,X_2,X_3,X_4$ . Therefore, conclusion of the significance of one predictor is conditional on all the other predictors in the model. The t-statistic given in (14) can only be used to test the significance of one predictor given all other predictors. It cannot be used to test joint significance or insignificance of two predictors.

For $j=1,2,\cdots,p$ , the following is a 100(1- $\alpha$ )% confidence interval for $\beta_j$ :

(15)

……….

$\displaystyle \hat{\beta_j} \pm t_{n-p-1,\frac{\alpha}{2}} \times \sqrt{s^2 \times (j+1) \text{st diagonal element of } (\bold X^{'} \bold X)^{-1}}$

where $c=t_{n-p-1,\frac{\alpha}{2}}$ is the value in the random variable T following a t-distribution with $n--p-1$ degrees of freedom such that $P(T \le c)=1-\frac{\alpha}{2}$ . In other words, $c=t_{\frac{\alpha}{2},n-p-1}$ is the $(100-\frac{\alpha}{2})$ th percentile of the t-distribution with $n-p-1$ degrees of freedom.

With the F-test and the individual t-test given in (14), we now have two different ways to test the significance of the regression. With the F-test, we can test the null hypothesis (9a) versus the alternative (9b). We can also test (9a) by carrying out $p$ individual t-tests $H_0:\beta_j=0$ for $j=1,2,\cdots,p$ simultaneously. Instead of carrying the F-test, why don’t we just perform the individual t-tests? If one of the t-tests has a small p-value, then we know that at least one of the $p$ predictors is related to the response. There is a problem with this approach especially if the number of predictors $p$ is large. The problem is that when carrying out these individual t-tests, it is almost certain that we will conclude that the regression is significant (i.e., at least one predictor is an important predictor of the response) when in reality the response is not related to the predictors at all.

As an illustration, let’s say $p=100$ . Suppose null hypothesis (9a) is true. This means that none of the 100 predictors is truly associated with the response. Suppose that we perform 100 individual t-tests on the 100 predictors according to (14). At 5% level of significance, about 5% of the 100 p-values will be below 0.05 by chance. In other words, we expect to see approximately five small p-values even though there is no relationship between the response and predictors. Assuming the 100 individual t-tests are independent, the probability of making a correct decision on each of the 100 tests is $0.95^{100}$ , which is practically zero (correct decision here is to not reject the null hypothesis). Thus, it is a near certainty that we will conclude some predictors are significant even though in reality none of them is significant predictor of the response. On the other hand, if we use the F-test testing null hypothesis (9a), there is 95% chance of making the correct decision, i.e., not rejecting the null hypothesis (9a). With the F-test, there is a 95% chance of making the correct decision while with the collection of $p$ many individual t-tests, there is a near 0% chance of making the correct decision.

Predictions

Before we discuss predictions, we comment on how to interpret the coefficients in the MLR model. We focus on predictors that are numeric. When the predictor $X_j$ is a numeric variable, we can interpret the coefficient $\beta_j$ as the expected change in the response $Y$ per unit change in the predictor $X_j$ , holding all other predictors fixed.

Given a vector $\bold X_0=(1,X_{01},X_{02},\cdots,X_{0n})^{'}$ of the predictor values, we would like to predict the corresponding response value, say $Y_0$ . As discussed here, the prediction $\hat{Y}_0$ has two sources of variability. The first source of the variability is the variability due to the random error term $\epsilon_i$ . The second source is due to the error in estimating the regression line at $\bold X_0$ . As a result, the following is the standard error of the prediction $Y_0$ .

(16)

……….

$\displaystyle \text{SE}_{\text{Pred}}=\sqrt{s^2 \ [1 + \bold X_0^{'} (\bold X^{'} \bold X)^{-1} \bold X_0 ]}$

The following is the 100(1- $\alpha$ )% prediction interval for $Y_0$ .

(17)

……….

$\displaystyle \begin{aligned}&\hat{Y_0} \pm t_{n-p-1,\frac{\alpha}{2}} \times \text{SE}_{\text{Pred}} \\&=\hat{Y_0} \pm t_{n-p-1,\frac{\alpha}{2}} \times \sqrt{s^2 [1+ \bold X_0^{'} (\bold X^{'} \bold X)^{-1} \bold X_0 ]} \end{aligned}$

We would like to point out that the prediction procedure should be carried out only when the hypothesis (9a) $H_0:\beta_1=\beta_2=\cdots=\beta_p=0$ is rejected. If the regression is not significant, it does not make sense to make prediction on a future response value. In fact, if the regression is not significant, the best predicted response value is $\overline{Y}$ .

….

This article was publised here in a companion site.

$\text{ }$

$\text{ }$
Dan Ma linear regression
Daniel Ma linear regression

Dan Ma simple linear regression
Daniel Ma simple linear regression

Dan Ma multiple linear regression
Daniel Ma multiple linear regression

$\copyright$ 2024 – Dan Ma

Posted October 9, 2024

Simple linear regression – example – ISLR advertising

Posted on October 8, 2024 by Dan Ma

We illustrate concepts and key calculations in simple linear regression using R. We use the Advertising dataset associated with the popular statistical learning textbook ISLR.

The Advertising Dataset

The Advertising dataset is a simulated dataset for teaching purpose, used in the text ISLR to illustrate linear regression. It has four variables, sales for a given product in 200 different markets along the advertising budgets for three different media: TV, radio and newspaper. The dataset is used in ISLR to illustrate multiple linear regression, more specifically whether there is an association between advertising and sales and if so, how the company selling this product can indirectly increase sales by adjusting the advertising budgets. In this example, we use the Advertising dataset to illustrate simple linear regression. Multiple linear regression on the advertising dataset will be handled separately. All calculations and graphical displays are produced in R.

Information in simple linear regression:

Univariate Data Exploration

The following R code loads the dataset, list the first several rows of the data and produces a numerical summary.

advert <- read.csv(url(“see link below”))
…………………
.head(advert)
…………………
.summary(advert).

link for the CSV file

The First Six Rows

…..	TV………	Radio…….	Newspaper..	Sales
1	230.1	37.8	69.2	22.1
2	44.5	39.3	45.1	10.4
3	17.2	45.9	69.3	9.3
4	151.5	41.3	58.5	18.5
5	180.8	10.8	58.4	12.9
6	8.7	48.9	75.0	7.2

The Summary of the Four Variables

…..	TV………	Radio…….	Newspaper..	Sales
Minimum	0.70	0.000	0.30	1.60
1st Quartile	74.38	9.975	12.75	10.38
Median	149.75	22.900	25.75	12.90
Mean	147.04	23.264	30.55	14.02
3rd Quartile	218.82	36.525	45.10	17.40
Maximum	296.40	49.600	114.00	27.00

The above table of results displays the five-number summary plus the mean for each of the four variables. The five-number summary includes minimum, first quartile (25th percentile), median, third-quartile (75th percentile) and the maximum. Examining such statistics will help us gain a general understanding of the key characteristics of the variables in the data, usually the first step in any predictive analytic exercise. The five-number summary tells us that sales ranges from 1.6 to 27 (1,600 to 27,000 units of the products). The TV advertising budget ranges from 0.7 to 296.40 ($700 to $296,400). The center of the distribution of sales is 12.9 (median) and 14.02 (mean). And so on.

The difference between the mean and the median of a numeric variable can provide an indication of skewness (or the lack of) of that variable. If the median and the mean of a variable are similar in size, that is an indication that the variable is more or less symmetric. For sales, the mean (14.02) and the median (12.9) are roughly equal, which is an indication of a symmetric distribution, or least one that has no clear sign of skewness. Similar observations can be made for TV and radio. For newspaper, the mean (30.55) and the median (25.75) are quite far apart with the mean substantially larger than the median, a sign of a right skewness in the distribution. The next step is to use graphical displays of these four variables to further confirm what we see in the numerical summary.

Since the four variables in the data are numeric variables, both histograms and boxplots can help us visualize their distributions. The following R code produces the histograms of the four numeric variables.

library(ggplot2)
library(gridExtra)
p1 <- ggplot(advert, aes(x=TV)) + geom_histogram()
p2 <- ggplot(advert, aes above imag(x=radio)) + geom_histogram()
p3 <- ggplot(advert, aes(x=newspaper)) + geom_histogram()
p4 <- ggplot(advert, aes(x=sales)) + geom_histogram()
grid.arrange(p1, p2, p3, p4, ncol = 2)

For sales, TV and radio, the histograms are roughly symmetric. There is no any clear sign of skewness in these three histograms. The histogram for newspaper does show a noticeable right skew. For a comparison, we now take a look at the boxplots. The following R code produces the boxplots of the four numeric variables.

p1 <- ggplot(advert, aes(y=TV)) + geom_boxplot()
p2 <- ggplot(advert, aes(y=radio)) + geom_boxplot()
p3 <- ggplot(advert, aes(y=newspaper)) + geom_boxplot()
p4 <- ggplot(advert, aes(y=sales)) + geom_boxplot()
grid.arrange(p1, p2, p3, p4, ncol = 2)

The box in a boxplot indicates the locations of the first quartile (the bottom edge of the box), median (the line inside the box) and third quartile (the upper edge of the box). If the line inside is roughly at the middle, that is an indication of a approximately symmetric distribution. For TV and radio, the median lines are roughly in the middle of the boxes. For sales, it does appear that the median line is closer to the bottom of the box than to the top of the box, indication of a slight right skew. However, the histogram of sales, though giving a tiny hint of a right skew, looks largely symmetric and bell-shaped. Thus, using a linear regression to predict sales seems appropriate. Note that a linear model requires that the random errors are normally distributed. Hence the response values should be normally distributed, at least approximately.

In the boxplot for newspaper, the median line is closer to the bottom edge than the top edge. Here, the right skew is more pronounced even though the boxplot gives the appearance of a less pronounced skew (the scale of the vertical axis may be distorting the picture). For example, the difference between the median and 1st quartile is 7, while the difference between the 3rd quartile and the median is 19.35, almost three times as large as the distance of 7. The same comparison for sales is 2.52 and 4.5. In comparison the skew in newspaper is more pronounced. For the discussion here, we keep the variable newspaper as is. In a practical application of linear regression, we can symmetrize newspaper by taking a log transformation.

Bivariate Data Exploration

Since we would like to examine how the predictor variables (TV, radio, newspaper) are related to the response (sales), we now explore the bivariate relationships. Since all four variables are numeric, correlation and scatterplots are appropriate tools to explore such relationships (or the lack of).

The following R code produces the scatterplots and correlation matrix.

p1 <- ggplot(advert, aes(x=TV,y=sales)) + geom_point() +
……….geom_smooth(method = “lm”, se = FALSE)
p2 <- ggplot(advert, aes(x=radio,y=sales)) + geom_point() +
……….geom_smooth(method = “lm”, se = FALSE)
p3 <- ggplot(advert, aes(x=newspaper,y=sales)) + geom_point() +
……….geom_smooth(method = “lm”, se = FALSE)
grid.arrange(p1, p2, p3, ncol = 2)
…………………
round(cor(AdData),4)

The following table shows the correlation matrix.

…….	TV……	Radio….	Newspaper..	Sales
TV	1.0000	0.0548	0.0566	0.7822
Radio	0.0548	1.0000	0.3541	0.5762
Newspaper	0.0566	0.3541	1.0000	0.2283
Sales	0.7822	0.5762	0.2283	1.0000

Each scatterplot includes a fitted regression line. In the scatterplot for TV versus sales, there is a clear positive association. The data points tend to cluster quite closely to the regression line, even though the data points with larger TV budgets are more spread out from the regression line. The scatterplot for radio versus sales also shows a fairly clear positive association, though the data points tend to deviate further from the regression line in comparison to the one for TV versus sales. In the scatterplot for newspaper versus sales, there is a hint of association, but seems quite weak. In the correlation matrix, we focus on either the last column or the bottom row. We see that the correlation between sales and TV is the strongest at 0.7822. The correlation between radio and sales is moderate (0.4237), confirming the moderate association we see in the scatterplot. The correlation between newspaper and sales is 0.1425, confirming the weak association we see in the scatterplot.

Model Fitting

For linear regression (simple or multiple), the function for model fitting is the lm() function (standing for linear models).

The following R code fits three simple linear models, one for advertising medium as the predictor.

SLR.TV <- lm(sales ~ TV, advert)
summary(SLR.TV)
…………………
SLR.Radio <- lm(sales ~ radio, advert)
summary(SLR.Radio)
…………………
SLR.Newspaper <- lm(sales ~ newspaper, advert)
summary(SLR.Newspaper)

The following three tables display the results of the three fitted models.

Table 1 – Results for Least Squares Model for the Regression of Sales on TV Advertising

…..	Coefficient…..	Std Error…………..	t-statistic…..	p-value
Intercept	7.032594	0.4578	15.36	<2e-16 ***
TV	0.047537	0.002691	17.67	<2e-16 ***

Accompanying Information (TV versus Sales)

Residual standard error: 3.259 on 198 degrees of freedom
Multiple R-squared: 0.6119, Adjusted R-squared: 0.6099
F-statistic: 312.1 on 1 and 198 DF, p-value: < 2.2e-16

Table 2 – Results for Least Squares Model for the Regression of Sales on Radio Advertising

…..	Coefficient…..	Std Error…………..	t-statistic…..	p-value
Intercept	9.31164	0.56290	16.542	<2e-16 ***
Radio	0.20250	0.02041	9.921	<2e-16 ***

Accompanying Information (Radio versus Sales)

Residual standard error: 4.275 on 198 degrees of freedom
Multiple R-squared: 0.332, Adjusted R-squared: 0.3287
F-statistic: 98.42 on 1 and 198 DF, p-value: < 2.2e-16

Table 3 – Results for Least Squares Model for the Regression of Sales on Newspaper Advertising

…..	Coefficient…..	Std Error…………..	t-statistic…..	p-value
Intercept	12.35141	0.62142	19.88	<2e-16 ***
Newspaper	0.05469	0.01658	3.30	0.00115 **

Accompanying Information (Newspaper versus Sales)

Residual standard error: 5.092 on 198 degrees of freedom
Multiple R-squared: 0.05212, Adjusted R-squared: 0.04733
F-statistic: 10.89 on 1 and 198 DF, p-value: 0.001148

The TV Model

Of the three potential predictors of sales, TV is the one that exhibits the strongest linear association with sales. We focus on the fitted regression model with TV as the predictor of sales (the results are in Table 1). The coefficient estimates for $\beta_0$ and $\beta_1$ are 7.0326 and 0.047537, respectively. According to this least squares regression model, in the absence of TV advertising, sales are 7,033 units. Furthermore, an additional $1,000 spent on TV advertising is associated, on average, with selling approximately 47.5 additional units of the products.

Table 1 shows that the coefficient estimates for intercept and TV are very large in comparison to their standard errors. As a result, the t-statistics are also large. The probabilities of observing such large t-statistics are the p-values (the last column in Table 1). These probabilities would be practically zero if the null hypotheses are true. Indeed the two p-values in Table 1 are virtually zero. Thus, we can conclude that the $\beta_0 \ne 0$ and $\beta_0 \ne 1$ . That is, we reject the null hypothesis $H_0: \beta_0 = 0$ and the null hypothesis $H_0: \beta_1 = 0$ . With the rejection of the null hypothesis $H_0: \beta_0 = 0$ , we can conclude that in the absence of TV expenditure, sales are non-zero. With the rejection of the null hypothesis $H_0: \beta_1 = 0$ , we can conclude that there is a relationship between expenditure in TV advertising and sales.

With the rejection of the two null hypotheses, the interpretations of the model results are accurate. That is, in the absence of TV advertising, sales are expected to be round 7,033 units and that an additional $1,000 in the expenditure in TV spending will result in an additional 47.5 units in sales. To further quantify these estimates, note that the 95% confidence interval for the intercept coefficient $\beta_0$ is [5.3125, 9.6309] while the 95% confidence interval for the slope coefficient $\beta_1$ is [0.0303, 0.0568]. Since the endpoints are all positive, we can conclude that $\beta_0$ and $\beta_1$ are positive.

According to Table 1, the residual standard error (RSE) is 3.259 or 3.26. Since RSE is in the same unit as the response (sales), this means that the actual sales in each market would be expected to deviate, on average, from the regression line by about 3,260 units. This means that if the true values of the unknown coefficients $\beta_0$ and $\beta_1$ were known, any prediction of sales on the basis of TV advertising expenditure would still be off by 3,260 units on average. Is this prediction error acceptable? It depends on the context. In this case, the prediction error of 3,260 leads to a percentage error of 23% (3260/1400) where 1,400 is the mean of sales.

According to Table 1, the $R^2$ statistic is 0.6119 or 0.62. This means that about 62% of the variation in sales is explained by a linear regression on TV advertising expenditure. This is only a moderate amount of variation explained. In this least squares regression model, we are only making use of the TV predictor. Very likely we need to bring in the other two predictors (radio and newspaper) to better explain sales.
….
….

This article was published here in a companion site.

Dan Ma Visualization in R
Daniel Ma Visualization in R
Dan Ma Data Exploration
Daniel Ma Data Exploration
Dan Ma linear regression
Daniel Ma linear regression
Dan Ma simple linear regression
Daniel Ma simple linear regression
Dan Ma multiple linear regression
Daniel Ma multiple linear regression

$\copyright$ 2024 – Dan Ma

Posted: October 8, 2024

Simple linear regression – making predictions

Posted on October 7, 2024 by Dan Ma

Suppose that the simple linear regression model has been fitted to the data and that through a test of hypotheses, we have some indication of a linear association between the predictor $X$ and the response $Y$ , i.e., the the test indicates that the slope coefficient $\beta_1$ is significantly different from zero. In addition to knowing how the response $Y$ is related to the predictor $X$ , we can also use the least squares estimates $\hat{\beta_0}$ and $\hat{\beta_1}$ to predict the response $Y$ given some value of interest for the predictor. We discuss such prediction procedure. We also show that there are two sources of uncertainty associated with this prediction procedure.

An illustration of the concepts discussed here is found here using the advertising dataset used in the popular text ISLR.

Two Sources of Variability

Given a value $X_0$ of the predictor that is of interest, we would like to predict the corresponding response value, say $Y_0$ . Suppose that $\hat{\beta_0}$ and $\hat{\beta_1}$ are the coefficient estimates of the unknown coefficients $\beta_0$ and $\beta_1$ , respectively (simple linear regression model fitting is discussed here). The following is the prediction of $Y_0$ .

(1)

……….

$\hat{Y_0}=\hat{\beta_0}+\hat{\beta_1} X_0=\overline{Y}+\hat{\beta_1} (X_0-\overline{X})$

In (1), the variables are in upper case to emphasize the random nature of the prediction $\hat{Y_0}$ . Because there is variability in $\hat{Y_0}$ (in fact there are the two sources of variability), it is more natural to express (1) in upper case letters. The following is the variance of the prediction $\hat{Y_0}$ .

(2)

……….

$\displaystyle \begin{aligned} \sigma_{\text{pred}}^2&=\text{Var}(Y_0 \lvert X_0) + \text{Var}(\hat{Y_0}) \\&=\ \ \sigma^2 \ + \ \ \biggl[\frac{1}{n}+\frac{(X_0-\overline{X})^2}{\Sigma (X_i-\overline{X})^2} \biggr] \sigma^2 \end{aligned}$

……………………..

Source 1

…………..

Source 2

The first source of the variability is the variability of the random error term $\epsilon$ , whose variance is $\sigma^2$ . In the simple linear regression model, the variance of the random error term can be expressed by $\text{Var}(Y_0 \lvert X_0)$ . Because of the inherent random error $\epsilon$ , we can never predict the response $Y_0$ precisely even if we know the true values of $\beta_0$ and $\beta_1$ .

The second source of the variability is due to the error in estimating the regression line at $X_0$ . This is the intrinsic variability of $\hat{\beta_0}$ and $\hat{\beta_1}$ , represented by $\text{Var}(\hat{Y_0})$ . The variance $\text{Var}(\hat{Y_0})$ is derived below.

The Prediction Procedure

The variance $\sigma^2$ , the variance of the random error term $\epsilon$ , is estimated by the $s^2$ , the square of the residual sum of squares (RSE, see here). The following gives the standard error of the prediction $\hat{Y_0}$ .

(3)

……….

$\displaystyle \text{SE}_{\text{Pred}}= s \ \sqrt{1 + \frac{1}{n}+\frac{(X_0-\overline{X})^2}{\Sigma (X_i-\overline{X})^2}}$

The 100(1- $\alpha$ )% prediction interval for $Y_0$ is given by the following:

(4)

……….

$\displaystyle \begin{aligned} &\ \ \hat{Y_0} \pm t_{n-2,\frac{\alpha}{2}} \times \text{SE}_{\text{Pred}} \\&=(\hat{\beta_0}+\hat{\beta_1} X_0) \pm t_{n-2,\frac{\alpha}{2}} \times s \ \sqrt{1 + \frac{1}{n}+\frac{(X_0-\overline{X})^2}{\Sigma (X_i-\overline{X})^2}} \end{aligned}$

Derivation

In this section, we derive the variance of $\hat{Y_0}$ , the second source of variation discussed in (2). We show the following.

(5)

……….

$\displaystyle \text{Var}(\hat{Y_0})=\biggl[\frac{1}{n}+\frac{(X_0-\overline{X})^2}{\Sigma (X_i-\overline{X})^2} \biggr] \sigma^2$

To start, we point out the following fact.

(6)

……….

$\displaystyle \text{Cov}(S_1,S_2)=(u_1v_1+u_2 v_2 + \cdots+u_n v_n) \sigma^2$

where

(6)………. $S_1=u_1 Y_1 + u_2 Y_2 +\cdots+u_n Y_n$

(6)………. $S_2=v_1 Y_1 + v_2 Y_2 +\cdots+v_n Y_n$

(6)………. $\text{Var}(Y_i)=\sigma^2$ …… $i=1,2,\cdots,n$

(6)………. $\text{Cov}(Y_i,Y_j)=0$ …… $i \ne j$

The result in (6) implies that when we calculate the covariance of two linear combinations of uncorrelated variables, we only need to sum the variances and can remove the covariance terms. To apply this result, let $u_i=\frac{1}{n}$ and $v_i=\frac{X_i-\overline{X}}{\Sigma (X_j-\overline{X})^2}$ for all $i=1,2,\cdots,n$ . Note that $\overline{Y}=u_1 Y_1+ \cdots+u_n Y_n$ . Furthermore, $\hat{\beta_1}=v_1 Y_1+ \cdots+v_n Y_n$ (see result (14) here). The variables $Y_i$ are independent and thus uncorrelated. Use (6) to derive $\text{Cov}(\overline{Y},\hat{\beta_1})$ .

$\displaystyle \begin{aligned}\text{Cov}(\overline{Y},\hat{\beta_1})&=(u_1v_1+u_2 v_2 + \cdots+u_n v_n) \sigma^2 \\&=\biggl[\frac{X_1-\overline{X}}{\Sigma (X_i-\overline{X})^2}+\cdots+\frac{X_n-\overline{X}}{\Sigma (X_i-\overline{X})^2} \biggr] \sigma^2 \\&=\frac{(\Sigma X_i)-n \overline{X}}{n \Sigma (X_i-\overline{X})^2} \sigma^2 \\&=\frac{n \overline{X}-n \overline{X}}{n \Sigma (X_i-\overline{X})^2} \sigma^2 \\&=0 \end{aligned}$

We now derive the second source of variation.

$\displaystyle \begin{aligned}\text{Var}(\hat{Y_0})&=\text{Var}(\hat{\beta_0}+\hat{\beta_1} X_0) \\&=\text{Var}[\hat{Y}+\hat{\beta_1} (X_0-\overline{X})] \\&=\text{Var}(\overline{Y})+(X_0-\overline{X})^2 \text{Var}(\hat{\beta_1})+2 (X_0-\overline{X}) \text{Cov}(\overline{Y},\hat{\beta_1}) \\&=\text{Var}(\overline{Y})+(X_0-\overline{X})^2 \text{Var}(\hat{\beta_1}) \\&=\frac{\sigma^2}{n}+(X_0-\overline{X})^2 \frac{\sigma^2}{\Sigma (X_i-\overline{X})^2} \\&=\biggl[ \frac{1}{n}+\frac{(X_0-\overline{X})^2}{\Sigma (X_i-\overline{X})^2} \biggr] \sigma^2\end{aligned}$

This article was originally published here in a companion site.

$\text{ }$

Dan Ma prediction interval
Daniel Ma prediction interval

Dan Ma linear regression
Daniel Ma linear regression

Dan Ma simple linear regression
Daniel Ma simple linear regression

Dan Ma multiple linear regression
Daniel Ma multiple linear regression

Posted: October 7, 2024

Simple linear regression – assessing the accuracy of the model

Posted on October 6, 2024 by Dan Ma

After a simple linear regression model is fitted (model fitting is discussed here), it is important to focus on the accuracy of the estimates $\hat{\beta_0}$ and $\hat{\beta_1}$ , where the former is the estimate of the unknown intercept coefficient $\beta_0$ and the latter is the estimate of the unknown slope coefficient $\beta_1$ in the simple linear regression model. For example, we would like to determine how close the estimates $\hat{\beta_0}$ and $\hat{\beta_1}$ are to the true values $\beta_0$ and $\beta_1$ . To accomplish this goal, we compute the standard errors of $\hat{\beta_0}$ and $\hat{\beta_1}$ in order to construct confidence intervals and draw inference on unknown $\beta_0$ and $\beta_1$ (discussed here). We now would like to focus on measures of global quality of model fit. For example, if a test of hypotheses tells us that the unknown slope coefficient $\beta_1$ is significantly different from zero, we would like to quantify the extent of the relationship. To this end, we discuss the notion of residual standard error along with the notion of the coefficient of determination $R^2$ . This necessitates a discussion on three different sums of squares along with a discussion of the ANOVA table and the F-test.

To facilitate the discussion, the following is the description of the training dataset used in the model fitting.

(1)

……….

$(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)$

There are $n$ observations in (1). For each observation, the following shows that fitted value of the simple linear regression model.

(2)

……….

$\hat{y_i}=\hat{\beta_0}+\hat{\beta_1} \ x_i$

where

(3)………. $\displaystyle \hat{\beta_1}=\frac{\Sigma_{i=1}^n (x_i-\overline{x})(y_i-\overline{y})}{\Sigma_{i=1}^n(x_i-\overline{x})^2}=\frac{\Sigma_{i=1}^n x_i y_i - n \ \overline{x} \ \overline{y}}{\Sigma_{i=1}^n x_i^2 - n (\overline{x})^2}$

…………. $\hat{\beta_0}=\overline{y}-\hat{\beta_1} \ \overline{x}$

Residual Standard Error

Under the simple linear regression model, each observation is associated with a random error term $\epsilon$ that has mean zero and a constant variance $\sigma^2$ . Because of the inherent uncertainty represented by the random error term, we would not be able to perfectly predict $Y$ from $X$ even if true values of $\beta_0$ and $\beta_1$ are known. Thus, we need an estimate of the standard deviation of the random error term $\epsilon$ to get an idea of whether the fitted simple linear regression model is a good fit to the data or not. The residual standard error (RSE) is an estimate of the standard deviation of the random error term $\epsilon$ in the simple linear regression model. In other words, RSE is an estimate of the $\sigma$ where $\sigma^2$ is the variance of the random error term $\epsilon$ .

(4)

……….

$\displaystyle \text{RSE} =\sqrt{\frac{\text{RSS}}{n-2}}$

…

$\displaystyle \text{RSS} =\Sigma (y_i-\hat{y_i})^2$

The quantity RSS is the residual sum of squares. RSE measures the average amount that the response will deviate from the true regression line. It represents a measure of the lack of fit of the simple linear regression model to the data. If each $y_i$ is close to the fitted value $\hat{y_i}$ , then the RSE will be small. Thus, the smaller the RSE, the better the model fit. On the other hand, a large RSE indicates that the model does not fit the data well. RSE provides an absolute measure of the lack of fit of the simple linear regression model. One drawback of using RSE is that it is measured in the same units as the response $Y$ . In many cases, it is not clear what is considered a good (small) RSE.

Sum of Squares

Unlike RSE, the $R^2$ statistic is a measure of the goodness of fit independent of the scale of the response $Y$ . Suppose we exclude the predictor $X$ from the modeling effort. Then intuitively, the sample mean $\overline{y}=\frac{y_1+\cdots+y_n}{n}$ is the best fitted value for the observed $Y_i$ . This can also be seen mathematically. If $X$ is not in the picture, then the model is $Y=\beta+\epsilon$ where $\beta$ is the unknown model coefficient and $\epsilon$ is the zero-mean random error term. We wish to find the estimate $\beta$ such that $\Sigma (Y_i-\beta)^2$ (the residual sum of squares RSS in this context) is minimized. It is straightforward to derive that this sum of squares is minimized at $\hat{\beta}=\overline{y}$ . In this case, $y_i-\overline{y}$ is the difference between a response value $Y_i$ and a fitted value under the model where only the response $Y$ is used.

In the simple linear regression model, the predictor $X$ is incorporated into the model and $\hat{y_i}=\hat{\beta_0}+\hat{\beta_1} \ x_i$ is the fitted value for the response value of $y_i$ . Then $y_i-\hat{y_i}$ is the difference between a response value and its fitted value under the simple linear regression model.

The difference $y_i-\overline{y}$ quantifies a typical size of the “error” when using this naive model of using the sample mean as the fitted value. Consider the sum of squares $\Sigma (y_i-\overline{y})^2$ , which is the sum of squared errors in the naive model. In the simple linear model, the difference $y_i-\hat{y_i}$ is the typical size of the error. Consider the sum of squares $\Sigma (y_i-\hat{y_i})^2$ , which is the sum of squared errors in using the simple linear regression model. These two sums of squares go by the following more common names.

(5)

……….

$\Sigma (y_i-\overline{y})^2$

………….

Total Sum of Squares (TSS)

(6)

……….

$\Sigma (y_i-\hat{y_i})^2$

…………

Residual Sum of Squares (RSS)

If the predictor $X$ is associated linearly with the response $Y$ , RSS is likely smaller than TSS. Then the difference TSS – RSS is the improvement of the simple linear regression model over the naive model. The larger this difference, the better the fit of the simple linear regression model. It can be shown that the difference TSS – RSS is equivalent to this sum of squares $\Sigma (\hat{Y_i}-\overline{Y})^2$ , which is known as regression sum of squares (Reg SS).

(7)

……….

$\Sigma (\hat{y_i}-\overline{y})^2$

…………

Regression Sum of Squares (RSS)

The three sums of squares are related by the following equation.

(8)

……….

$\Sigma (y_i-\hat{y_i})^2=\Sigma (y_i-\hat{y_i})^2+\Sigma (\hat{y_i}-\overline{y})^2$

(8)

…………….

TSS

…………….

RSS

…………….

Reg SS

For the proof of (8), see below. The value of Reg SS as a proportion of TSS defines the $R^2$ statistic, which is also called the coefficient of determination.

(9)

……….

$\displaystyle R^2=\frac{\text{Reg SS}}{\text{TSS}}=1-\frac{\text{RSS}}{\text{TSS}}$

The total sum of squares (TSS) measures the total variance in the response $Y$ . It can be thought of as the amount of variability inherent in the response before the regression is performed. On the other hand, RSS measures the amount of variability that is left unexplained after performing the regression. Therefore, Reg SS measures the amount of variability in the response $Y$ that is explained by performing the linear regression. As a result, the $R^2$ statistic measures the proportion of the variability in the response $Y$ that can be explained by including the predictor $X$ in the regression model.

By definition, an $R^2$ statistic is always between 0 and 1. An $R^2$ statistic that is close to 1 indicates that a large proportion of the variability in the response $Y$ has been explained by the linear regression model. An $R^2$ statistic that is near 0 indicates that the linear regression model did not explain much of the variability in the response $Y$ .

Overall, the $R^2$ statistic is a measure of the linear relationship between the single predictor $X$ and the response $Y$ . Though we discuss the $R^2$ statistic in the context of simple linear regression, the concept is applicable to multiple linear regression where there are more than one predictor that can potentially explain the response $Y$ . However, in the simple linear regression setting, the $R^2$ statistic is the same as the square of the correlation $r$ , which is defined by the following.

(10)

……….

$\displaystyle r=\frac{\Sigma (x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\Sigma (x_i-\overline{x})^2} \sqrt{\Sigma (y_i-\overline{y})^2}}$

The correlation $r$ is also a measure of linear relationship between the predictor $X$ and the response $Y$ . However, correlation cannot be extended to the multiple linear regression setting. Correlation can only quantify the association between a single predictor $X$ and the response. It cannot quantify the association between multiple predictors and the response.

ANOVA Table

The preceding discussion on the sums of squares can be concisely summarized in an ANOVA table. For simple linear regression, the ANOVA table is as follows.

Source of Variation…..	Sum of Squares…..	df…………..	Mean Square…..	F-statistic
$\text{Regression}$	$\text{Reg SS}$	$1$	$\text{Reg SS} /1$	$\displaystyle \frac{\text{Reg SS}/1}{s^2}$
$\text{Error}$	$\text{RSS}$	$n-2$	$s^2=\text{RSS}/(n-2)$
$\text{Total}$	$\text{TSS}$	$n-1$

The first column in the table indicates the source of the variation. The second column in the table displays the sums of squares. They are calculated in (6) and (7). Note that in the second column, the bottom line is the sum of the first two lines. The third column displays the degrees of freedom. For simple linear regression, the df for RSS is 1 while the df for Reg SS is $n-2$ . Again, the bottom row is the sum of the first two rows. The fourth column is for Mean Square, which is the sum of squares divided by the df. The last column shows the F-statistic, which is the ratio of the mean square of Reg SS over the mean square of RSS. The F-statistic provides a way to test the significance of the linear regression model in addition to the t-test discussed here.

The F-Test

For a simple linear regression model, we test the null hypothesis $H_0:\beta_1=0$ versus the alternative hypothesis $H_1:\beta_1 \ne 0$ . We perform the test using the F-statistic described in the above ANOVA table, which has degrees of freedom 1 and $n-2$ .

We use the F-statistic to evaluate whether Reg SS is sufficiently large enough for us to believe that there is a linear association between $X$ and $Y$ . Under the null hypothesis $H_0$ , the F-statistic has an F-distribution with 1 and $n-2$ degrees of freedom, with an expected value close to 1. An F-statistic that is sufficiently large (as supported by a small p-value) provides evidence for rejecting the null hypothesis and in favor of the alternative hypothesis, that the simple linear regression model is useful for explaining the variation in the response $Y$ .

An illustration of the concepts discussed here is found here using the advertising dataset used in the popular text ISLR.

The Proof of (8)

Start with $y_i-\overline{y}=(y_i-\hat{y_i})+(\hat{y_i}-\overline{y})$ Then take the sums of squares of both sides.

$\displaystyle \begin{aligned} \Sigma_{i=1}^n (y_i-\overline{y})^2&=\Sigma_{i=1}^n (y_i-\hat{y_i})^2+ \Sigma_{i=1}^n (\hat{y_i}-\overline{y})^2 + 2 \ \Sigma_{i=1}^n (y_i-\hat{y_i}) \cdot (\hat{y_i}-\overline{y}) \\&=\Sigma_{i=1}^n (y_i-\hat{y_i})^2+ \Sigma_{i=1}^n (\hat{y_i}-\overline{y})^2 \end{aligned}$

Note that the quantity $2 \ \Sigma_{i=1}^n (y_i-\hat{y_i}) \cdot (\hat{y_i}-\overline{y})$ works out to be zero. The following gives the derivation.

$\displaystyle \begin{aligned} \Sigma_{i=1}^n (y_i-\hat{y_i}) \cdot (\hat{y_i}-\overline{y})&=\Sigma_{i=1}^n (y_i-\hat{\beta_0}-\hat{\beta_1} x_i)(\hat{\beta_0}+\hat{\beta_1} x_i-\overline{y}) \\&=\Sigma_{i=1}^n [y_i-(\overline{y}-\hat{\beta_1} \overline{x}]-\hat{\beta_1} x_i) (y_i-(\overline{y}-\hat{\beta_1} \overline{x}+\hat{\beta_1} x_i - \overline{y}) \\&=\Sigma_{i=1}^n [y_i-\overline{y}-\hat{\beta_1} (x_i-\overline{x})] [\hat{\beta_1} (x_i-\overline{x})] \\&=\hat{\beta_1} \Sigma_{i=1}^n (x_i-\overline{x})(y_i-\overline{y})-\hat{\beta_1}^2 \Sigma (x_i-\overline{x})^2=0 \\& \end{aligned}$

Note that the last line is zero after we plug the quantity in (3) into $\hat{\beta_1}$ .

This article was originally published here in a companion site.

$\text{ }$

Dan Ma coefficient of determination
Daniel Ma coefficient of determination

Dan Ma R-squared
Daniel Ma R-squared

Dan Ma F-test
Daniel Ma F-test

Dan Ma linear regression
Daniel Ma linear regression

Dan Ma simple linear regression
Daniel Ma simple linear regression

Dan Ma multiple linear regression
Daniel Ma multiple linear regression

Posted: October 6, 2024

Simple linear regression – assessing the accuracy of the coefficient estimates

Posted on October 5, 2024 by Dan Ma

We discuss the issue of assessing the accuracy of the coefficient estimates in simple linear regression. First, a recap of the previous discussion on simple linear regression here.

A Recap

We examine the relationship between a response variable $Y$ and a single predictor $X$ . For the purpose of this discussion, assume both are numeric variables. We assume that the relationship between $Y$ and $X$ takes the form of $Y=f(X)+\epsilon$ where $f$ is some unknown function containing systematic information that the predictor $X$ provides about the response $Y$ and $\epsilon$ is a random error term with mean zero. In simple linear regression, the function $f$ providing the systematic information takes the form of a linear function $f(X)=\beta_0+\beta_1 X$ . We consider the following modeling scenario

(1)

……….

$Y=\beta_0+\beta_1 X + \epsilon$

where $\beta_0$ is the intercept coefficient, $\beta_1$ is the slope coefficient and $\epsilon$ is the zero-mean random error term. The random error terms across different values of $X$ are independent. The equation in (1) is called a population regression line. Our goal is to obtain the sample (or estimated) regression line.

(2)

……….

$\hat{y}=\hat{\beta_0}+\hat{\beta_1} x$

where $\hat{\beta_0}$ is the coefficient estimate of the intercept coefficient $\beta_0$ and $\hat{\beta_1}$ is the coefficient estimate of the slope coefficient $\beta_1$ . The following coefficient estimates are derived in this here.

(3a)

……….

$\displaystyle \hat{\beta_1}=\frac{\Sigma_{i=1}^n (x_i-\overline{x})(y_i-\overline{y})}{\Sigma_{i=1}^n(x_i-\overline{x})^2}$

….

$\displaystyle \hat{\beta_1}=\frac{\Sigma_{i=1}^n x_i y_i - n \ \overline{x} \ \overline{y}}{\Sigma_{i=1}^n x_i^2 - n (\overline{x})^2}$

(3b)………. $\displaystyle \hat{\beta_0}=\overline{y}-\hat{\beta_1} \ \overline{x}$

The coefficient estimates in (3a) and (3b) are based on the sample data $(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)$ . The regression line as described by (3a) and (3b) is an optimal estimator of the population regression line in (1) in the sense that $\hat{\beta_0}$ and $\hat{\beta_1}$ satisfy the least squares criterion. In other words, the coefficient estimates $\hat{\beta_0}$ and $\hat{\beta_1}$ are obtained in such a way that the sum of squared errors (called the residual sum of squares or RSS) is minimized (see here for more information). When the coefficient estimates $\hat{\beta_0}$ and $\hat{\beta_1}$ are obtained according to the least squares criterion, the regression line (2) is called the least squares regression line and the resulting coefficient estimates in (3a) and (3b) are called the least squares estimates. The least squares slope estimate $\hat{\beta_1}$ is particularly important since it quantifies the directional influence of the predictor $X$ on the response variable $Y$ .

The Sampling Distributions

With the coefficient estimates $\hat{\beta_0}$ and $\hat{\beta_1}$ established, we would like to get a sense of how close $\hat{\beta_0}$ and $\hat{\beta_1}$ are to the true values of the population parameters $\beta_0$ and $\beta_1$ . To this end, we would like to discuss some behind the scene math, more specifically the sampling distributions of $\hat{\beta_0}$ and $\hat{\beta_1}$ .

Consider $\hat{\beta_1}$ , the least squares estimator of the true value of the unknown slope coefficient in the simple linear model. Note that $\hat{\beta_1}$ is constructed from training data and hence is random in nature in the sense that using a different training dataset will result in a different value of $\hat{\beta_1}$ . Therefore, $\hat{\beta_1}$ can be considered a random variable. This means that it makes sense to discuss the probability distribution of $\hat{\beta_1}$ . This probability distribution is called the sampling distribution of $\hat{\beta_1}$ . As a result, notions that are usually associated with a probability distribution, e.g., notions such as expectation, variance, bias and standard error, make sense for $\hat{\beta_1}$ . We can say the same for the sampling distribution of $\hat{\beta_0}$ .

We now define the notion of bias. Let $\hat{\theta}$ be an estimator of the true population parameter $\theta$ (but a hypothetical one). The bias of the estimator $\hat{\theta}$ is the difference $E(\hat{\theta})-\theta$ . If this difference is zero, i.e., $E(\theta)=\theta$ , then $\hat{\theta}$ is said to be an unbiased estimator of the parameter $\theta$ . Otherwise $\hat{\theta}$ is said to be a biased estimator of $\theta$ . If $\hat{\theta}$ is an unbiased estimator of $\theta$ , then it means that the estimator $\hat{\theta}$ does not systematically over or under estimate the true population mean $\theta$ .

One handy example of an unbiased estimator is the sample mean $\hat{\mu}=(x_1+\cdots+x_n)/n$ . The sample mean is an unbiased estimator of the population mean since $E(\hat{\mu})=\mu$ . It follows that the sample mean $\hat{\mu}$ does not systematically over or under estimate the true population mean $\mu$ .

The property of unbiasedness holds true for the least squares estimators $\hat{\beta_0}$ and $\hat{\beta_1}$ given in (3a) and (3b) above. We show below that $E(\hat{\beta_1})=\hat{\beta_1}$ and $E(\hat{\beta_0})=\hat{\beta_0}$ . If we estimate $\beta_0$ and $\beta_1$ with one particular dataset, the estimate will not be exactly equal to the true values of $\beta_0$ and $\beta_1$ and can substantially deviate from the true value. However, if we make the estimates using a large number of datasets, then the average of the estimates would be spot on.

Another useful notion is that of standard error. First, we derive $Var(\hat{\beta_0})$ and $Var(\hat{\beta_1})$ , the variances of $\hat{\beta_0}$ and $\hat{\beta_1}$ , respectively. The standard deviations $\sqrt{Var(\hat{\beta_0})}$ and $\sqrt{Var(\hat{\beta_1})}$ will give us a sense of how close the estimators $\hat{\beta_0}$ and $\hat{\beta_1}$ are to the true values $\beta_0$ and $\beta_1$ . However, the standard deviations contain $\sigma$ , an unknown coefficient, where $\sigma^2$ is the variance of the random error term in the simple linear regression model (1). We can estimate the standard deviations by estimating $\sigma^2$ using data, the resulting estimated standard deviations are called standard errors.

We derive below the standard errors of $\hat{\beta_0}$ and $\hat{\beta_1}$ . Given a training dataset, we obtain the point estimates for $\hat{\beta_0}$ and $\hat{\beta_1}$ . Just having the point estimates is not informative since the estimates would vary from sample to sample. That’s why we need to accompany the point estimates with the standard errors to quantify the uncertainty. The smaller of their standard errors, the closer the coefficient estimates $\hat{\beta_0}$ and $\hat{\beta_1}$ will be to the true values of $\beta_0$ and $\beta_1$ . In other words, the smaller the standard errors of $\hat{\beta_0}$ and $\hat{\beta_1}$ , the more likely that the estimates, produced from training data, will be accurate representations of the true values of $\hat{\beta_0}$ and $\hat{\beta_1}$ . To achieve the goal of assessing the accuracy of the least squares estimates and drawing inference about the intercept coefficient $\beta_0$ and the slope coefficient $\beta_1$ , the remainder of this article is to specify the sampling distributions of $\hat{\beta_0}$ and $\hat{\beta_1}$ .

The Underlying Model Assumptions

To derive the variances and the standard errors, we need to know the sampling distributions of $\hat{\beta_0}$ and $\hat{\beta_1}$ . The simple linear regression model is described in (1) above. Since we are dealing with sampling distribution, the model in (1) needs to be updated with the observed data in mind. Hence, we have the following model.

(4)

……….

$Y_i=\beta_0+\beta_1 X_i + \epsilon_i \ \ \ \ \ i=1,2,\cdots,n$

The model reflects the $n$ data points in the training data. The $n$ instances of $Y_i$ are the realizations of the response variable $Y$ and the $n$ instances of $X_i$ are the realizations of the predictor variable $X$ . The $n$ instances of $\epsilon_i$ are the realizations of the random error term $\epsilon$ . The mathematical model can be stated either in terms of the random error terms or the conditional distributions of the response $Y_i$ given $X_i$ .

Model Assumptions for Simple Linear Regression
….1. The random error terms $\epsilon_i$ have means zero.
….2. The random error terms $\epsilon_i$ have constant variance $\sigma^2$ , i.e., it is not a function
……..of $X$ .
….3. The random error terms $\epsilon_i$ are independent.
….4. The random error terms $\epsilon_i$ follow a normal distribution.

….

Model Assumptions for Simple Linear Regression
….1. $E(Y_i \lvert X_i)=\beta_0+\beta_1 X_i$ .
….2. $Var(Y_i \lvert X_i)=\sigma^2$ .
….3. The conditional random variables $Y_i \lvert X_i$ are independent.
….4. The conditional random variables $Y_i \lvert X_i$ follow a normal distribution.

The assumptions in the preceding box say that the conditional mean of the response $Y$ conditioning on a particular value of the predictor $X$ is a linear function of $X$ and that such conditional distributions have constant variance and are independent. Furthermore, the conditional random variables of $Y$ given $X$ are normally distributed. These assumptions indicate that the response $Y$ and the predictor $X$ are modeled by a bivariate normal distribution (see here).

….

Unbiasedness and Standard Errors

Both $\hat{\beta_0}$ and $\hat{\beta_1}$ are unbiased estimators of the intercept coefficient $\beta_0$ and the slope coefficient $\beta_1$ , respectively. That is, $E(\hat{\beta_0})=\beta_0$ and $E(\hat{\beta_1})=\beta_1$ . The variances of $\hat{\beta_0}$ and $\hat{\beta_1}$ are shown to be the following.

(5a)

……….

$\displaystyle Var(\hat{\beta_1})=\frac{\sigma^2}{\Sigma (X_i-\overline{X})^2}$

(5b)………. $\displaystyle Var(\hat{\beta_0})=\sigma^2 \ \biggl[\frac{1}{n}+\frac{\overline{X}^2}{\Sigma (X_i-\overline{X})^2} \biggr]=\sigma^2 \ \biggl[\frac{\Sigma X_i^2}{n \ \Sigma (X_i-\overline{X})^2} \biggr]$

For the derivation of the property of unbiasedness and the variances, see the section below called “Derivation – Unbiasedness and Standard Errors.”

The standard deviation of the intercept coefficient estimator $\hat{\beta_0}$ is $\sqrt{Var(\hat{\beta_0})}$ . Similarly, the standard deviation of the slope coefficient estimator $\hat{\beta_1}$ is $\sqrt{Var(\hat{\beta_1})}$ . For emphasis, we display the two standard deviations below.

(6a)

……….

$\displaystyle \sqrt{Var(\hat{\beta_1})}=\frac{\sigma}{\sqrt{\Sigma (X_i-\overline{X})^2}}$

(6b)………. $\displaystyle \sqrt{Var(\hat{\beta_0})}=\sigma \ \sqrt{ \frac{1}{n}+\frac{\overline{X}^2}{\Sigma (X_i-\overline{X})^2} }=\sigma \ \sqrt{\frac{\Sigma X_i^2}{n \ \Sigma (X_i-\overline{X})^2}}$

Note that in (5a), (5b), (6a) and (6b), it is customary to express the predictor variable in lower case x since the values of the predictor variable are fixed (they are realized values of $X$ ).

The two standard deviations measure how close the two estimators are to the true values of the two model coefficients. We have some observations. One is that for $\hat{\beta_1}$ , the standard deviation is smaller when the the predictor values $X_i$ are more spread out. Intuitively, this means that we have more leverage to estimate the slope coefficient when the predictor values are more spread out. Furthermore, if $\overline{X}=0$ , then the standard deviation of $\hat{\beta_0}$ would be the same as the standard deviation of the sample mean $\hat{\mu}$ . When this happens, $\hat{\beta_0}=\overline{Y}$ .

In (5a) and (5b) as well as (6a) and (6b), we see that the variances and the standard deviations contain an unknown parameter $\sigma^2$ and $\sigma$ . In general, $\sigma^2$ or $\sigma$ is unknown, but can be estimated from data. The estimate of $\sigma$ is called the residual standard error and is given below.

(7)

……….

$\displaystyle RSE=\sqrt{\frac{RSS}{n-2}}$

…

$RSS=\Sigma (y_i-\hat{y_i})^2$

RSS is called the residual sum of squares. Thus the estimated standard errors are:

(8a)

……….

$\displaystyle SE(\hat{\beta_1})=\frac{s}{\sqrt{\Sigma (X_i-\overline{X})^2}}$

(8b)………. $\displaystyle SE(\hat{\beta_0})=s \ \sqrt{ \frac{1}{n}+\frac{\overline{X}^2}{\Sigma (X_i-\overline{X})^2} }=s \ \sqrt{\frac{\Sigma X_i^2}{n \ \Sigma (X_i-\overline{X})^2}}$

where s = RSE as defined in (7).

More on Sampling Distributions of the Estimators

Before we can draw inference on $\hat{\beta_0}$ and $\hat{\beta_1}$ , we need to nail down the sampling distributions of these estimators.

Based on the discussion in the derivation section below, both the estimators $\hat{\beta_0}$ and $\hat{\beta_1}$ are linear combinations of independent normal random variables. To further facilitate the discussion, here’s the expression of $\hat{\beta_0}$ and $\hat{\beta_1}$ as linear combinations of $Y_i$ .

….

……….

$\displaystyle \hat{\beta_1}=\Sigma w_i Y_i$

…

$w_i=\frac{X_i-\overline{X}}{\Sigma (X_i-\overline{X})^2}$

….

….………. $\displaystyle \hat{\beta_0}=\Sigma v_i Y_i$ …where $v_i=[\frac{1}{n}-\overline{X} w_i ]$ ….See (15) below.

Each $Y_i$ is a conditional distribution conditional on $X_i$ . By the model assumptions above, they are independent normal distributions. Since sums of independent normal random variables are also normally distributed, the two estimators in questions have the following normal distributions.

(9a)

……….

$\displaystyle \hat{\beta_0} \sim N \biggl[\beta_0, \frac{\Sigma X_i^2}{n \ \Sigma (X_i-\overline{X})^2} \ \sigma^2 \biggr]$

(9b)………. $\displaystyle \hat{\beta_1} \sim N \biggl[\beta_1, \frac{\sigma^2}{n \ \Sigma (X_i-\overline{X})^2} \biggr]$

The notation $U \sim N(\mu,\sigma^2)$ means that the random variable $U$ follows a normal distribution with mean $\mu$ and variance $\sigma^2$ . In particular, $N(0,1)$ refers to the standard normal distribution (mean 0 and variance 1). Both (9a) and (9b) can also be translated using the standard normal distribution.

(10a)

……….

$\displaystyle \frac{\hat{\beta_0}-\beta_0}{\sqrt{\frac{\Sigma X_i^2}{n \ \Sigma (X_i-\overline{X})^2}} \ \sigma} \sim N(0,1)$

(10b)………. $\displaystyle \frac{\hat{\beta_1}-\beta_1}{\frac{\sigma}{\sqrt{\Sigma (X_i-\overline{X})^2}}} \sim N(0,1)$

In practice, $\sigma^2$ is unknown (the variance of the random error term in the simple linear regression model). If the underlying mathematical assumptions in the SLR are correct, then $\sigma$ can be estimated by s = RSE. When $\sigma$ is replaced by s in (10a) and (10b), the resulting ratios are no longer normal distributions and are instead t-distributions with $n-2$ degrees of freedom, the same number of degrees of freedom as the RSS. We have the following information about the t-distribution.

(11a)

……….

$\displaystyle t=\frac{\hat{\beta_0}-\beta_0}{\sqrt{\frac{\Sigma X_i^2}{n \ \Sigma (X_i-\overline{X})^2}} \ s} \sim t(n-2)$

(11b)………. $\displaystyle t=\frac{\hat{\beta_1}-\beta_1}{\frac{s}{\sqrt{\Sigma (X_i-\overline{X})^2}}} \sim t(n-2)$

(11a) and (11b) can be used to construct confidence intervals and hypothesis testing for the unknown coefficients $\beta_0$ and $\beta_1$ . This is discussed in the next section.

Inference on the Slope and Intercept Estimators

We show how to construct confidence intervals and perform hypothesis testing on the model parameters $\beta_0$ and $\beta_1$ . The t-distributions specified in (11a) and (11b) is used.

……………..

$100(1-\alpha)$

$\beta_0$

(12a)

……….

$\displaystyle \hat{\beta_0} \ \pm \ t_{\frac{\alpha}{2},n-2} \cdot SE(\hat{\beta_0})$

…………….. $100(1-\alpha)$ % confidence interval for $\beta_1$
(12b)………. $\displaystyle \hat{\beta_1} \ \pm \ t_{\frac{\alpha}{2},n-2} \cdot SE(\hat{\beta_1})$

Note that $c=t_{\frac{\alpha}{2},n-2}$ is the value in the random variable $T$ following a t-distribution with $n-2$ degrees of freedom such that $P(T \le c)=1-\frac{\alpha}{2}$ . In other words, $t_{\frac{\alpha}{2},n-2}$ is the $(100-\frac{\alpha}{2})$ th percentile of the t-distribution with $n-2$ degrees of freedom.

We now turn to the topic of hypothesis testing. For simple linear regression, the most common hypothesis test involves testing the following null hypothesis $H_0$ versus the following alternative hypothesis $H_1$ .

(13a)

……….

$\displaystyle H_0:$

$X$

$Y$

(13b)

……….

$\displaystyle H_1:$

$X$

$Y$

Mathematically, this corresponds to testing the following hypotheses.

(13c)

……….

$\displaystyle H_0:\beta_1=0$

(13d)

……….

$\displaystyle H_1:\beta_1 \ne 0$

The justification for (13c) versus (13d) is that if $\beta_1=0$ , then the model (1) reduces to $Y=\beta_0+\epsilon$ , which means that there is no association between $X$ and $Y$ . If we reject (13c) in favor of (13d), we believe that the predictor $X$ has an effect on the response $Y$ and the effect is quantified by the estimate $\hat{\beta_1}$ . Under the assumption that the null hypothesis $H_0:\beta_1=0$ holds, the statistic $t$ in (11b) follows a t-distribution with $n-2$ degrees of freedom.

The Reasoning Process in Tests of Hypotheses

The starting point in the reasoning process in the hypothesis testing is the null hypothesis. In this case, we assume that there is no relationship between the predictor $X$ and the response $Y$ .
We then evaluate the observed data. In this case, the observed data looks like this: $(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)$ . From the observed data, we calculate $\hat{\beta_1}$ according to (3a). Then compute the t-statistic $t$ according to (11b) with $\beta_1=0$ . The question to ask is this. What would be the likelihood of observing this $t$ if there is no relationship between $X$ and $Y$ ?In other words, we calculate the p-value. The p-value is the probability of observing a t-statistic as extreme as or more extreme than the observed $t$ .
If the calculated p-value is small (typically less than 0.05 or less than 0.01), we reject the initial belief of there being no relationship between $X$ and $Y$ (null hypothesis), rather than believing that the observed t-statistic is due to random chance. A small p-value provides strong evidence for rejecting the null hypothesis $H_0$ in favor of the alternative $H_1$ . If the p-value is not small, then we do not reject the null hypothesis $H_0$ . A larger p-value does not provide sufficient evidence for us to believe that the null hypothesis is in jeopardy. In this example, a larger p-value would not provide sufficient evidence to not believe that there is no relationship between $X$ and $Y$ .

We now elaborate on the p-value, which is the probability of observing another value of the t-statistic as extreme as or more extreme than the observed $t$ . More specifically, the p-value is the probability (according to t-statistic context)

The p-value

$\displaystyle \begin{aligned} P \biggl[ \lvert t(n-2) \lvert > \lvert t \lvert \biggr]&=P \biggl[ t(n-2)<-\lvert t \lvert \biggr]+P \biggl[t(n-2) > \lvert t \lvert \biggr] \\&=2 \times P \biggl[t(n-2) > \lvert t \lvert \biggr] \end{aligned}$

The p-value is a measure, on the scale from 0 to 1, of the strength of the evidence against $H_0$ in favor of $H_1$ . The closer the p-value is to 0, the stronger evidence against the null hypothesis $H_0$ . When the p-value is small (typically less than 0.05 or less than 0.01), we reject the null hypothesis $H_0:\beta_1=0$ and in favor of the alternative hypothesis $H_1:\beta_1 \ne 0$ . Otherwise we do not reject the null hypothesis.

Let’s consider an example. Let’s say the sample size $n$ is 30 (28 degrees of freedom) and we derive an estimate of $\hat{\beta_1}$ such that the value of $t$ is 2.1. Looking up in a t-table, the p-value is less than 0.05. Using a TI-84 calculator, the p-value is 0.02243. This probability is calculated assuming the null hypothesis (no relationship between $X$ and $Y$ ). An observed t-statistic of 2.1 is very unlikely to occur if there is indeed no relationship between the predictor $X$ and the response $Y$ . In this scenario, we have evidence to reject the null hypothesis and to support the notion that there is a relationship between $X$ and $Y$ .

Let’s consider another scenario where the estimate $\hat{\beta_1}$ results in a t-statistic of $t$ = 1.3 (sample size is also 30). The p-value is now 0.102. This p-value does not provide sufficient evidence to reject the null hypothesis of there being no association between $X$ and $Y$ .

The above discussion is for a two-sided test, $H_0:\beta_1=0$ versus $H_1:\beta_1 \ne 0$ . If we wish to perform a one-sided test, then compute the one-sided p-value and reach the conclusion accordingly. The test of (13c) versus (13d) with $\beta_1=0$ is the most common test for simple linear regression since it is natural to want to know whether there is a real linear relationship between the predictor and the response. If we wish to test $H_0:\beta_1=b$ versus $H_1:\beta_1 \ne b$ for $b \ne 0$ , we simply calculate the t-statistic of (11b) accordingly and evaluate the resulting p-value.

An illustration of the concepts discussed here is found here using the advertising dataset used in the popular text ISLR.

Derivation – Unbiasedness and Standard Errors

In this section, we derive the property of unbiasedness of slope estimator $\hat{\beta_1}$ and the intercept estimator $\hat{\beta_0}$ , along with their variances. We first focus on the sampling distribution of $\hat{\beta_1}$ , the estimator of the slope coefficient. Starting with (3a), we have the following derivation.

(14)

……….

$\displaystyle \begin{aligned} \hat{\beta_1}&=\frac{\Sigma (X_i-\overline{X}) (Y_i-\overline{Y})}{\Sigma (X_i-\overline{X})^2} \\&=\frac{\Sigma (X_i-\overline{X}) Y_i}{\Sigma (X_i-\overline{X})^2} \\&=\Sigma \biggl[\frac{X_i-\overline{X}}{\Sigma (X_i-\overline{X})^2} \biggr] Y_i =\Sigma w_i Y_i \end{aligned}$

where $w_i=\frac{X_i-\overline{X}}{\Sigma (X_i-\overline{X})^2}$ . Thus, the estimator $\hat{\beta_1}$ is a linear combination of $n$ many independent random variables. Note that $\Sigma w_i=0$ and that $\Sigma w_i X_i=1$ . The following shows that $\hat{\beta_1}$ is unbiased. (Note that $Y_i$ in $\Sigma w_i Y_i$ is a conditional random variable. Thus, $E(Y_i)$ below is a conditional expectation).

$\displaystyle \begin{aligned} E(\hat{\beta_1})&=\Sigma w_i E(Y_i) \\&=\Sigma w_i (\beta_0+\beta_1 X_i) \\&=\beta_0 \Sigma w_i + \beta_1 \Sigma w_i X_i \\&=\beta_1 \end{aligned}$

Thus, $\hat{\beta_1}$ is an unbiased estimator of the slope coefficient $\beta_1$ in the simple linear regression model. The following gives the variance of $\hat{\beta_1}$ .

$\displaystyle \begin{aligned}Var(\hat{\beta_1})&=\Sigma w_i^2 \ Var(Y_i) \\&=\sigma^2 \ \Sigma w_i^2 \\&=\sigma^2 \ \frac{\Sigma (X_i-\overline{X})^2}{[ \Sigma (X_i-\overline{X})^2]^2} \\&=\frac{\sigma^2}{ \Sigma (X_i-\overline{X})^2} \end{aligned}$

We now derive the results for the intercept estimator $\hat{\beta_0}$ . The following shows that the estimator $\hat{\beta_0}$ is also a sum of independent normal random variables.

(15)

……….

$\displaystyle \begin{aligned}\hat{\beta_0}&=\overline{Y}-\hat{\beta_1} \overline{X} \\&=\frac{1}{n} \ \Sigma Y_i-\overline{X} \ \Sigma w_i Y_i \ \ \text{where } w_i=\frac{X_i-\overline{X}}{\Sigma (X_i-\overline{X})^2} \\&=\Sigma \biggl[\frac{1}{n}- \overline{X} w_i \biggr] Y_i \\&=\Sigma v_i Y_i \ \ \text{where } v_i=[\frac{1}{n}-\overline{X} \ w_i] \end{aligned}$

To see the unbiasedness, note that $\overline{Y}=\beta_0+\beta_1 \overline{X}+\overline{e}$ where $\overline{e}=\frac{\Sigma e_i}{n}$ . Consider the following derivation.

$\displaystyle \begin{aligned} \hat{\beta_0}&=\overline{Y}-\hat{\beta_1} \overline{X} \\&=\beta_0+\beta_1 \overline{X}+\overline{e}-\hat{\beta_1} \overline{X} \\&=\beta_0+(\beta_1-\hat{\beta_1}) \overline{X}+\overline{e} \end{aligned}$

Note that $E[(\beta_1-\hat{\beta_1})]=0$ since $\hat{\beta_1}$ is unbiased. Also note that $E(\overline{e})=0$ . The following derivation shows that $\hat{\beta_0}$ is an unbiased estimator of the intercept coefficient $\beta_0$ .

$\displaystyle \begin{aligned} E(\hat{\beta_0})&=E(\beta_0+(\beta_1-\hat{\beta_1}) \overline{X}+\overline{e})\\&=E(\beta_0)+E[(\beta_1-\hat{\beta_1}) \overline{X}]+E(\overline{e}) \\&=\beta_0+E[(\beta_1-\hat{\beta_1})] \overline{X}\\&=\beta_0 \end{aligned}$

We now derive the variance of $\hat{\beta_0}$ . Since $\hat{\beta_0}=\overline{Y}-\hat{\beta_1} \overline{X}$ , we have $Var(\hat{\beta_0})=Var(\overline{Y})-2 \overline{X} Cov(\overline{Y},\hat{\beta_1})+\overline{X}^2 Var(\hat{\beta_1})$ . The following gives the derivation of each quantity.

$\displaystyle \begin{aligned} Var(\overline{Y})&=\frac{1}{n^2} Var(\Sigma Y_i) \\&=\frac{1}{n^2} \Sigma Var(Y_i) \\&=\frac{1}{n^2} \Sigma \sigma^2 \\&=\frac{\sigma^2}{n} \end{aligned}$

$\displaystyle \begin{aligned}Cov(\overline{Y},\hat{\beta_1})&=Cov(\overline{Y},\Sigma w_i Y_i) \text{ where } w_i=\frac{X_i-\overline{X}}{\Sigma (X_i-\overline{X})^2}\\&=\displaystyle \Sigma \ [w_i Cov(\overline{Y},Y_i)] \\&=\displaystyle \Sigma \ \biggl[w_i \frac{\sigma^2}{n} \biggr] \\&=\frac{\sigma^2}{n \Sigma (X_i-\overline{X})^2} \ \Sigma (X_i-\overline{X})=0 \text{ since } \Sigma (X_i-\overline{X})=0 \end{aligned}$

$\displaystyle \begin{aligned}\overline{X}^2 Var(\hat{\beta_1})&=\sigma^2 \ \frac{\overline{X}^2}{\Sigma (X_i-\overline{X})^2} \end{aligned}$

Plugging in each quantity, the following gives $Var(\hat{\beta_0})$ .

$\displaystyle \begin{aligned}Var(\hat{\beta_0})&=\frac{\sigma^2}{n}+\sigma^2 \ \frac{\overline{X}^2}{\Sigma (X_i-\overline{X})^2} \\&=\sigma^2 \biggl[ \frac{1}{n}+\frac{\overline{X}}{\Sigma (X_i-\overline{X})^2}\biggr] \\&=\sigma^2 \biggl[ \frac{\Sigma X_i^2}{n \ \Sigma (X_i-\overline{X})^2}\biggr] \end{aligned}$

This article was originally published here in a companion site.

$\text{ }$

$\text{ }$
Dan Ma linear regression
Daniel Ma linear regression

Dan Ma simple linear regression
Daniel Ma simple linear regression

Dan Ma multiple linear regression
Daniel Ma multiple linear regression

Posted: October 5, 2024

Math and Stats for AIML and Data Science

The behind the scene math and stats in AI, machine learning and data science

Generalized F-test

Ordinary least squares estimators are BLUE

Statistical significance in SLR but not in MLR

Multiple linear regression – partial correlation

Multiple linear regression – adjusted R-squared

Multiple linear regression – model fitting and basic topics

Simple linear regression – example – ISLR advertising

Simple linear regression – making predictions

Simple linear regression – assessing the accuracy of the model

Simple linear regression – assessing the accuracy of the coefficient estimates