Remedial Measures - Fixing Assumption Violations

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  Buttercup   F  2.3   9.6
## 2     Apollo   M  2.8  10.2
## 3       Finn   M  3.1  11.5
## 4      Tabby   M  2.7   9.6
## 5   Fuzzball   M  3.5  15.7
## 6      Maple   F  2.3   7.9
## 7      Chloe   F  2.0   9.5
## 8      Titan   M  2.8  13.5
## 9      Turbo   M  3.2  11.9
## 10     Groot   M  3.1  12.1

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

Q-Q Plot for the residuals

We left off last section looking at the Q-Q plot and determined that it’s not clear that the residuals are Normal:

# Adding the residuals
cats2 <- 
  augment(
    x = cats_lm,
    data = cats
  )

# Q-Q plot
gg_qq_rescat <- 
  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'
  ) 

gg_qq_rescat

We could remove the outlier and try again:

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

# Refitting the model
cats_no_lm <- lm(heart ~ body, cats_no_outlier)

# Residuals for the no outlier data
cats_no_outlier <- 
  cats_no_outlier |> 
  augment(x = cats_no_lm)

# QQ Plot
gg_qq_rescat + cats_no_outlier

Still has the wiggle, even without the outlier :(

So what can we do?

Remedial Measure: Transforming \(Y\)

If we find that there is non-constant spread in the residual plot (\(Var(Y) \ne \sigma^2\)) or if the residuals are not Normal, we can try transforming the response variable.

There are many, many, many transformations we can apply to \(Y\), we’ll look at a subset called Power Transformations:

\[\tilde{Y}_i = Y_i^{\lambda}\]

The one exception is if \(\lambda = 0\) because \(Y_i^0 = 1\), regardless of \(Y_i\).

If \(\lambda = 0\), then \(\tilde{Y}_i = \log(Y_i)\)

Putting the pieces together, we get:

\[ \tilde{Y}_i = \begin{cases} Y_i^{\lambda} & \text{if } \lambda \ne 0 \\ \log(Y_i) & \text{if } \lambda = 0 \\ \end{cases} \]

This transformation above was created by George Box and David Cox and is called a Box-Cox transformation.

The \(\lambda\) value that minimizes a certain function is the value of \(\lambda\) that we use, sort of.

Instead of using the exact best value of \(\lambda\), we ‘round’ it to the nearest power that has an interpretation to it.

IE, if \(\lambda = 0.4326\), we’ll use \(\lambda = 0.5\)

If \(\lambda = -0.12364\), we’ll use \(\lambda = 0\)

Note: \(\log\) is the natural log, not \(\log_{10}\). In statistics, we default to the natural log unless specified otherwise!

Box-Cox Transformation

We can quickly apply a Box-Cox transformation using the boxcox(model) function from the MASS package.

I don’t load the MASS package because it also has a select() and filter() function that overrides those same functions from the dplyr package. Instead, I use MASS::func().

cats_bc <- MASS::boxcox(cats_lm)

The graph above shows the best choice of \(\lambda\) and a 95% confidence interval. We want to use a ‘good’ choice of \(\lambda\) from the choices in the interval.

It looks like both 0 and 0.5 are in the interval, and since the middle line is closest to 0, we’ll use it. If we find that a log transformation doesn’t work, we could try a square root transformation.

cats <- 
  cats |> 
  mutate(heart_log = log(heart))

tibble(cats)

## # A tibble: 144 × 5
##    name    sex    body heart heart_log
##    <chr>   <chr> <dbl> <dbl>     <dbl>
##  1 Luna    F       2     7        1.95
##  2 Bella   F       2     7.4      2.00
##  3 Chloe   F       2     9.5      2.25
##  4 Nala    F       2.1   7.2      1.97
##  5 Ginger  F       2.1   7.3      1.99
##  6 Cleo    F       2.1   7.6      2.03
##  7 Mittens F       2.1   8.1      2.09
##  8 Lucy    F       2.1   8.2      2.10
##  9 Daisy   F       2.1   8.3      2.12
## 10 Zoe     F       2.1   8.5      2.14
## # ℹ 134 more rows

Note: All the heart weights are positive (duh). BUT if \(Y\) has negative values, we need to add a constant \(C\) to \(Y\) so that \(\min(Y + C)>0\) otherwise most transformations won’t work.

Checking the transformed heart weight

Now let’s redo the analysis with the log of heart weight:

cats_log_lm <- lm(heart_log ~ body, cats)

# Creating the QQ-Plot
gg_qq_rescat + 
  labs(title = "Q-Q plot of the residuals", subtitle = "Log Transformation of Y") + 
  augment(
    cats_log_lm,
    cats
  )

It looks better, but still has a little of a wiggle in the middle and the right tail deviates from the line.

Redo the transformation part, but use \(\sqrt{Y}\) instead!