The fortify() function

This data frame contains the predicted values (.fitted) as well as the residuals (.resid) and standardised residuals (.stdresid). There are other bits of information there, such as the leverage (.hat) of each point, but these will not be needed at the level of these tutorials.

fort_df <- fortify(linear_model)

You can see an excerpt of the fortified data below.

head(fort_df)

Residuals vs Fitted

The residuals vs fitted plot is a diagnostic that can be used to assist checking two assumptions:

1. The zero mean assumption: the mean of the residuals must be zero at all points along this graph.
2. The homogeneity of variances (homoscedasticity) assumption: the variance (spread) of the residuals must be approximately equal across the whole plot.

The residuals vs fitted plot can be produced using the fortified data frame from the previous section.

ggplot(data = fort_df, aes(x = .fitted, y = .resid)) +
  geom_point(alpha = 0.3) +
  geom_smooth(se = FALSE) + # this makes it easier to check the zero mean assumption
  labs(x = "Fitted values", y = "Residuals") +
  theme_bw()

Note that both the x and y axis in this plot have the same units as the y variable in your linear regression model. However, we do not usually put the units on the axis labels of this graph.

The Quantile-Quantile (QQ) Plot

The normal quantile-quantile plot, also called a qq-plot, allows us to check whether the residuals from a linear regression model are normally distributed.

ggplot(data = fort_df, aes(sample = .stdresid)) +
  stat_qq() + # use stat_qq to plot the qq-plot
  geom_abline(slope = 1, intercept = 0) + # draws the one-to-one comparison line 
  theme_bw()

The residuals are considered normally distributed if the points on this plot fall along the straight line with an intercept of 0 and a slope of 1.

It is good practice to update the axis labels and make them more information, like so.

ggplot(data = fort_df, aes(sample = .stdresid)) +
  stat_qq() + # use stat_qq to plot the qq-plot
  geom_abline(slope = 1, intercept = 0) + # draws the one-to-one comparison line 
  labs(x = "q values from standard normal",
       y = "q values from standardised residuals") +
  theme_bw()

Comparing Nested Models with an \(F\)-Test: the anova() function

The anova() function can be used to perform an \(F\)-test with the null hypothesis that some linear model explains the same variability (for the population, not the data) as another linear model that is nested inside it (i.e., some parameters are set to zero).

Below, we compare the model: \[{\sf mpg} = \beta_0 + \beta_1 \cdot {\sf wt} + \beta_2\cdot{\sf drat} + \beta_3\cdot {\sf wt}\cdot {\sf drat} \] against the model \[{\sf mpg} = \beta_0 + \beta_1 \cdot {\sf wt}.\] using variables from the mtcars dataset.

model1 <- lm(data = mtcars, mpg ~ wt * drat)
model2 <- lm(data = mtcars, mpg ~ wt)
anova(model2, model1)

It is important to note that the two models must be nested for the interpretation of the \(p\)-value to match the description of the alternative hypothesis above.

An advantage of using an \(F\)-test is that it can compare models with a large difference in the number of terms. This can help avoid errors due to multiple comparisons: even if the null hypothesis is true, 5% of the time one will reject it at the \(\alpha=0.05\) level (so one wants to avoid doing too many tests if possible, recall the jellybean example cartoon.