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? 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 is 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) \]
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.
fit <- lm(amnt_earnings_med_10y ~ cost_avg, data = college_scorecard)
augment(fit) |>
ggplot(aes(x = .fitted, y = .resid)) +
geom_point(size = 0.7, alpha = 0.5) +
geom_smooth()
augment(fit) |>
ggplot(aes(x = cost_avg, y = .resid)) +
geom_point(size = 0.7, alpha = 0.5)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, we say that the errors are called homoskedastic. Meanwhile, the errors are called heteroskedastic if they have non-constant variance.[[need to add influence on inference and/or predictions]]
Like for 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 variance of the residuals to be approximately constant across the fitted values. [[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, [[weighted regression]]. 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
Assumption (4) requires that the errors are normally distributed. If this assumption holds, then the sampling distribution of our parameter estimates will be normally distributed; so, we can use standard procedures for obtain valid confidence intervals and \(p\)-values [[need to explain more what these are]]. On the other hand, if 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. However, an important caveat to this is that if the number of observations is sufficiently large [[add heuristic of what this is]], the sampling distribution of our parameter estimates will be approximately normally distribution even if the errors are not normally distributed [[footnote to CLT]]. Therefore, assumption (4) is “optional” 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 our sample.
We can assess whether the errors are normally distributed by plotting a histogram of the residuals. Figure 4.3 shows an example of the histogram of the residuals, where the error terms are normally distributed. We can see that the residuals appear to be roughly normally distributed.
We can formally test whether the residuals follow a normal distribution using the Shapiro-Wilk test. The null hypothesis (\(H_0\)) of the Shapiro-Wilk test is the data is normally distributed, and the alternative hypothesis (\(H_1\)) is that the data are not normally distributed. If we obtain a \(p\)-value less than our pre-specified cut-off (e.g., \(p < 0.05\)), then we reject the null hypothesis that the data are normally distributed. If the \(p\)-value is greater than our pre-specified cut-off we fail to reject the null hypothesis that the data are normally distributed.
y.hat <- predict(fit)
e <- df$y - y.hat
shapiro.test(e)$p.value[1] 0.9099408[[shapiro wilk test example in R]]
Another option to check for non-normal errors is to use a QQ-plot (i.e., a Quantile-Quantile plot). A QQ-plot uses the model’s standardized residuals, which are the residuals divided by their standard deviation (i.e., square root of the variance). Recall by assumption:
\[ e_i \sim \]
\[ e_{\text{standardized}, i} = \frac{e_i}{\sigma} \]
There are called standardized residuals because they have mean zero and standard deviation one; that is:
\[ e_{\text{standardized}, i} \sim \mathcal{N}(0, 1) \]
where we plot the observed quantiles of the standardized residuals against the theoretical quantiles of a normal distribution with mean zero and variance one [[need to better introduce standardized residuals]]. If the errors are normally distributed, we expect observed quantiles of standardized residuals to approximately follow the theoretical quantiles of a normal distribution with mean zero and variance one. Figure 4.4 shows the QQ-plot, in which the residuals are simulated to be normally distributed. While the observed standardized residuals quantiles appear to deviate a bit from the theoretical quantiles of the normal distribution (particularly in the upper and lower quantiles), these deviations are relatively small. This indicates that the assumption of the errors being normally distributed is satisfied.
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.
| 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\):1
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() Recall that the standard error of an estimate is the square-root of its variance.↩︎