Introduction: Fitting the Linear Model - Cats example

Loading the data from the MASS package

# Getting the data
cats <- MASS::cats

Plotting the data

Any good examination starts with exploratory data analysis. Let’s create the scatterplot between body weight and heart weight:

gg_cats <- 
  ggplot(
    data = cats,
    mapping = aes(
      x = Bwt,
      y = Hwt
    )
  ) + 
  # Adding the points to the graph:
  geom_point(color = 'steelblue') + 
  # Adding the correlation between body and heart weight
  annotate(
    geom = 'text',
    label = paste0('r = ', round(cor(cats$Bwt, cats$Hwt), 3)),
    x = 2.2,
    y = 17.5,
    size = 5,
    color = 'red'
  ) +
  # Changing the theme
  theme_bw() + 
  # Adding a title and labeling the axes
  labs(
    title = 'Cats: Body weight vs Heart weight',
    x = 'Body Weight (kg)',
    y = 'Heart Weight (g)'
  )


gg_cats

Stating the model:

Since the association between body weight (\(X\)) and heart weight (\(Y\)) is linear, we can model \(Y\) using a linear function of \(X\):

\[Y_i= \beta_0 + \beta_1 X_i + \varepsilon_i\]

\[E(Y_i|X_i) = \beta_0 + \beta_1 X_i\]

Note: For the remainder of this example, \(Y_i|X_i\) will be shortened just as \(Y_i\) unless stated otherwise

We apply the distributional assumption on the random error, \(\varepsilon_i\):

\[\varepsilon_i \sim N(0, \sigma^2)\]

which implies a distribution for \(Y_i\):

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

Fitting the model

Since the \(\beta\)’s are unknown population parameters, we need to estimate them using the sample data. So how do formulate their estimators?

\[\hat{\beta}_0 = b_0\]

\[\hat{\beta}_1 = b_1\]

The sample model is:

\[\hat{y}_i = b_0 + b_1 x_i\]

and

\[y_i = b_0 + b_1 x_i + e_i\]

where \(e_i\) is called the residual, which is a prediction of the random error, \(\varepsilon\):

\[\hat{\varepsilon}_i = e_i\]

We calculate the residual using the difference between the observed value, \(y_i\), and the estimated value, \(\hat{y}_i\):

\[e_i = y_i - \hat{y}_i\]

What is the best fitting line?

Estimating the line:

We can add the line of best fit to the scatterplot created earlier (gg_cats) by adding geom_smooth(method = 'lm'):

gg_cats_SLR <- 
  gg_cats + 
  geom_smooth(
    method = 'lm',
    formula = y ~ x,  
    se = F,           # This removes a confidence band we'll look at later
    color = 'orchid'
  )

gg_cats_SLR

The red line above uses the sample model:

\[\hat{y}_i = b_0 + b_1 x_i\]

How do we find that line? That is to say, what are the estimators for \(b_0\) and \(b_1\)?

Since we’ll be using the line of best fit for the remainder of the script, let’s find \((b_0, b_1)\) and discuss how they’re calculated.

The lm() function

Most functions in R that create a model use a formula type object to define the explanatory and response variables. The lm() function (short for linear model) is no exception!

The code chunk below will fit the SLR model for us:

cats_lm <- lm(formula = Hwt ~ Bwt, data = cats)

cats_lm

## 
## Call:
## lm(formula = Hwt ~ Bwt, data = cats)
## 
## Coefficients:
## (Intercept)          Bwt  
##     -0.3567       4.0341

The formula is coded as response ~ explanatory, which reads as ‘The response variable by the explanatory variable’, or Heart Weight by Body Weight.

The intercept is \(b_0 = -0.36\) and \(b_1 = 4.03\). Why are those the ‘best’?

Which line is the best line?

Why is the purple line better than the added red line below?

gg_cats_SLR + 
  # Adding an arbitrary line:
  geom_abline(
    slope = 3, 
    intercept = mean(cats$Hwt) - 3*mean(cats$Bwt),
    color = 'red',
    linetype = 'dashed',
    linewidth = 1
  )

Visually we can tell that the red line is better than the blue line. But how do we show that, statistically speaking?

We define some criteria that the line needs to meet to be the ‘best line’:

1. The predicted value of \(Y\) when \(X\) is at the mean is the mean of \(Y\):

\[\bar{y} = b_0 + b_1 \bar{x} \]

or, in other words, the line must pass through the coordinate \((\bar{x}, \bar{y})\):

gg_cats_SLR + 
  # Adding an arbitrary line:
  geom_abline(
    slope = 3, 
    intercept = mean(cats$Hwt) - 3*mean(cats$Bwt),
    color = 'red',
    linetype = 'dashed',
    linewidth = 1
  ) + 
  # Adding the 'center' of the plot
  geom_vline(
    mapping = aes(xintercept = mean(Bwt)),
    linewidth = 1,
    color = 'orange'
  ) + 
  geom_hline(
    mapping = aes(yintercept = mean(Hwt)),
    linewidth = 1,
    color = 'steelblue'
  )

Both the red and purple line pass through \((\bar{x}, \bar{y})\), so what other criteria does the best line need to meet?

We want the line to be as close as possible to the points in the data. One initially obvious choice is to try to minimize the average residual, the vertical distance a point is from the line:

# Creating a sample to calculate the residuals:
cats_res_sample <- 
  cats |> 
  slice(c(3, 40, 100, 130, 144)) 

# Adding the predicted heart weight:
cats_res_sample$Hwt_hat <- predict(cats_lm, newdata = cats_res_sample)

cats_res_sample

##   Sex Bwt  Hwt   Hwt_hat
## 1   F 2.0  9.5  7.711463
## 2   F 2.7  8.5 10.535307
## 3   M 3.0 10.0 11.745526
## 4   M 3.4 14.4 13.359151
## 5   M 3.9 20.5 15.376182

# Visualizing the residuals for the 5 sample cats
ggplot(
  data = cats_res_sample,
  mapping = aes(
    x = Bwt,
    y = Hwt
  )
) + 
  # Adding the 'center' of the plot
  geom_vline(
    data = cats,
    mapping = aes(xintercept = mean(Bwt)),
    linewidth = 1,
    color = 'orange'
  ) + 
  geom_hline(
    data = cats,
    mapping = aes(yintercept = mean(Hwt)),
    linewidth = 1,
    color = 'steelblue'
  ) + 
  geom_segment(
    mapping = aes(
      xend = Bwt,
      yend = Hwt_hat
    ), 
    linewidth = 1,
    color = 'orchid',
    linetype = 'dashed'
  ) + 
  geom_smooth(
    data = cats,
    se = F,
    method = 'lm',
    formula = y ~ x,
    color = 'orchid'
  ) + 
  geom_point(size = 3) 

Let’s calculate the average residual for our best fitting line. You can ‘extract’ the residuals using cats_lm$res

round(cats_lm$residuals[1:10], 3)

##      1      2      3      4      5      6      7      8      9     10 
## -0.711 -0.311  1.789 -0.915 -0.815 -0.515 -0.015  0.085  0.185  0.385

Let’s calculate the mean of the residuals, \(\bar{e}\):

mean(cats_lm$residuals)

## [1] -3.58901e-17

which is computationally 0.

How about our other line? What is the average residual of the red line:

\[\widetilde{Hwt} = 2.46 + 3 \times Bwt\]

Note: We use \(\widetilde{Hwt}\) to indicate that it is different that the predicted heart weight from the best model, \(\widehat{Hwt}\)

cats |>
  # Adding the predicted heart weight using the red line and the residual
  mutate(
    hwt_pred = mean(Hwt) - 3*mean(Bwt) + 3 * Bwt,
    hwt_res = Hwt - hwt_pred
  ) |>
  # Averaging the residuals:
  summarize(
    res_avg = mean(hwt_res)
  )

##         res_avg
## 1 -1.782109e-15

which is also computationally 0! So both lines have an average residual of 0, therefore using that criteria would indicate that both are equally good, despite what our eyes told us earlier.

In fact, any line that passes through \((\bar{x}, \bar{y})\) will have an average residual of 0. That sounds weird at first, but you can demonstrate it mathematically:

If a line passes through \((\bar{x}, \bar{y})\), then for any slope, \(b_1\),

\[b_0 = \bar{y} - b_1\bar{x}\]

which gives us the residual as:

\[e_i = y - (\bar{y} - b_1\bar{x} + b_1x_i) \]

\[e_i = (y - \bar{y}) - b_1(x_i - \bar{x})\]

Instead, we minimize the sum (or mean) of the squared residuals!

The Sums of Squares Criterion and the Normal Equations

First, we’ll define the function \(Q\):

\[Q = \sum_{i = 1}^n [Y_i - (\beta_0 - \beta_1X_i)]^2\]

We want to find the values of \((\beta_0, \beta_1)\) that minimize \(Q\)!

Let’s compare the two models we have so far:

cats |>
  # Adding the predicted heart weight using the red line and the residual
  mutate(
    hwt_hat  = cats_lm$fit,                        # Preds from the best line
    hwt_pred = mean(Hwt) - 3*mean(Bwt) + 3 * Bwt   # Preds from other line
  ) |>
  # Averaging the residuals:
  summarize(
    best_SSE = sum((Hwt - hwt_hat)^2),
    not_as_good_SSE = sum((Hwt - hwt_pred)^2)
  )

##   best_SSE not_as_good_SSE
## 1 299.5331        335.5464

Notice that the sums of the squares of the residuals is smaller!

How do we find the values that minimize \(Q\)? Do we try every combination possible and keep the one that makes the sums of squares the smallest?

While we theoretically could, there is a much more efficient way to do that: calculus!

We’ll do the math on the board, but we can take the derivative of \(Q\) with respect to the two \(\beta\)’s individually:

\[\frac{\partial Q}{\partial \beta_0} = \sum_{i = 1}^n (-2)(Y_i - \beta_0 - \beta_1 X_i) \]

\[\frac{\partial Q}{\partial \beta_1} = \sum_{i = 1}^n (-2X_i)(Y_i - \beta_0 - \beta_1 X_i) \]

To find the values of \((\beta_0, \beta_1)\) that minimize Q, we take the partial derivatives, set them equal to 0, then do a little simplification, we can create what we call The Normal Equations:

\[\sum_{i = 1}^n (Y_i - b_0 - b_1 X_i) = 0\]

\[\sum_{i = 1}^n (Y_i - b_0 - b_1 X_i)X_i = 0\]

We can solve the system of the Normal Equations to find the estimators of our model parameters:

\[b_1 = \frac{\sum_{i = 1}^n(X_i - \bar{X})(Y_i - \bar{Y})}{\sum_{i = 1}^n(X_i - \bar{X})^2}= \frac{S_{XY}}{S_{XX}}\] where \(S_{XY}\) is the estimated covariance between \(X, Y\) and \(S_{XX}\) is the estimated variance of \(X\)

and

\[b_0 = \bar{Y} - b_1\bar{X}\]

There are a couple of alternative ways we can write \(b_1\) that are easier to work with:

1. \[b_1 = \frac{(\sum_{i = 1}^n X_iY_i) - n\bar{X}\bar{Y}}{(\sum_{i = 1}^n X_i^2) - n\bar{X}^2}\]

which is often easier to use when calculating expected values and variances of the estimates, which we’ll do later!

2. \[b_1 = r\frac{s_Y}{s_X}\]

where \(r\) is the correlation between \(X, Y\), \(s_X\) is the estimated standard deviation of \(X\), and the same for \(Y\)

The fitted model:

Now that we have the functions for how to calculate \(b_0\) and \(b_1\), we can write our model as:

\[\hat{Y_i} = \bar{Y} + r\frac{s_Y}{s_X}(X_i - \bar{X})\]

Calculating the estimates ‘by hand’

Next, let’s calculate the estimates, \(b_0\) and \(b_1\) using the equations above

cats |>
  mutate(
    dev_x = Bwt - mean(Bwt),  # X - X_bar
    dev_y = Hwt - mean(Hwt)   # Y - Y_bar
  ) |>
  
  # Calculating b1 and b0
  summarize(
    # Slope using the 'raw' data
    b1_manual1  = sum(dev_x * dev_y) / sum(dev_x^2),
    b1_manual2  = (sum(Bwt*Hwt) - n() * mean(Bwt)*mean(Hwt))/(sum(Bwt^2) - n() * mean(Bwt)^2),
    b1_formula1 = cov(Bwt, Hwt) / var(Bwt),            # S_XY / S_XX
    b1_formula2 = cor(Bwt, Hwt) * sd(Hwt) / sd(Bwt),   # r (s_Y / s_X)
    
    # Calculating b0
    b0_manual = mean(Hwt) - b1_manual1 * mean(Bwt)
  )

##   b1_manual1 b1_manual2 b1_formula1 b1_formula2  b0_manual
## 1   4.034063   4.034063    4.034063    4.034063 -0.3566624

and we can compare it to the estimates from lm() using cats_lm$coef:

cats_lm$coef

## (Intercept)         Bwt 
##  -0.3566624   4.0340627

They’re the same!