4 Linear Models

4.1 Setup

Simple linear regression describes the linear relationship between the outcome \(Y\) and a covariate of interest \(x\) (for simplicity, we will assume that the \(x\)’s are fixed and the \(Y\)’s are random variables). Under simple linear regression, we assume that for all units \(i = 1, \dots, n\):

\[ \begin{aligned} Y_i = \beta_0 + \beta_1 x_i + \epsilon_i \quad \text{with }\epsilon_i \sim \text{Normal}(0, \sigma^2) \end{aligned} \tag{4.1}\]

We can extend the simple linear regression model given in Equation 4.1 in various ways. For example, we might want to include different powers of covariate \(x\) (e.g., \(x, x^2, x^3, \dots\)) in the linear model. This is known as polynomial regression, which we discuss in more detail in Chapter 6.

Another natural extension, which we focus on in this chapter, is to include multiple covariates \(\mathbf{x} = (x_1, \dots, x_p)\), rather than just a single covariate, in the linear model. This is commonly referred to as multiple linear regression, which assumes that for all units \(i = 1, \dots, n\):

\[ Y_i = \beta_0 + \beta_1 x_{1i} + \cdots + \beta_p x_{pi} + \epsilon_i \quad \text{with } \epsilon_i \sim \text{Normal}(0, \sigma^2) \tag{4.2}\]

Notice that we can view simple linear regression as a special case of multiple linear regression, where the number of covariates \(p = 1\). For this reason, multiple linear regression is often just referred to as linear regression, and we use the term linear regression throughout this chapter.

4.2 Model Assumptions

Per Equation 4.2, linear regression makes the following assumptions:

The error terms have mean zero: \(E[\epsilon_i] = 0\), or equivalently, the mean of \(Y_i\) is a linear function of \(\mathbf{x}_i\): \[ E[Y_i] = \beta_0 + \beta_1 x_{1i} + \dots + \beta_p x_{pi}. \]
The error terms have constant variance: \(\text{Var}[\epsilon_i] = \sigma^2\).
The error terms are independent of each other. That is, having information about an error term \(\epsilon_i\) does not tell us something about another error term \(\epsilon_j\).
The error terms are normally distributed: \(\epsilon_i \sim \text{Normal}(0, \sigma^2)\)

Important: we do not directly know the linear regression parameters: \(\beta_0, \beta_1, \dots, \beta_p\), and \(\sigma^2\); instead, we estimate them using a (random) sample from the population:

\[ (\mathbf{x}_1, y_1), (\mathbf{x}_2, y_2), \dots, (\mathbf{x}_n, y_n), \]

where \(\mathbf{x}_i = (x_{i1}, x_{i2}, \dots, x_{in})\) represents the observed covariates for unit \(i\), and \(y_i\) denotes the observed outcome for unit \(i\).

To distinguish the true linear regression parameters from our estimates, we use \(b_0, b_1, \dots, b_p\) to denote our parameter estimates. As a result, our fitted linear model is:

\[ \hat{y}_i = b_0 + b_1 x_{1i} + \dots + b_p x_{pi}, \tag{4.3}\]

where \(\hat{y}_i\) is the predicted outcome, or fitted value for unit \(i\).

4.3 Model Diagnostics

Model diagnostics are necessary to assess whether the linear regression assumptions in Section 4.2 are satisfied and, consequently, whether linear regression is a good fit for our data.

An important tool for model diagnostics are the model residuals. Residuals are the observed error terms, which are the sample version of the error terms in Equation 4.2. For unit \(i\), the residual is defined as:

\[ \begin{aligned} e_i &= y_i - \hat{y}_i \\ &= y_i - (b_0 + b_1 x_{1i} + \dots + b_p x_{pi}) \end{aligned} \tag{4.4}\]

4.3.1 Mean Model Mispecification

Assumption (1) requires that the mean model of the outcome is linear with respect to the included covariates:

\[ E[Y_i] = \beta_0 + \beta_1 x_{1i} + \dots + \beta_p x_{pi}. \]

If this assumption is violated, then the linear model is not a good representation of the data-generating process, and the mean model of the outcome is misspecified.

There are two primary ways in which the mean model of the outcome may be misspecified in linear regression. First, some covariates may have a nonlinear relationship with the outcome that is not properly captured by the specified linear model. Second, the linear model may omit one or more important covariates that are associated with the outcome.

Nonlinearity: We can assess whether there is nonlinear relationship between an included covariate and the outcome by plotting the residuals against the fitted values. If the mean model is correctly specified, the residuals should be randomly scattered around zero for every fitted value, with no clear pattern.

To illustrate the difference between the fitted vs. residual consider the following toy example, where the true relationship between our outcome \(Y\) and covariate \(x\) is quadratic; that is:

\[ Y_i = 2 + 3 x_{i} + x^2_{i} + \epsilon_i, \quad \text{with } \epsilon_i \sim \text{Normal}(0, \sigma^2) \]

suppressMessages(library(tidyverse))
set.seed(4747)
error <- rnorm(n = 100, mean = 0, sd = 2)
x <- runif(n = 100, min = 1, max = 10)

df1 <- tibble(
  x = x,
  y = 2 + 3*x + x**2 + error
  ) # save this to data folder?

However, let’s suppose we fit a linear model with only \(x_i\), excluding the squared-term \(x_i^2\). Then, our linear model would be misspecified, and there would be a clear (quadratic) trend in the fitted vs. residual plot.

library(broom)
misspecified_fit <- lm(y ~ x, data = df1)

augment(misspecified_fit) |>
  ggplot(aes(x = .fitted, y = .resid)) +
  geom_point(size = 0.7) +
  geom_hline(yintercept = 0, color = 'grey20', linetype = "dashed") +
  geom_smooth() +
  labs(x = "Fitted", y = "Residual") +
  theme_bw()

In contrast, when our linear model is correctly specified, with both \(x_i\) and \(x_i^2\) included, the residuals are randomly scattered around zero for each fitted value.

Missing covariates: Our linear model may also be misspecified because we are missing are one or more important covariates. To diagnose this, we can plot the values of any potential covariates not included in the linear model against the residuals to determine if that predictor should be added.

For instance, suppose we have two covariates \(\mathbf{x}_1\) and \(\mathbf{x}_2\), but only include \(\mathbf{x}_1\) in our linear model.

\[ Y_i = 2 + 3 x_{1i} + 4x_{2i} + \epsilon_i, \quad \text{with } \epsilon_i \sim \text{Normal}(0, \sigma^2) \]

Figure 4.1 shows that the points in the plot of \(\mathbf{x}_2\) against the residuals have a linear trend. This suggests that our mean model is incorrectly specified; in particular, \(\mathbf{x}_2\) should be included in our linear model.

4.3.2 Non-Constant Error Variance

Assumption (2) of linear regression is that the error terms have constant variance. If the error variance is constant across the fitted values, the errors are called homoskedastic. Meanwhile, the errors are called heteroskedastic if they have non-constant variance across the fitted values.

Like with mean model misspecification, we can use the fitted values versus residuals plot to determine whether the constant error variance assumption is satisfied. Under constant error variance, we expect for the residuals to have approximately the same variance for each given fitted value. [[explain figure]]

There are various ways to correct for error heteroskedasticity. One approach is to transform the outcome to stablize the variance.

[[add how this would be done and what are some covariates of this]]. Alternatively, we could use a weighted regression if we can estimate the error variance \(\sigma^2_i\) for each unit \(i = 1, \dots, n\). A third approach is to use the sandwich estimator to estimate the variance of our parameter estimates. We show how to compute the sandwich estimator in R in [[ADD SECTION]].

4.3.3 Non-Normality of Errors

The fourth assumption is that the errors are normally distributed. If normality of errors holds, then the sampling distribution of our parameter estimates will be normally distributed. In turn, we can use standard procedures for obtain valid confidence intervals and \(p\)-values. On the other hand, if only the Assumptions (1)-(3) hold but Assumption (4) does not, the linear model is still perfectly valid; we just cannot trust the estimated uncertainties on the coefficients and associated \(p\)-values. In other words inference is affected, but not prediction when Assumption (4) is violated.

However, an important caveat to this is that if the number of observations is sufficiently large (a common heuristic is \(n > 30\) observations), the sampling distribution of our parameter estimates will be approximately normally distributed even if the errors are not perfectly normally distributed.¹ Therefore, a violation of Assumption (4) is not concerning if (i) we have a small number of observations but only care about prediction, not inference, or (ii) we have a sufficiently large number of observations in the sample such that the sampling distributions of our parameter parameter will be approximately normal anyway.

We can assess whether the errors are normally distributed by plotting a histogram of the residuals. For example,Figure 4.3 shows a histogram of the residuals, in a linear model where the error terms are normally distributed.

From the histogram, we can see that the residuals appear to be roughly normally distributed. Another way to visually inspect if the normality of error assumption holds is via a Quantile-Quantile (QQ) plot. In a QQ-plot, we compare the quantiles of the standard normal distribution against the quantiles of the model’s standardized residuals (i.e., the residuals divided by square root of the variance): \[ r_i = \frac{e_i}{\sqrt{\sigma^2}} \]

If the errors are normally distributed, then \(r_i \approx \mathcal{N}(0, 1)\). Hence, the quantiles of the standard normal distribution should match the quantiles of the standardized residuals. For example, Figure 4.4 shows the QQ-plot, in which the errors are simulated to follow a normal distribution. While the observed standardized residuals quantiles deviate a bit from the theoretical quantiles of the normal distribution (particularly in the upper and lower quantiles), these deviations are relatively small. This suggests that the assumption of the errors being normally distributed is satisfied.

Alternatively, we can also perform a Shapiro-Wilk test to assess whether the residuals follow a normal distribution. For the Shapiro-Wilk test, the null hypothesis (\(H_0\)) of the is that the data is normally distributed, and the alternative hypothesis (\(H_1\)) is that the data is not normally distributed. Thus, if we obtain a \(p\)-value less than our pre-specified cut-off (e.g., \(p < 0.05\)), then we would reject the null hypothesis that the data are normally distributed. In contrast, we fail to reject the null hypothesis that the data are normally distributed if the \(p\)-value is greater than our pre-specified cut-off. So, in the example below, we fail to reject the null hypothesis that the errors are normally distributed at an \(\alpha\)-level of 0.05.

set.seed(0913)
normal_errors <- rnorm(n = 100, mean = 0, sd = 1)
shapiro.test(normal_errors)$p.value

[1] 0.2902862

4.4 Outliers

4.5 Multiollinearity

4.5.1 Variance Inflation Factor

The variance inflation factor, or vif, is the factor by which the estimated variance for a coefficient is inflated due to multicollinearity. If a vif value is high, then the coefficient estimate for that variable obtained from the lm() function can be much more uncertain than R indicates.

Example: Suppose we have the following fitted model:

\[ \hat{Y} = 5 + 4 x_1 - 2 x_2, \]

then, \(b_0 = 5, b_1 = 4\), and \(b_2 = -2\). Let the estimated standard errors for \(b_1\) and \(b_2\) both be 2, and the vifs for \(b_1\) and \(b_2\) be \(4\) and \(9\), respectively. This information is summarized in Table 4.1.

Table 4.1: My Caption

Coefficient estimate	Standard Error estimate	vif
4	2	4
-2	2	9

Since the vif value is the factor by which the estimated variance is inflated, multiplying the estimated standard error by the square-root of the vif value gives us the actual standard error estimates for \(b_1\) and \(b_2\):²

The actual standard error estimate for \(b_1\): \(2 \times \sqrt{4} = 4\)
The actual standard error estimate for \(b_2\): \(2 \times \sqrt{9} = 6\).

A typical rule-of-thumb is to worry if vif values are above 5 (more conservative) or 10 (less conservative). One mitigation strategy is to remove one variable at a time until the vif values for the remaining covariates fall below 5 or 10. However, this can lead to [[post-selection problems?]]. Moreover, removing covariates can adversely affect the linear model’s predictive abilities. Thus, common alternative approaches include performing ridge regression or principal components regression. [[I think it would be best to include more details of these in other chapters]].

5 Example in R

For example, suppose we are interested in understanding the association of [[ADD]] with [[add]] and [[add]] for American colleges and universities, and fit the following linear model:

\[ \begin{aligned} \text{median earnings after graduation}_i = \beta_0 &+ \beta_1 \cdot (\text{undergraduate enrollment})_i \\ &+ \beta_2 \cdot (\text{average net cost})_i + \epsilon_i \end{aligned} \]

suppressMessages(library(tidyverse))
suppressMessages(library(collegeScorecard))

college_scorecard <- school |>
  left_join(scorecard, by = join_by(id)) |>
  filter(academic_year == "2020-21",
  deg_predominant == "Bachelor") |>
  select(name, amnt_earnings_med_10y, n_undergrads, cost_avg) |>
  drop_na()

For more details, see the Central Limit Theorem.↩︎
Recall that the standard error of an estimate is the square-root of its variance.↩︎