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).

Loading packages

Here we load the MASS package, which is a very large collection of data sets and functions, and the ISLR2 package, which includes the data sets associated with this book:

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:ISLR2':
## 
##     Boston
## The following object is masked from 'package:dplyr':
## 
##     select
library(ISLR2)

If you receive an error message when loading any of these libraries, it likely indicates that the corresponding library has not yet been installed on your system. Some libraries, such as MASS, come with R and do not need to be separately installed on your computer.

Installing packages

However, other packages, such as ISLR2, must be downloaded the first time they are used. This can be done directly from within R (using the buttons/menus of an IDE or install.packages() function). For example, on a Windows system, select the Install package option under the Packages tab. After you select any mirror site, a list of available packages will appear. Simply select the package you wish to install and R will automatically download the package. Alternatively, this can be done at the R command line via install.packages("ISLR2").

An installation only needs to be done the first time you use a package. However, the library() function must be called within each R session (to load the package).

Predictions with confidence vs. prediction intervals

The predict() function can be used to produce confidence intervals and prediction intervals for the prediction of medv for a given value of lstat:

# values of interest: 
my_vals <- c(5, 10, 15)

# predictions with confidence intervals:
predict(lm.fit, data.frame(lstat = my_vals),
    interval = "confidence")
##        fit      lwr      upr
## 1 29.80359 29.00741 30.59978
## 2 25.05335 24.47413 25.63256
## 3 20.30310 19.73159 20.87461

# predictions with prediction intervals: 
predict(lm.fit, data.frame(lstat = my_vals),
    interval = "prediction")
##        fit       lwr      upr
## 1 29.80359 17.565675 42.04151
## 2 25.05335 12.827626 37.27907
## 3 20.30310  8.077742 32.52846

For instance, the 95,% confidence interval associated with a lstat value of 10 is \((24.47, 25.63)\), and the 95,% prediction interval is \((12.828, 37.28)\). As expected, the confidence and prediction intervals are centered around the same value (a predicted value of \(25.05\) for medv when lstat equals 10), but the prediction intervals are substantially wider.

Visualizations

We will now plot medv and lstat along with the least squares regression line using the plot() and abline() functions:

plot(x = lstat, y = medv, 
     main = "Scatterplot of response by predictor", 
     pch = 20, col = adjustcolor("black", .33))
abline(lm.fit, lwd = 2, col = "deepskyblue")

There is some evidence for non-linearity in the relationship between lstat and medv (see below).

Plotting lines

The abline() function can be used to draw any line, not just the least squares regression line. To draw a line with intercept a and slope b, we type abline(a, b). Below we experiment with some additional settings for plotting lines and points. The lwd = 3 argument causes the width of the regression line to be increased by a factor of 3; this also works for the plot() and lines() functions. We can also use the pch option to create different plotting symbols:

# Basic plot: 
plot(lstat, medv)

# Lines:
# abline(lm.fit, lwd = 3)
abline(lm.fit, lwd = 3, col = "deeppink")

# Point symbols and colors:
plot(lstat, medv, 
     pch = 20, col = adjustcolor("darkblue", alpha.f = .33))

# plot(lstat, medv, pch = 20)
plot(lstat, medv, pch = "+", col = adjustcolor("deeppink", alpha.f = .67))

# Note: Different pch values:
# plot(1:20, 1:20, pch = 1:20)

Diagnostic plots

Next we examine some diagnostic plots, several of which were discussed in Section 3.3.3. Four diagnostic plots are automatically produced by applying the plot() function directly to the output object obtained by lm(). This command will normally produce one plot at a time, and hitting Enter will generate the next plot. To view all four plots together, we can use the par() and mfrow() functions, which instruct R to split the display screen into separate panels, so that multiple plots can be viewed simultaneously. For example,par(mfrow = c(2, 2)) divides the plotting region into a \(2 \times 2\) grid of panels:

par(mfrow = c(2, 2))

# 4 standard diagnostic plots:  
plot(lm.fit, 
     pch = 20, col = adjustcolor("black", .33))  

par(opar)  # restore original par

Residual plots

Alternatively, we can compute the residuals from a linear regression fit using the residuals() function. The function rstudent() will return the studentized residuals, and we can use this function to plot the residuals against the fitted values.

# Residual plot:
plot(x = predict(lm.fit), y = residuals(lm.fit), 
     main = "Residual plot", 
     pch = 20, col = adjustcolor("black", .33))

# Studentized residual plot: 
plot(x = predict(lm.fit), y = rstudent(lm.fit), 
     main = "Studentized residuals", 
     pch = 20, col = adjustcolor("black", .33))

On the basis of the residual plots, there is some evidence of non-linearity.

Leverage statistics

Leverage statistics can be computed for any number of predictors using the hatvalues() function:

# Compute leverage statistics: 
hatvalues(lm.fit)[1:6]  # first 6 values
##           1           2           3           4           5           6 
## 0.004262518 0.002455527 0.004863680 0.005639779 0.004058706 0.004127513

# Plotting leverage for all values:
plot(hatvalues(lm.fit), 
     main = "Leverage diagnostics", 
     pch = 20, col = adjustcolor("black", .33)) 

# Identify max value:
which.max(hatvalues(lm.fit))
## 375 
## 375

The which.max() function identifies the index of the largest element of a vector. In this case, it tells us which observation has the maximum leverage statistic.

Checking for collinearity

To check for potential collinearity between predictors, the vif() function of the car package can be used to compute variance inflation factors.¹ Most VIF’s are low to moderate for this data:

# install.packages(car)  # if car has not been installed yet
library(car)
vif(lm.fit)
##     crim       zn    indus     chas      nox       rm      age      dis 
## 1.792192 2.298758 3.991596 1.073995 4.393720 1.933744 3.100826 3.955945 
##      rad      tax  ptratio    black    lstat 
## 7.484496 9.008554 1.799084 1.348521 2.941491

Excluding variables and updating models

What if we would like to perform a regression using all of the variables but one? For example, in the above regression output, age has a high \(p\)-value. So we may wish to run a regression excluding this predictor. The following syntax results in a regression using all predictors except age:

lm.fit_2 <- lm(medv ~ . - age, data = Boston)
summary(lm.fit_2)
## 
## Call:
## lm(formula = medv ~ . - age, data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.6054  -2.7313  -0.5188   1.7601  26.2243 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.436927   5.080119   7.172 2.72e-12 ***
## crim         -0.108006   0.032832  -3.290 0.001075 ** 
## zn            0.046334   0.013613   3.404 0.000719 ***
## indus         0.020562   0.061433   0.335 0.737989    
## chas          2.689026   0.859598   3.128 0.001863 ** 
## nox         -17.713540   3.679308  -4.814 1.97e-06 ***
## rm            3.814394   0.408480   9.338  < 2e-16 ***
## dis          -1.478612   0.190611  -7.757 5.03e-14 ***
## rad           0.305786   0.066089   4.627 4.75e-06 ***
## tax          -0.012329   0.003755  -3.283 0.001099 ** 
## ptratio      -0.952211   0.130294  -7.308 1.10e-12 ***
## black         0.009321   0.002678   3.481 0.000544 ***
## lstat        -0.523852   0.047625 -10.999  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.74 on 493 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7343 
## F-statistic: 117.3 on 12 and 493 DF,  p-value: < 2.2e-16

Alternatively, the update() function can be used on an existing model fit:

lm.fit_3 <- update(lm.fit, ~ . - age)
summary(lm.fit_3)
## 
## Call:
## lm(formula = medv ~ crim + zn + indus + chas + nox + rm + dis + 
##     rad + tax + ptratio + black + lstat, data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.6054  -2.7313  -0.5188   1.7601  26.2243 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.436927   5.080119   7.172 2.72e-12 ***
## crim         -0.108006   0.032832  -3.290 0.001075 ** 
## zn            0.046334   0.013613   3.404 0.000719 ***
## indus         0.020562   0.061433   0.335 0.737989    
## chas          2.689026   0.859598   3.128 0.001863 ** 
## nox         -17.713540   3.679308  -4.814 1.97e-06 ***
## rm            3.814394   0.408480   9.338  < 2e-16 ***
## dis          -1.478612   0.190611  -7.757 5.03e-14 ***
## rad           0.305786   0.066089   4.627 4.75e-06 ***
## tax          -0.012329   0.003755  -3.283 0.001099 ** 
## ptratio      -0.952211   0.130294  -7.308 1.10e-12 ***
## black         0.009321   0.002678   3.481 0.000544 ***
## lstat        -0.523852   0.047625 -10.999  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.74 on 493 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7343 
## F-statistic: 117.3 on 12 and 493 DF,  p-value: < 2.2e-16

Comparing models

We can use the anova() function to quantify the extent to which the quadratic fit is superior to the linear fit:

lm.fit1 <- lm(medv ~ lstat)  # simple LR model 
anova(lm.fit1, lm.fit2)      # comparing 2 models
## Analysis of Variance Table
## 
## Model 1: medv ~ lstat
## Model 2: medv ~ lstat + I(lstat^2)
##   Res.Df   RSS Df Sum of Sq     F    Pr(>F)    
## 1    504 19472                                 
## 2    503 15347  1    4125.1 135.2 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here Model 1 represents the linear submodel containing only one predictor, lstat, while Model 2 corresponds to the larger quadratic model that has two predictors, lstat and lstat^2. The anova() function performs a hypothesis test comparing the two models. The null hypothesis is that the two models fit the data equally well, and the alternative hypothesis is that the full model is superior. Here the \(F\)-statistic is \(135\) and the associated \(p\)-value is virtually zero. This provides very clear evidence that the model containing the predictors lstat and lstat^2 is superior to the model that only contains the predictor lstat. This is not surprising, since earlier we saw evidence for non-linearity in the relationship between medv and lstat.

Diagnostic plots

To show diagnostic plots for a model lm.fit2, we type:

opar <- par(no.readonly = TRUE)  # store original par
par(mfrow = c(2, 2))  # 2 rows and 2 cols

plot(lm.fit2, 
     pch = 20, col = adjustcolor("black", .25))

par(opar)  # restore original par

Interpretation:

• When the lstat^2 term is included in the model, there is little discernible pattern in the residuals.

Higher-order polynomial models

In order to create a cubic fit, we could include a predictor of the form I(X^3). However, this approach can start to get cumbersome for higher-order polynomials. A better approach involves using the poly() function to create the polynomial within lm(). For example, the following command produces a fifth-order polynomial fit:

lm.fit5 <- lm(medv ~ poly(lstat, 5))
summary(lm.fit5)
## 
## Call:
## lm(formula = medv ~ poly(lstat, 5))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.5433  -3.1039  -0.7052   2.0844  27.1153 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       22.5328     0.2318  97.197  < 2e-16 ***
## poly(lstat, 5)1 -152.4595     5.2148 -29.236  < 2e-16 ***
## poly(lstat, 5)2   64.2272     5.2148  12.316  < 2e-16 ***
## poly(lstat, 5)3  -27.0511     5.2148  -5.187 3.10e-07 ***
## poly(lstat, 5)4   25.4517     5.2148   4.881 1.42e-06 ***
## poly(lstat, 5)5  -19.2524     5.2148  -3.692 0.000247 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.215 on 500 degrees of freedom
## Multiple R-squared:  0.6817, Adjusted R-squared:  0.6785 
## F-statistic: 214.2 on 5 and 500 DF,  p-value: < 2.2e-16

This suggests that including additional polynomial terms, up to fifth order, leads to an improvement in the model fit! However, further investigation of the data reveals that no polynomial terms beyond fifth order have significant \(p\)-values in a regression fit.

By default, the poly() function orthogonalizes the predictors: this means that the features output by this function are not simply a sequence of powers of the argument. However, a linear model applied to the output of the poly() function will have the same fitted values as a linear model applied to the raw polynomials (although the coefficient estimates, standard errors, and \(p\)-values will differ). In order to obtain the raw polynomials from the poly() function, the argument raw = TRUE must be used.

Other transformations

Of course, we are in no way restricted to using polynomial transformations of the predictors. Here we try a log transformation:

summary(lm(medv ~ log(rm), data = Boston))
## 
## Call:
## lm(formula = medv ~ log(rm), data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.487  -2.875  -0.104   2.837  39.816 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -76.488      5.028  -15.21   <2e-16 ***
## log(rm)       54.055      2.739   19.73   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.915 on 504 degrees of freedom
## Multiple R-squared:  0.4358, Adjusted R-squared:  0.4347 
## F-statistic: 389.3 on 1 and 504 DF,  p-value: < 2.2e-16

(For details on linear and non-linear transformations, see Gelman & Hill, 2007, Chapter 4.)

Defining contrasts

The contrasts() function returns the coding that R uses for the dummy variables:

# attach(Carseats)
contrasts(Carseats$ShelveLoc)
##        Good Medium
## Bad       0      0
## Good      1      0
## Medium    0      1

Use ?contrasts to learn about other contrasts, and how to set them.

R has created a ShelveLocGood dummy variable that takes on a value of 1 if the shelving location is good, and 0 otherwise. It has also created a ShelveLocMedium dummy variable that equals 1 if the shelving location is medium, and 0 otherwise. A bad shelving location corresponds to a zero for each of the two dummy variables.

Interpretation: The fact that the coefficient for ShelveLocGood in the regression output is positive indicates that a good shelving location is associated with high sales (relative to a bad location). And ShelveLocMedium has a smaller positive coefficient, indicating that a medium shelving location is associated with higher sales than a bad shelving location, but lower sales than a good shelving location.