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

Confidence Interval for \(E(Y|X)\) formula

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)
) |> gt::gt()

stats	value
pred Y when X = 3	11.746
Critical Value	1.977
MSE	2.109
S_XX	33.680
(X - X_bar)^2	0.076
Standard Error	0.139

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

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

round(pred_Y_h_CI, 2)

## lower upper 
## 11.47 12.02

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

Calculating the confidence interval using predict()

Instead of calculating the 95% confidence interval by hand, we can use predict() with the following arguments:

1. object = the lm object we created
2. newdata = the data frame that has the new values we want to predict
   • The columns of the data frame have to be the same as the explanatory variables used to build the model
3. interval = 'confidence' will have it create a confidence interval
4. level = can be used to determine the confidence level (defaults to 0.95).

# Body weight for the CI
new_cats <- 
  data.frame(
    name = c('Donut', 'Ferdinand', 'Salem'),
    body = c(2.5, 3, 3.5)
  )

# Adding the predicted values and interval to the new_cats data
new_cats <- 
  new_cats |> 
  mutate(
    predict(
    object = cats_lm,
    newdata = new_cats,
    interval = 'confidence',
    level = 0.99
  ) |> data.frame()
  )
  
new_cats

##        name body       fit       lwr      upr
## 1     Donut  2.5  9.728494  9.380351 10.07664
## 2 Ferdinand  3.0 11.745526 11.381561 12.10949
## 3     Salem  3.5 13.762557 13.164889 14.36022

Plotting the confidence band

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!

Prediction Interval Formula

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!

Prediction Interval by hand:

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 weighs 3 kg 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 and use predict() to make the interval:

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


# We can find the interval using predict() 
# like we did with the confidence interval, but interval = 'prediction'
cats_PI_band <- 
  cats_PI_band |> 
  mutate(
    predict(
      object = cats_lm,
      newdata = cats_PI_band,
      interval = 'prediction',
      level = 0.95
    ) |> data.frame()
  )

cats_PI_band |>
  round(2)

## # A tibble: 1,901 × 4
##     body   fit   lwr   upr
##    <dbl> <dbl> <dbl> <dbl>
##  1  2     7.71  4.81  10.6
##  2  2     7.72  4.81  10.6
##  3  2     7.72  4.82  10.6
##  4  2     7.72  4.82  10.6
##  5  2     7.73  4.82  10.6
##  6  2     7.73  4.83  10.6
##  7  2.01  7.74  4.83  10.6
##  8  2.01  7.74  4.84  10.6
##  9  2.01  7.74  4.84  10.6
## 10  2.01  7.75  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 = lwr,
      ymax = upr,
      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'
  )