Diagnostic Plots and Remedial Measures

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    Sage   F  2.3  10.1
## 2   Sassy   F  2.1   8.7
## 3    Dash   M  3.2  11.6
## 4   Ivory   F  2.6   8.7
## 5   Ebony   F  2.6  10.1
## 6  Marble   M  2.7   9.6
## 7   Vader   M  3.0  13.3
## 8  Nebula   M  2.7  12.0
## 9   Gizmo   M  2.4   9.3
## 10   Lynx   M  3.4  12.8

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 1: Residual Plot

A residual plot is just like our scatter plot, but instead of \(Y\) on the y-axis, we replace it with the residuals: \(e_i = y_i - \hat{y}_i\)

It’s easier to see an violation of our assumptions using the residual plot rather than the scatter plot

Creating the residual plot

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 Gus      M       3.7  11     14.6  -3.57  0.0353    1.43 0.114        -2.50 
##  2 Domino   M       2.6   9.4   10.1  -0.732 0.00740   1.46 0.000953     -0.506
##  3 Chloe    F       2     9.5    7.71  1.79  0.0225    1.45 0.0178        1.25 
##  4 Fuzzball M       3.5  15.7   13.8   1.94  0.0248    1.45 0.0232        1.35 
##  5 Willow   F       2.3   7.3    8.92 -1.62  0.0123    1.45 0.00784      -1.12 
##  6 Maple    F       2.3   7.9    8.92 -1.02  0.0123    1.45 0.00311      -0.708
##  7 Cocoa    F       2.7   8.5   10.5  -2.04  0.00696   1.45 0.00693      -1.41 
##  8 Vader    M       3    13.3   11.7   1.55  0.00921   1.45 0.00538       1.08 
##  9 Rusty    M       3.8  16.8   15.0   1.83  0.0413    1.45 0.0356        1.28 
## 10 Gracie   F       2.2   9.7    8.52  1.18  0.0151    1.45 0.00515       0.820

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 residual plot:

gg_cat_resplot <- 
  ggplot(
    data = cats2,
    mapping = aes(x = body, y = .resid)
  ) + 
  geom_point() +
  geom_hline(yintercept = 0, color = 'red') +
  theme_bw() + 
  labs(
    x = 'Body Weight (kg)', 
    y = 'Heart Weight Residuals',
    title = 'Residual plot of cats heart weight by body weight'
  ) 

gg_cat_resplot

Overall, it looks pretty good!

There is one thing we want to look at that is not in the LINE acronym: outliers!

There looks to be one cat that has a larger, positive residual:

gg_cat_resplot +
  geom_text(
    data = cats2 |> filter(.resid > 4), # Only getting the potential outlier
    mapping = aes(label = name),
    color = 'red',
    nudge_y = -0.3
  )

From a subjective stand point, it looks like Ranger might be an outlier. How can we be more sure?

We can convert the residuals to the standardized residuals:

\[e_i^* = \frac{e_i}{\sqrt{MSE}}\]

We can get the \(\sqrt{MSE}\) using summary(cats_lm)$sigma

cats_rMSE <- summary(cats_lm)$sigma
cats_rMSE

## [1] 1.452373

# Let's check we get the same thing:
cats2 |> 
  summarize(  # sum(e^2) / (n - 2)
    sigma2 = sum(.resid^2) / (n() - 2)
  ) |> 
  mutate(sigma = sqrt(sigma2))

## # A tibble: 1 × 2
##   sigma2 sigma
##    <dbl> <dbl>
## 1   2.11  1.45

Hooray, we get the same thing when we do it by hand!

Alternatively, you can use glance() from the broom package:

glance(cats_lm)$sigma

## [1] 1.452373

Multiple ways to ‘skin’ the cat to get the \(\sqrt{MSE}\).

Residual Plot for the standardized/studentized residuals

Let’s redo our residual plot, but now using the standardized residuals.

But why do we care about the standardized residuals? They’re equivalent to a z-score, and we can use our rule of 3 to help determine if an observation is an outlier:

\[|e_i^*| > 3 \rightarrow \text{Outlier}\]

Why 3? As long as the errors are Normal, the probability of getting a standardized residual above 3 is about 3 in 1000. So we’ll use that as our barometer if something is a residual.

Note: augment() creates a column called .std.resid, which is a type of standardized residual, but it is not the same as the one defined above. It uses something called Cook’s Distance in the calculation, and we’ll look at what Cook’s Distance is in a later chapter.

# Adding a column of the standardized residuals
cats2 <- 
  cats2 |> 
  mutate(
    stan_res = .resid / cats_rMSE
  )

# Creating the plot
gg_cats_stan_resplot <- 
  ggplot(
    data = cats2,
    mapping = aes(x = body, y = stan_res)
  ) + 
  geom_point() +
  geom_hline(yintercept = 0, color = 'red') +
  theme_bw() + 
  labs(
    x = 'Body Weight (kg)', 
    y = 'Standardized Residuals',
    title = 'Residual plot of cats heart weight by body weight'
  ) +
  # Adding tick marks at -3, -2, -1, 0, 1, 2, 3 on the y-axis
  scale_y_continuous(breaks = (-3):3,
                     minor_breaks = NULL)

gg_cats_stan_resplot

To help see if an observation falls outside of the (-3, 3) range, we can add dotted lines using geom_hline(), like we did with the line at 0:

gg_cats_stan_resplot +
  geom_hline(
    yintercept = c(-3, 3),
    color = 'red',
    linetype = 2 # makes the lines dashed instead of solid
  )

It does look like Ranger is an outlier!

So what do we do?

Let’s remove him from the data, refit the model, and see what changes:

Remedial Measure: Removing the Outlier

We’ll use filter() to remove the outlier and start over:

# Removing the outlier
cats_no_outlier <- 
  cats |> 
  filter(heart < 20)

# Next: Fit the model
cats_no_lm <- lm(heart ~ body, cats_no_outlier)

# Checking the model
tidy(cats_no_lm)

## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    0.118     0.674     0.175 8.61e- 1
## 2 body           3.85      0.244    15.7   7.62e-33

It does look like the slope changed quite a bit, from 4.04 to 3.84. The hypothesis test does still show very strong evidence of a linear association between body and heart weight.

Let’s compare the fit statistics:

bind_rows(
  .id = 'data',
  'with outlier' = glance(cats_lm), 
  'without outlier' = glance(cats_no_lm)
)

## # A tibble: 2 × 13
##   data       r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC
##   <chr>          <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl>
## 1 with outl…     0.647         0.644  1.45      260. 6.97e-34     1  -257.  520.
## 2 without o…     0.637         0.635  1.39      248. 7.62e-33     1  -249.  504.
## # ℹ 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>, nobs <int>

The fit statistics barely changed! So while Ranger is an outlier, he doesn’t seem to be too much of an influential outlier.

Let’s compare the lines in the scatter plot:

ggplot(
  data = cats,
  mapping = aes(x = body, y = heart)
) + 
  geom_point() + 
  geom_abline(
    intercept = cats_lm$coef[1], 
    slope = cats_lm$coef[2],
    color = 'blue'
  ) + 
  geom_abline(
    intercept = cats_no_lm$coef[1], 
    slope = cats_no_lm$coef[2],
    color = 'red'
  )

Visually speaking, the two lines are very, very similar.

If we check the resulting residual plot:

# Calculating the standardized residuals for the no outlier data set
cats_no_outlier <- 
  cats_no_outlier |> 
  mutate(
    pred_heart = cats_no_lm$fit,
    res_heart = heart - pred_heart,
    stan_res = res_heart / summary(cats_no_lm)$sigma
  )
  
# Creating the residual plot for the cats without Ranger
gg_cats_stan_resplot +
  geom_hline(
    yintercept = c(-3, 3),
    linetype = 2,
    color = 'red'
  ) + 
  labs(title = 'Residual Plot without Ranger: the cat with the large heart') +
  # Replacing the data set with the one without any outliers
  cats_no_outlier

The residual plot looks great! No clear trend (Linear), there’s no funnel pattern (Equal Spread).