Inference for the Response Variable - Cats example

Loading the data from the MASS package and cleaning the names

Checking the data

cats |>
  slice_sample(n = 10)

##    body heart
## 1   2.1   7.6
## 2   2.9  10.1
## 3   2.5  11.0
## 4   3.6  15.0
## 5   2.9  11.8
## 6   3.2  12.3
## 7   2.2   8.7
## 8   2.5   8.8
## 9   2.2  11.0
## 10  3.1  12.1

Scatter plot for body vs heart:

We should start with a plot of body vs heart weight:

gg_cats <- 
  ggplot(
    data = cats,
    mapping = aes(
      x = body,
      y = heart
    )
  ) + 
  geom_point() + 
  theme_bw() + 
  labs(
    title = 'Cats: Body weight vs heart weight',
    subtitle = 'LOESS line added',
    x = 'Body',
    y = 'Heart'
  ) + 
  # Adding the units to the axes
  scale_x_continuous(
    labels = scales::label_number(suffix = ' kg')
  ) + 
  scale_y_continuous(
    labels = scales::label_number(suffix = ' g')
  )

gg_cats +
  geom_smooth(
    method = 'loess',
    formula = y ~ x,
    se = F
  )

Fitting the linear model:

While we’ve been calculating the slope ‘by hand’, we typically use a function to do it for us, like lm(y ~ x, data = ...) and we can get the summary stats from summary(model)

cats_lm <- lm(heart ~ body, data = cats)

summary(cats_lm)

## 
## Call:
## lm(formula = heart ~ body, data = cats)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5694 -0.9634 -0.0921  1.0426  5.1238 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.3567     0.6923  -0.515    0.607    
## body          4.0341     0.2503  16.119   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.452 on 142 degrees of freedom
## Multiple R-squared:  0.6466, Adjusted R-squared:  0.6441 
## F-statistic: 259.8 on 1 and 142 DF,  p-value: < 2.2e-16

The broom package has some useful functions that are popular as they return objects that are easier to work with

• tidy(model) returns the coefficient table as a data frame

• glance(model) returns a 1 row data frame with the fit statistics of the model

  • We’ll get to that later

• augment_columns(model, data) will add the predicted response, residuals, and other useful stats to the data frame

For now, we’ll just work with tidy()

cats_lm_table <- broom::tidy(cats_lm)

cats_lm_table

## # 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

# Finding the MSE:
sum(cats_lm$residuals^2) / (nrow(cats) - 2)

## [1] 2.109388

# Finding S_XX
sum((cats$body - mean(cats$body))^2)

## [1] 33.67972

Inference for a mean of the response

What if I want to know the average heart weight for a cat that weighs 3 kilograms?

Instead of wanting to know how the heart weight changes as body weight increases, I’m interested in the mean for a particular body weight.

From the slides, we can get a confidence interval for the mean when \(X = X_h\)

\[\hat{Y}_h \pm t(1-\alpha/2; n - 2) \sqrt{MSE\left[\frac{1}{n} + \frac{(X_h - \bar{X})^2}{S_{XX}} \right]}\]

We’ll create the confidence interval ‘by hand’ first:

Confidence interval for \(X_h = 3\) by hand

Start by calculating the different statistics we need:

1. \(\hat{Y}_h\)

2. \(t(1-\alpha/2; n - 2)\)

3. \(MSE\)

4. \(S_{XX}\)

5 \((X_h - \bar{X})^2\)

X_h <- 3
n <- nrow(cats)

# 1) Pred Y       Intercept                     slope
Y_hat_h <-   cats_lm$coefficients[1] + cats_lm$coefficients[2] * X_h


# 2) Critical value
crit_val <- qt(1 - 0.05/2, df = n - 2)

# 3) MSE
MSE_cats <- sum(cats_lm$residuals^2) / (n-2)

# 4) S_XX
S_XX <- sum((cats$body-mean(cats$body))^2)

# 5) X_h squared deviance
dev_X_h <- (X_h - mean(cats$body))^2

# The standard error using 3, 4, and 5
SE_h <- sqrt(MSE_cats * (1/n + dev_X_h/S_XX))

# Display results
tribble(
              ~ stats,                ~ value,
  'pred Y when X = 3',      round(Y_hat_h, 3),
     'Critical Value',     round(crit_val, 3),
                'MSE',     round(MSE_cats, 3),
               'S_XX',         round(S_XX, 3),
      '(X - X_bar)^2',      round(dev_X_h, 3),
     'Standard Error',         round(SE_h, 3)
)

## # A tibble: 6 × 2
##   stats              value
##   <chr>              <dbl>
## 1 pred Y when X = 3 11.7  
## 2 Critical Value     1.98 
## 3 MSE                2.11 
## 4 S_XX              33.7  
## 5 (X - X_bar)^2      0.076
## 6 Standard Error     0.139

Now that we have all the pieces, we can create the confidence interval:

pred_Y_h_CI <- 
  c('lower' = Y_hat_h - crit_val * sqrt(MSE_cats *(1/n + dev_X_h/S_XX)),
    'upper' = Y_hat_h + crit_val * sqrt(MSE_cats *(1/n + dev_X_h/S_XX))
    )

pred_Y_h_CI

## lower.(Intercept) upper.(Intercept) 
##          11.46995          12.02110

We can be 95% confident that the average heart weight of cats that weigh 3 kg is between 11.5 g to 12 g.

We can also visualize the confidence interval for the predictions using geom_smooth(method = 'lm')

gg_cats +
  geom_smooth(
    method = 'lm',
    color = 'steelblue',
    fill = 'steelblue',
    formula = y ~ x,
    level = 0.95       # Confidence level
  ) + 
  labs(subtitle = 'Linear regression line with 95% confidence band')

The shaded area is the confidence interval for the average response of \(Y\) given the value of \(X\).

Prediction intervals

The confidence intervals in the previous section tries to estimate the mean response, aka, the average of all cats that weigh 3 kg.

What if we don’t want the expected value (mean) of the response, but instead to know the range where 90%, 95%, or 99% of the population will be for a given value of \(X\)?

Instead, we want to make a prediction interval. Like with the confidence interval, we can do it ‘by hand’, it just takes a small adjustment!

From the slides, we found the variance for \(Y_h\) to be:

\[Var(Y_h) = Var(b_0 + b_1 X_h + \varepsilon_h)\]

Looks similar to the variance for \(\hat{Y}_h\), but now we also have \(\varepsilon_h\). Weirdly, \(b_0\) and \(b_1\) are independent of \(\varepsilon_h\) (but NOT each other!), we can split the variance term above:

\[Var(Y_h) = Var(b_0 + b_1 X_h) + Var(\varepsilon_h)\]

From previous results and our model assumptions, we get:

\[Var(Y_h) = \sigma^2\left[\frac{1}{n} + \frac{(X_h - \bar{X})^2}{S_{XX}} \right] + \sigma^2\]

which simplifies down to:

\[Var(Y_h) = \sigma^2\left[1 +\frac{1}{n} + \frac{(X_h - \bar{X})^2}{S_{XX}} \right]\]

The standard error is then:

\[SE(Y_h) = \sqrt{MSE\left[1 +\frac{1}{n} + \frac{(X_h - \bar{X})^2}{S_{XX}} \right]}\]

The only difference is the \(1\) inside the square root of the standard error!

Let’s update the SE for the prediction interval:

SE_h_PI <- sqrt(MSE_cats * (1 + 1/n + dev_X_h/S_XX))
round(SE_h_PI, 2)

## [1] 1.46

While we’re “only” adding 1 to go from the CI to PI, the standard error increases over 10 times!

Prediction interval for \(X_h = 3\)

c('lower' = as.numeric(Y_hat_h) - crit_val * SE_h_PI,
  'upper' = as.numeric(Y_hat_h) + crit_val * SE_h_PI) |>
  round(2)

## lower upper 
##  8.86 14.63

We are 95% confident that a cat that weights 3kg will have a heart weight between 8.9 to 14.6 grams.

Prediction band for cat heart weights

Like what we saw for the confidence interval, we can create a prediction band around the line of best fit.

Unfortunately, there isn’t a quick way to create the prediction band like what we had with geom_smooth(...). We have to ‘manually’ find the lower and upper band for the range of body weights.

We’ll create the band for the range of body weights in increments of 0.01 kg:

# Creating the lower and upper range of the prediction band:
cats_PI_band <- 
  tibble(
    # Range of body weight by 0.01 kg
    body = seq(min(cats$body), max(cats$body), by = 0.001),
    
    # Y_hat for each X:
    Y_hat = cats_lm$coef[1] + cats_lm$coef[2] * body,
    
    # Calculating (X - X_bar)^2
    body_dev2 = (body - mean(cats$body))^2,
    
    # Calculating the SE for each X
    SE = sqrt(MSE_cats * (1 + 1/n + body_dev2 / S_XX)),
    
    # Lower boundary:
    lower_95_PI = Y_hat - crit_val * SE,
    
    # upper boundary:
    upper_95_PI = Y_hat + crit_val * SE
  )

cats_PI_band |>
  round(2)

## # A tibble: 1,901 × 6
##     body Y_hat body_dev2    SE lower_95_PI upper_95_PI
##    <dbl> <dbl>     <dbl> <dbl>       <dbl>       <dbl>
##  1  2     7.71      0.52  1.47        4.81        10.6
##  2  2     7.72      0.52  1.47        4.81        10.6
##  3  2     7.72      0.52  1.47        4.82        10.6
##  4  2     7.72      0.52  1.47        4.82        10.6
##  5  2     7.73      0.52  1.47        4.82        10.6
##  6  2     7.73      0.52  1.47        4.83        10.6
##  7  2.01  7.74      0.51  1.47        4.83        10.6
##  8  2.01  7.74      0.51  1.47        4.84        10.6
##  9  2.01  7.74      0.51  1.47        4.84        10.6
## 10  2.01  7.75      0.51  1.47        4.85        10.6
## # ℹ 1,891 more rows

Now we’ll add it to the scatterplot gg_cats using geom_ribbon(...):

gg_cats +
  # Changing the subtitle
  labs(subtitle = 'Linear regression line with 95% confidence and prediction bands') + 
  # Adding the prediction band
  geom_ribbon(
    data = cats_PI_band,
    mapping = aes(
      x = body,
      ymin = lower_95_PI,
      ymax = upper_95_PI,
      y = NULL            # Since y is mapped in gg_cats, we need to remove it 
    ),
    alpha = 0.5,
    color = 'white',
    fill = 'steelblue'
  ) + 
 
  
  # Adding the prediction line and confidence band
  geom_smooth(
    method = 'lm',
    formula = y ~ x,
    color = 'steelblue',
    fill = 'steelblue'
  )