Diagnostic Plots: Q-Q Plot

Cats Data

Like we have with the previous R examples, we’ll start by loading the packages and getting our data:

library(tidyverse); library(broom)
cats <- read.csv('https://raw.githubusercontent.com/Shammalamala/STA2410/refs/heads/main/data/Ch3/cats_with_names.csv')

slice_sample(cats, n = 10)

##       name sex body heart
## 1     Dash   M  3.2  11.6
## 2  Jasmine   F  2.3   9.0
## 3     Lily   F  2.2  11.0
## 4  Biscuit   M  2.5   9.3
## 5     Thor   M  2.9  10.6
## 6  Nibbles   M  3.5  12.9
## 7   Stripe   M  2.7   9.0
## 8    Comet   M  2.7   9.8
## 9  Charlie   M  2.4   7.9
## 10   Sassy   F  2.1   8.7

Fitting the linear model

We’ll fit our linear model using the built-in function, lm() and save the results as cats_lm.

cats_lm <- lm(heart ~ body, cats)

# View the resulting coefficients using tidy() from the broom package:
tidy(cats_lm) 

## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   -0.357     0.692    -0.515 6.07e- 1
## 2 body           4.03      0.250    16.1   6.97e-34

The p-value in the table above is only reliable IF the other pieces of our model are true. The complete model for linear regression is:

\[Y_i \sim N(\beta_0 + \beta_1 X_i, \sigma^2)\]

If any part of the model is not true, the p-value can be small, even when \(\beta_1 = 0\). So we need to check the other assumptions are correct.

Assumptions:

The assumptions for linear regression can be remembered with the acronym LINE:

• Linear association: The relationship between \(X\) and \(Y\) needs to be linear

• Independent: The rows in the data need to be independent

  • The cats can’t be related.

• Normal errors: The residuals should appear to be approximately Normal

• Equal Spread: We need constant variance so we can estimate a single variance, \(\sigma^2\)

  • Occasionally called homoscedasticity

We can check the LIE part of our acronym with a Residual Plot and the Normality assumption with a Q-Q Plot

Diagnostic Plot 2: Histogram and Q-Q Plot

Creating the histogram

If we want to look at if a variable is Normal, we can look at a histogram. The difference is that instead of looking at a histogram of \(Y\), we create a histogram for the residuals: \(e_i = y_i - \hat{y}_i\)

First, we need to get our residuals. We can do it by hand with cats\$res <- cats\$heart - cats_lm\$fit or simpler cats\$res <- cats_lm\$res. Alternatively, as seen below, we can use the augment() function from the broom package:

cats2 <- 
  augment(
    x = cats_lm,  # The linear model we'll get the residuals from
    data = cats   # The data to calculate the fitted and residuals
  )

slice_sample(cats2, n = 10)

## # A tibble: 10 × 10
##    name     sex    body heart .fitted  .resid    .hat .sigma  .cooksd .std.resid
##    <chr>    <chr> <dbl> <dbl>   <dbl>   <dbl>   <dbl>  <dbl>    <dbl>      <dbl>
##  1 Pounce   M       3.5  15.6   13.8   1.84   0.0248    1.45  2.09e-2     1.28  
##  2 Marigold F       2.3   9.7    8.92  0.778  0.0123    1.46  1.81e-3     0.539 
##  3 Dash     M       3.2  11.6   12.6  -0.952  0.0137    1.46  3.02e-3    -0.660 
##  4 Rosie    F       2.3   9      8.92  0.0783 0.0123    1.46  1.83e-5     0.0543
##  5 Tiger    M       2.2   7.6    8.52 -0.918  0.0151    1.46  3.11e-3    -0.637 
##  6 Angel    F       2.2   8.7    8.52  0.182  0.0151    1.46  1.22e-4     0.126 
##  7 Merlin   M       2.9  11.3   11.3  -0.0421 0.00787   1.46  3.36e-6    -0.0291
##  8 Lynx     M       3.4  12.8   13.4  -0.559  0.0205    1.46  1.59e-3    -0.389 
##  9 Felix    M       2.3   9.6    8.92  0.678  0.0123    1.46  1.37e-3     0.470 
## 10 Chloe    F       2     9.5    7.71  1.79   0.0225    1.45  1.78e-2     1.25

All the columns augment() added have a . before their name to make it easier to ID which columns it added.

Now we can use cats2 to create our histogram:

gg_cat_reshist <- 
  ggplot(
    data = cats2,
    mapping = aes(x = .resid)
  ) + 
  geom_histogram(
    breaks = seq(min(cats2$.resid), max(cats2$.resid)+0.2, by = 0.4),
    fill = 'steelblue',
    color = 'white'
  ) +
  theme_bw() + 
  labs(
    x = 'Heart Weight Residuals',
    title = 'Histogram of Residuals of cats heart weight by body weight'
  ) +
  scale_y_continuous(expand = c(0, 0, 0.05, 0))

gg_cat_reshist

The histogram looks okay. Has a rough bellshape to it.

Histograms can do well in identifying the shape a variable takes, but if have a specific distribution we want to check, we can use a more specific tool: the Q-Q Plot

Graph 3: Q-Q Plot

Creating a Q-Q plot by hand is complicated, but pretty straight forward in R.

If we’re using the ggplot2 framework, we can create our Q-Q plot similar to how we made the histogram with a couple of important changes:

• x = .resid becomes sample = .resid

• geom_histogram() becomes geom_qq()

• add a geom_qq_line() to add the line the points should follow

ggplot(
  data = cats2,
  mapping = aes(sample = .resid)
) + 
  # Distribution defaults to normal (qnorm), but included as a demo
  geom_qq_line(
    distribution = qnorm,
    color = 'red'
  ) + 
  geom_qq(distribution = qnorm) +

  theme_bw() + 
  labs(
    x = 'Heart Weight Residuals',
    title = 'Q-Q Plot of Residuals of cats heart weight by body weight'
  ) 

It has a bit of a wiggle to it, especially in the middle. So maybe the residuals aren’t Normal :(

While the Q-Q plot is typically better for helping assess Normality, there is still subjectivity and grey area that can happen.

So what do we do?

Look for part 3!