True relation vs. least-square regression line

Figure 3.3 (p. 64) compares the true relationship between \(X\) and \(Y\) with the estimated coefficients of simple LR for a simulated data set:

• Left: The red line represents the true relationship, \(Y = f(X) = 2 + 3X\), which is known as the population regression line. The blue line is the least squares line; it is the least squares estimate for \(f(X)\) based on the observed data, shown in black.

• Right: The population regression line is again shown in red, and the least squares line in dark blue. In light blue, ten least squares lines are shown, each computed on the basis of a separate random set of observations. Each least squares line is different, but on average, the least squares lines are quite close to the population regression line.

Just like many sample means \(\hat{\mu}\) provide unbiased estimates of the true population mean \(\mu\)), the estimated coefficients \(\hat{\beta_0}\) and \(\hat{\beta_1}\) provide unbiased estimates of the true regression line in the population. When drawing many samples, the average of the estimates approximate the true population values.

Measures of deviation

To quantify the deviation between population values and sample-based estimates, we can compute standard errors (SE), residual standard errors (RSE), and confidence intervals (CI).

For simple LR, the 95%-confidence interval for the population (or true) coefficients (\(beta_0\) and \(beta_1\)) can be computed from the estimated (or sample-based) coefficients (\(\hat{\beta_0}\) and \(\hat{\beta_1}\)) and their corresponding standard errors as:

\[ \hat{\beta_0} \pm 2 \cdot SE(\hat{\beta_0}) \] \[ \hat{\beta_1} \pm 2 \cdot SE(\hat{\beta_1}) \]

Hypothesis testing

Endowed with estimated coefficients and standard errors, we can perform hypothesis tests on the coefficients:

• \(H_a\): There is a relationship between \(X\) and \(Y\): \(\beta_1 \neq 0\)
• \(H_0\): There is no relationship between \(X\) and \(Y\): \(\beta_1 = 0\)

Using a \(t\)-distribution (to assess the likelihood of an estimated coefficient’s deviation from 0) yields a \(p\)-value: The probability that the observed effect was caused exclusively by chance (i.e., \(P(data | H_0)\)) — which is not the probability that \(H_a\) is true!

For our model slm_1 (from above):

# Simple linear regression model: ---- 
slm_1 <- lm(sales ~ TV, data = ad) 

# Summary of the slm_1 model (above):
summary(slm_1)
## 
## Call:
## lm(formula = sales ~ TV, data = ad)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.3860 -1.9545 -0.1913  2.0671  7.2124 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.032594   0.457843   15.36   <2e-16 ***
## TV          0.047537   0.002691   17.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 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

As the SE-values for both coefficients are very small, the estimated (or sample-based) coefficients (\(\hat{\beta_0}\) and \(\hat{\beta_1}\)) deviate from 0.

This is tested by the \(t\)-statistic (and corresponding \(p\)-values).

Rejecting \(H_0\) for both \(\beta_0\) and \(\beta_1\) implies:

• \(\beta_0 \neq 0\): There are sales in the absence of TV advertising.
• \(\beta_1 \neq 0\): There is a relationship between TV advertisements and sales.

Residual standard error (RSE)

The residual standard error (RSE) estimates the standard deviation of the irreducible error \(\epsilon\) to quantify the extent to which our predictions would miss on average if we knew the true LR line (i.e., the population values \(\beta_0\) and \(\beta_1\)). The residual standard error (RSE) is based on the residual sums of squares (RSS, see the formulas on p. 69):

\[\text{RSS} = \sum_{i=1}^{n}{(y_i - \hat{y_i})^2}\]

and computed as:

\[\text{RSE} = \sqrt{\frac{1}{n-2} \cdot \text{RSS}}\]

As RSE quantifies the lack of fit of a model in units of \(Y\), lower values are better (ideally, \(RSE \rightarrow 0\)). In our above example model slm_1, \(RSE = 3.26\) and \(Y\) represents sales in units of 1,000. Thus, even if our model was perfectly correct, our predictions would deviate from the true values by approximately 3,260 units, on average.

Proportion of variance explained (\(R^2\))

Whereas RSE provides an absolute measure of model fit (in units of the criterion \(Y\)), the \(R^2\) statistic provides an alternative measure of fit as the proportion of non-erroreous squared deviations, or the proportion of variance explained:

\[R^2 = \frac{\text{TSS} - \text{RSS}}{\text{TSS}}\] with

\[\text{TSS} = \sum_{i=1}^{n}{(y_i - \bar{y_i})^2}\]

(Note that TSS is a special case of RSS, with only a single predicted value \(\bar{y_i}\).)

In contrast to error measures, \(R^2\) ranges from \(0\) to \(1\), with higher values indicating better model fits (ideally, \(R^2 \rightarrow 1\)).

Given a model, \(R^2\) measures the proportion of variability in \(Y\) that can be explained by \(X\). An \(R^2\) value close to 1 indicates that a large proportion of the variability in the response variable is explained by the regression model.

Note:

• TSS is a special case of RSS (with only )

• Although \(R^2\) is always between 0 and 1, it can still be challenging to determine what constitutes a “good” \(R^2\) value. For domains in which the true model is highly linear and \(\epsilon\) is small, we can obtain \(R^2 \approx 1\). In social sciences, relationships can be non-linear and subject to many other influences, even small \(R^2 > 0\) can be meaningful.

• The \(R^2\) is related to the correlation \(r = \text{Cor}(X, Y)\). In simple LR (with only 1 predictor), \(R^2 = r^2\).

Accuracy measures in R

We can easily define R functions for the total sum of squares (TSS), residual sum of squares (RSS), the residual standard error (RSE), and the \(R^2\)-statistic, when providing both the true and predicted values as arguments:

# Measures of model accuracy: 
# (see @JamesEtAl2021, p. 69ff.): 


# Total sum of squares (TSS): ---- 
TSS <- function(true_vals, pred_vals = NA){
  
  # Special case: 
  if (is.na(pred_vals)){
    pred_vals <- mean(true_vals)
  }
  
  sum((true_vals - pred_vals)^2)
}


# Residual sum of squares (RSS): ---- 
RSS <- function(true_vals, pred_vals){
  
  sum((true_vals - pred_vals)^2)
}


# Residual standard error (RSE): ----
RSE <- function(true_vals, pred_vals){
  
  n <- length(true_vals)  
  
  rss <- RSS(true_vals, pred_vals)  # defined above
  
  sqrt(1/(n - 2) * rss)
}


# R^2 statistic: ----
Rsq <- function(true_vals, pred_vals){
  
  tss <- TSS(true_vals)  # defined above
  
  rss <- RSS(true_vals, pred_vals)  # defined above
  
  (tss - rss)/tss
}

ad 1. Overall effect vs. contributions of individual predictors

Is at least one of the predictors \(X_1,\ X_2,\ \ldots,\ X_p\) useful in predicting the response \(Y\)?

A large \(F\)-statistic (\(F > 1\), taking into account the \(df\) implied by \(n\) and \(p\)) implies that at least one of the predictors is successful in predicing the outcome/response variable.

In contrast to this overall effect, the \(t\)-statistic and the corresponding \(p\)-value of each coefficient informs about whether each predictor is related to the outcome/response when adjusting for the other predictors.

???: Note that \({t_{j}}^2 = F_j\) for the partial effect of each predictor \(X_j\) (see James et al., 2021, p. 77, Footnote 7).

Good question: Given individual \(p\)-values for each coefficient, why do we need the overall \(F\)-statistic?
Consider an example with many predictors (or multiple tests): For instance, assume a multiple LR model with \(p = 100\) predictors, but
\(H_0: \beta_1 = \beta_2 =\ \ldots\ = \beta_p = 0\) is true, so no variable is truly associated with the response. In this situation, about \(5\)% of the \(p\)-values associated with each coefficient will be below \(0.05\) by chance (and — simply by definition — “statistically significant”). In other words, we expect to see approximately five significant \(p\)-values in the absence of any true association between the predictors and the response. By contrast, the \(F\)-statistic does not suffer from the problem of multiple tests because it adjusts for the number of predictors \(p\).

ad 2. Deciding on important variables

Do all the predictors help to explain \(Y\) , or is only a subset of the predictors useful?

The task of determining which predictors are associated with the response, in order to fit a single model involving only those predictors, is referred to as variable selection (see Chapter 6). With \(p\) predictors, we can create \(2^p\) different models containing subsets of predictor variables (not yet considering interactions). Various statistics allow judging the quality of a model (including AIC, BIC).

Three classical approaches for solving the variable selection task are:

1. Forward selection: Begin with the null model and add the predictor that results in a model with the lowest RSS. Proceed until some stopping criterion is reached.

2. Backward selection: Begin with the full model and remove the predictor that has the highest \(p\) value. Proceed until some stopping criterion is reached.

3. Mixed selection: Begin with forward selection, but also remove predictors when their \(p\)-values increases above some criterion. Continue until all variables remaining in the model have a sufficiently low \(p\)-value.

Note:

• Backward selection cannot be used when \(p > n\).

• Forward selection can always be used (even when \(p > n\)), but is a greedy approach (i.e., might include variables early that later become redundant). Mixed selection can remedy this.

ad 3. Model fit

How well does the model fit the data?

Two of the most common numerical measures of model fit are the residual standard error (RSE) and the fraction of variance explained (\(R^2\)).

Note:

• Relation between \(R^2\) and correlation: In simple LR (with only 1 predictor), we had \(R^2 = r^2\). In multiple LR, we have \(R^2 = \text{Cor}(Y, \hat{Y})^2\) (i.e., multiple LR maximizes the correlation between response \(Y\) and the predictions of the fitted linear model \(\hat{Y}\)).

• Number of predictors \(p\), \(R^2\), vs. model quality: \(R^2\) always increases when adding variables to a model, even if those variables are only weakly associated with the response. This is due to the fact that adding another variable always results in a decrease in the RSS on the training data, but not necessarily the testing data.

• Graphical summaries can reveal problems with a model that are not visible from numerical statistics (e.g., non-linearity).

Example: Figure 3.5 (p. 81) shows a non-linear relationship in data for a multiple LR with 2 predictors. Such a non-linear pattern often suggests a synergy or interaction effect.

ad 4. Predictions

Given a set of predictor values, what response value should we predict, and how accurate is our prediction?

For a specific multiple LR model and a set of values for the predictors \(X_1, X_2,\ \ldots, X_p\), it is straightforward to predict the response \(\hat{Y}\). There are three sorts of uncertainty associated with this prediction:

1. Reducible error: The coefficients \(\hat{\beta_0}, \hat{\beta_1}, \ldots, \hat{\beta_p}\) are only estimations for the true population values \(\beta_0, \beta_1,\ \ldots, \beta_p\).

2. Model bias: Even the best linear model \(f(X)\) must not be the true relationship (which may be non-linear).

3. Irreducible error: Even the true population model would be off by random error/noise \(\epsilon\).

Distinction between confidence interval (average \(Y\) over many observations) vs. prediction interval (uncertainty surrounding \(Y\) for a particular observation).

• Confidence intervals express sampling uncertainty in quantities estimated from many data points. The more data, the less sampling uncertainty, and hence the thinner the interval.

• Prediction intervals, on top of the sampling uncertainty, also express uncertainty around a single value, which makes them wider than the confidence intervals.

See blog post at Confidence intervals vs prediction intervals (2021-11-08).

1. Non-additive relationships

Removing the additive assumption of standard LR models: Including interaction effects.

Example: Two quantitative predictors

Figure 3.5 (p. 81) suggested a non-linear relationship in the Advertising data for a multiple LR with 2 predictors (TV and radio). Such a non-linear pattern often suggests a synergy or interaction effect between both predictors:

# Multiple linear regression model (2 predictors + interaction): ---- 
mlm_6 <- lm(sales ~ TV + radio + TV:radio, data = ad)
mlm_6 <- lm(sales ~ TV * radio, data = ad)  # same model (main effects + interaction)

# Summary mLR:
summary(mlm_6)
## 
## Call:
## lm(formula = sales ~ TV * radio, data = ad)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3366 -0.4028  0.1831  0.5948  1.5246 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.750e+00  2.479e-01  27.233   <2e-16 ***
## TV          1.910e-02  1.504e-03  12.699   <2e-16 ***
## radio       2.886e-02  8.905e-03   3.241   0.0014 ** 
## TV:radio    1.086e-03  5.242e-05  20.727   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9435 on 196 degrees of freedom
## Multiple R-squared:  0.9678, Adjusted R-squared:  0.9673 
## F-statistic:  1963 on 3 and 196 DF,  p-value: < 2.2e-16

# Without interaction: ---- 
mlm_x <- lm(sales ~ TV + radio, data = ad)
# summary(mlm_x)

Notes:

• Testing for interactions relaxes the additivity assumption insofar as multiplicative relations are included (see explanation and examples on p. 88).

• The new \(R^2 = 96.7\%\), whereas it was only \(R^2 = 89.6\%\) without the interaction.
  Thus,68.3% of the remaining variance in sales was explained by including the interaction term.

• The hierarchical principle: If we include an interaction in a model, we should also include the main effects, even if the \(p\)-values associated with their coefficients are not significant. (In mlm_6, both the two main effects and their interaction are significant.)

Example: One qualitative and one quantitative predictor

An interaction between a qualitative variable and a quantitative variable has a particularly nice interpretation (p. 89 ). Consider the Credit data set from above, and suppose that we wish to predict balance using the income (quantitative) and student (qualitative) variables:

# head(cred)

# (a) Linear model with a quantitative and qualitative predictor (NO interaction):
mlm_7a <- lm(Balance ~ Income + Student, data = cred)
summary(mlm_7a)
## 
## Call:
## lm(formula = Balance ~ Income + Student, data = cred)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -762.37 -331.38  -45.04  323.60  818.28 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 211.1430    32.4572   6.505 2.34e-10 ***
## Income        5.9843     0.5566  10.751  < 2e-16 ***
## StudentYes  382.6705    65.3108   5.859 9.78e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 391.8 on 397 degrees of freedom
## Multiple R-squared:  0.2775, Adjusted R-squared:  0.2738 
## F-statistic: 76.22 on 2 and 397 DF,  p-value: < 2.2e-16

# (b) Linear model with a quantitative and qualitative predictor (WITH interaction):
mlm_7b <- lm(Balance ~ Income + Student + Income:Student, data = cred)
summary(mlm_7b)
## 
## Call:
## lm(formula = Balance ~ Income + Student + Income:Student, data = cred)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -773.39 -325.70  -41.13  321.65  814.04 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       200.6232    33.6984   5.953 5.79e-09 ***
## Income              6.2182     0.5921  10.502  < 2e-16 ***
## StudentYes        476.6758   104.3512   4.568 6.59e-06 ***
## Income:StudentYes  -1.9992     1.7313  -1.155    0.249    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 391.6 on 396 degrees of freedom
## Multiple R-squared:  0.2799, Adjusted R-squared:  0.2744 
## F-statistic:  51.3 on 3 and 396 DF,  p-value: < 2.2e-16

When expressed graphically, including the interaction results in a more flexible model that allows not only for different intercepts, but also different slopes for both qualitative groups:

Notes:

• Whereas including main effects allow only for parallel lines with different intercepts (for each level of the qualitative predictor), including interaction term additionally allows for different slopes (i.e., different impacts of the quantitative predictor on different qualitative groups).

• See derivation of line equations (i.e., values of intercept and slope) from coefficients (James et al., 2021, p. 90)

# Scatterplot: 
plot(x = cred$Income, y = cred$Balance, 
     main = "Balance by income (including interaction with student/no student)", 
     pch = 20, col = ifelse(cred$Student == "Yes", adjustcolor("deepskyblue", .33), adjustcolor("deeppink", .33)),  
     xlab = "Income", ylab = "Balance", xlim = c(0, 200), ylim = c(0, 1600))

# From model:
cf <- coefficients(mlm_7b)

# Line equations (with interaction of categorical predictor): 
l_inc_no_stud <- function(x){ cf[1] + cf[2] * x }
l_inc_student <- function(x){ (cf[1] + cf[3]) + (cf[2] + cf[4]) * x }

curve(expr = l_inc_student, lwd = 2, col = "deepskyblue", add = TRUE)
curve(expr = l_inc_no_stud, lwd = 2, col = "deeppink", add = TRUE)

Recreating Figure 3.7b: 1 quantitative and 1 qualitative predictor with interaction.

2. Non-linear relationships

Removing the linear assumption of standard LR models: Using polynomial regression.

The LR model assumes a linear relationship between the predictors and response. When the true relationship between the response and the predictors is non-linear, a simple way to allow for incorporating incorporating non-linear associations into a linear model is to include transformed versions of the predictors.

Extending the linear model to accommodate non-linear relationships is known as polynomial regression, as we include polynomial functions of the predictors in a LR model.

Example of polynomial regression: The points in Figure 3.8 (p. 91) seem to have a quadratic shape, suggesting that a model of the form

\(\text{mpg} = \beta_0 + (\beta_1 \cdot \text{hp}) + (\beta_2 \cdot \text{hp}^2) + \epsilon\)

may provide a better fit.

# Data: 
Auto <- ISLR2::Auto

# Variables:
hp  <- Auto$horsepower  # predictor
mpg <- Auto$mpg  # outcome

# Scatterplot:
plot(x = hp, y = mpg,
     main = "Auto data: mpg by horsepower", xlab = "horsepower", ylab = "mpg", 
     pch = 20, cex = 1, col = adjustcolor("black", .25))

# (a) Linear model with one predictor (hp):
mlm_8a <- lm(mpg ~ hp)
summary(mlm_8a)
## 
## Call:
## lm(formula = mpg ~ hp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.5710  -3.2592  -0.3435   2.7630  16.9240 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 39.935861   0.717499   55.66   <2e-16 ***
## hp          -0.157845   0.006446  -24.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.906 on 390 degrees of freedom
## Multiple R-squared:  0.6059, Adjusted R-squared:  0.6049 
## F-statistic: 599.7 on 1 and 390 DF,  p-value: < 2.2e-16

abline(mlm_8a, lwd = 2, col = "deepskyblue")

# (b) Polynomial LR model with one predictor and its square (as 2nd predictor):

hp_sq <- hp^2  # squared predictor 
mlm_8b <- lm(mpg ~ hp + hp_sq)

mlm_8b <- lm(mpg ~ hp + I(hp^2))  # Note: I() inhibits the interpretation of ^ 
                                  # (which has special meaning in an R formula)
summary(mlm_8b)
## 
## Call:
## lm(formula = mpg ~ hp + I(hp^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.7135  -2.5943  -0.0859   2.2868  15.8961 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 56.9000997  1.8004268   31.60   <2e-16 ***
## hp          -0.4661896  0.0311246  -14.98   <2e-16 ***
## I(hp^2)      0.0012305  0.0001221   10.08   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.374 on 389 degrees of freedom
## Multiple R-squared:  0.6876, Adjusted R-squared:  0.686 
## F-statistic:   428 on 2 and 389 DF,  p-value: < 2.2e-16

curve(expr = coef(mlm_8b)[1] + coef(mlm_8b)[2] * x + coef(mlm_8b)[3] * x^2, 
      lwd = 2, col = "deeppink", add = TRUE)

Figure 3.8

Notes:

• The fit of the polynomial model (mlm_8b) is superior: \(R^2 = 0.686\), compared to \(0.605\) for the linear fit (mlm_08a). Also, the \(p\)-value for the coefficient of the quadratic term is highly significant.

• In this equation, the value of mpg is a non-linear function of hp. However, the model is still a linear model — a multiple LR model with two predictors \(X_1 = \text{hp}\) and \(X_2 = \text{hp}^2\).

• We could extend the polynomial LR approach by including terms for \(\text{hp}^3\), \(\text{hp}^4\), etc. Although the resulting curves would provide closer fits, they are also prone to overfitting the model training data.

See Chapter 7 for further non-linear extensions of the linear model (e.g., basis functions, splines, local regression, generalized additive models).

ad 1. Non-linearity of the response-predictor relationships

A LR model assumes a linear (i.e., straight-line) relationship between predictors and response. If the true relationship is not linear, all conclusions drawn from a linear model are suspect.

Diagnosis:

A useful graphical tool for identifying non-linearity are residual plots:

Re-creating Figure 3.9 (p. 93):

• Simple LR model: Plotting the residual errors \(e_i = y_i − \hat{y_i}\) as a function of the the predictor \(x_i\):

# Data: 
Auto <- ISLR2::Auto

# Variables:
hp  <- Auto$horsepower  # predictor
mpg <- Auto$mpg  # outcome

# (a) Linear model with one predictor (hp):
mlm_8a <- lm(mpg ~ hp)
# summary(mlm_8a)

# Scatterplot of residuals:
mpg_residuals <- mpg - predict(mlm_8a)
plot(x = hp, y = mpg_residuals,
     main = "Residuals: mpg by horsepower", xlab = "horsepower", ylab = "mpg", 
     pch = 20, cex = 1, col = adjustcolor("black", .25))
abline(0, 0, lty = 2, lwd = .5)

# identify(x = hp, y = mpg_residuals, plot = TRUE)

# Note: Default residual plot in 
# plot(mlm_8a)

• Polynomial LR model: Plotting residual errors \(e_i = y_i − \hat{y_i}\) vs. the predictor \(x_i\):

# (b) Polynomial LR model with one predictor and its square (as 2nd predictor):
hp_sq <- hp^2  # squared predictor 
mlm_8b <- lm(mpg ~ hp + hp_sq)
# summary(mlm_8b)

# Scatterplot of residuals:
mpg_residuals_2 <- mpg - predict(mlm_8b)
plot(x = hp, y = mpg_residuals_2,
     main = "Residuals: mpg by horsepower + hp^2", xlab = "hp", ylab = "mpg", 
     pch = 20, cex = 1, col = adjustcolor("black", .25))
abline(0, 0, lty = 2, lwd = .5)

# identify(x = hp, y = mpg_residuals_2, plot = TRUE)

# Note: Default residual plot in 
# plot(mlm_8b)

• Multiple LR model: Plotting the residual errors \(e_i = y_i − \hat{y_i}\) vs. the predicted (or fitted) values .

Interpretation:

• The presence of a pattern (e.g., U-shape) may indicate a problem with some aspect of the linear model.

• If the linearity assumption is violated, a simple approach for improving a model can be to use a non-linear transformation of the predictor(s) (e.g., log(\(X\)), \(\sqrt{X}\), or \(X^2\)).

ad 2. Correlation of error terms

A LR model assumes that the error terms \(e_i = y_i − \hat{y_i}\) are uncorrelated. Violating this assumption results in unjustified confidence (e.g., larger confidence intervals).

Examples:

• If we accidentally used each of \(n\) data points twice (i.e., created a model that considered \(2n\) data points), our confidence intervals would overestimate the truly warranted confidence by a factor of \(\sqrt{2}\).

• Correlated error terms frequently occur in the context of time series data, which often exhibit positively correlated errors of adjacent data points.

Diagnosis:

Plot the residual errors from a LR model as a function of time (see Figure 3.10, p. 95):

Interpretation:

• When viewing residual errors as a function of time, we should not see any time-related trends or patterns. By contrast, tracking in the residuals (i.e., adjacent residuals having similar signs or values) suggest correlated error terms.

Examples beyond time series:

• Correlation between errors can occur when other relevant variables influence the relationship between predictor(s) and response variables. For instance, a study may aim to predict individuals’ height from their weight. If clusters in the data (e.g., families, communities, countries) consider individuals on similar diets, their error terms may be correlated. Avoiding such correlations is critical factor of (ensuring the internal validity of an) experimental design.

ad 3. Non-constant variance of error terms

A LR model assumes that error terms have a constant variance \(\text{Var}(\epsilon_i) = \sigma^2\) (aka. heteroscedasticity).

However, variances of the error terms may increase (or decrease) with the value of the response.

Diagnosis:

• A funnel shape in the residual plot indicates non-constant variances in the errors (heteroscedasticity).

• A concave transformation of the outcome/response variable \(Y\) (e.g., using \(\text{log}(Y)\) or \(\sqrt{Y}\)) can restore homoscedasticity (see Figure 3.11, p. 96):

• When the variance of response values \(Y_i\) is known, we can fit a LR model with weighted least squares, with weights proportional to the inverse variances. LR software usually allows for observation weights.

ad 4. Outliers

Definition: An outlier is a data point for which \(y_i\) differs substantially from the value \(\hat{y_i}\) predicted by the model (i.e., the value’s residual error is large). This can — but does not have to — indicate errors (e.g., in measuring or coding the original data) and should be investigated carefully.

Even when the model coefficients (and thus the best-fitting regression line) may change very little, removing outliers may dramatically improve the model fit (by reducing \(RSE\) and increasing \(R^2\)).

Example of Figure 3.12 (p. 97):

Interpretation:

• Left: The least squares regression line is shown in red, and the regression line after removing the outlier is shown in blue.

• Center: The residual plot clearly identifies the outlier.

• Right: The outlier has a studentized residual of \(+6\); typically we expect values between \(−3\) and \(+3\).

Diagnosis:

• Residual plots can be used to identify outliers. To standardize, we can plot studentized residuals that divide each residual error \(e_i\) by its estimated standard error. Observations whose studentized residuals are greater than \(|3|\) are possible outliers.

ad 5. High-leverage points

Definition: Whereas outliers are observations for which the response \(y_i\) is unusual given the predictor \(x_i\). By contrast, observations with high leverage have an unusual predictor value for \(x_i\).

Example: Observation 41 in the left-hand panel of Figure 3.13 (p. 98) has high leverage, in that the predictor value for this observation is large relative to the other observations:

Interpretation:

• Left: Observation 41 is a high leverage point, while 20 is (perhaps) an outlier. The red line is the fit to all the data, and the blue line is the fit with observation 41 removed.

• Center: The red observation is not unusual in terms of its \(X_1\) value or its \(X_2\) value, but still falls outside the bulk of the data, and hence has high leverage.

• Right: Observation 41 has a high leverage and a high residual.

Overall, a few high leverage observations can have a sizable impact on the estimated regression line.

Diagnosis:

• In order to quantify an observation’s leverage, we compute the leverage statistic: If a given observation has a leverage statistic that greatly exceeds \((p + 1)/n\), then we may suspect that the corresponding point has high leverage.

• In Figure 3.13 (p. 98), observation 41 stands out as having a very high leverage statistic as well as a high studentized residual. In other words, it is an outlier as well as a high leverage observation. By contrast, observation 20 has low leverage (i.e., not much influence on the model/line of best fit).

ad 6. Collinearity

Definition: Collinearity refers to the situation in which two or more predictor variables are closely related to one another.

When two or more predictors are collinear (i.e., related to each other), it is difficult or impossible to separate their individual effects on the response. In other words, when predictors tend to increase or decrease together, their \(\beta\)-weights can fluctuate widely, as we cannot separately measure their relation to the response.

Example: Two scatterplot for predictors in the Credit dataset (without and with collinearity between predictors):

Interpretation:

• Left: A plot of age versus limit. These two variables are not collinear.

• Right: A plot of rating versus limit. There is high collinearity.

Example: Two contour plots of RSS (without and with collinearity between predictors):

Interpretation:

• Left: A contour plot of RSS for the regression of balance onto age and limit. In the absence of collinearity between predictors, the minimum value is well defined.

• Right: A contour plot of RSS for the regression of balance onto rating and limit. In the presence of collinearity between predictors, there are many pairs (\(\beta_\text{Limit}\), \(\beta_\text{Rating}\)) with a similar RSS value.

Consequences of collinearity:

• Collinearity yields a great deal of uncertainty in the coefficient estimates.

• The standard error of coefficient estimates \(\hat{\beta_j}\) grows, resulting in a declining \(t\)-statistic, and lower power (for rejecting \(H_0\)).

Example: See Table 3.11 (p. 102) for two models without vs. with collinearity.

Diagnosis:

• To detect collinearity between pairs of predictors, inspect the correlation matrix of the predictors.

• To detect multicollinearity (i.e., relationships between three or more predictors), compute the variance inflation factor (VIF) for each predictor. VIF values are \(\geq 1\): VIF\((\hat{\beta_j}) \geq 5\) suggest problematic amounts of collinearity.

Remedy:

1. Drop one or more problematic predictors from the model.

2. Combine collinear variables into a single predictor (e.g., the average of standardized versions of the individual predictors).