How to use this sheet

Full worked answers are in the solutions sheet.

Everything below uses one dataset, cars, which ships with R. It records the speed of a car in miles per hour and the distance it took to stop, in feet, for 50 cars.

data(cars)
str(cars)
## 'data.frame':    50 obs. of  2 variables:
##  $ speed: num  4 4 7 7 8 9 10 10 10 11 ...
##  $ dist : num  2 10 4 22 16 10 18 26 34 17 ...
head(cars, 3)
##   speed dist
## 1     4    2
## 2     4   10
## 3     7    4

We will take \(x\) to be speed and \(Y\) to be dist, and carry this one dataset from raw numbers to the estimates, the fitted line, and how far the data sit from it.

ggplot(cars, aes(speed, dist)) +
  geom_point(size = 1.4, colour = "grey25") +
  geom_smooth(method = "lm", se = FALSE, colour = "#1F5FA9", linewidth = 0.7) +
  labs(x = "Speed (mph)", y = "Stopping distance (ft)")

1. The toolkit

This section is reference material, not a lesson. Skim it, then use it in Section 2.

1.1 Sums of squares

The three quantities the whole course runs on:

x <- cars$speed
y <- cars$dist
n <- length(x)

Sxx <- sum((x - mean(x))^2)
Syy <- sum((y - mean(y))^2)
Sxy <- sum((x - mean(x)) * (y - mean(y)))

c(n = n, Sxx = Sxx, Syy = Syy, Sxy = Sxy)
##        n      Sxx      Syy      Sxy 
##    50.00  1370.00 32538.98  5387.40

R has no built-in Sxx. What it has is var and cov, which divide by \(n - 1\):

c(var_x = var(x), Sxx_over_n_minus_1 = Sxx / (n - 1))
##              var_x Sxx_over_n_minus_1 
##           27.95918           27.95918
c(cov_xy = cov(x, y), Sxy_over_n_minus_1 = Sxy / (n - 1))
##             cov_xy Sxy_over_n_minus_1 
##           109.9469           109.9469

1.2 Fitting

fit <- lm(dist ~ speed, data = cars)
coef(fit)
## (Intercept)       speed 
##  -17.579095    3.932409

summary(fit) prints a table you will meet again all term. Only the first column, Estimate, is used in this sheet: the other three are the subject of next week’s lectures, so ignore them for now.

summary(fit)$coefficients
##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) -17.579095  6.7584402 -2.601058 1.231882e-02
## speed         3.932409  0.4155128  9.463990 1.489836e-12
summary(fit)$sigma
## [1] 15.37959

The second line is \(\hat\sigma\), which Stage 3 builds from scratch.

1.3 Fitted values and residuals

yhat <- fitted(fit)
e    <- residuals(fit)
head(round(cbind(x, y, yhat, e), 3), 3)
##   x  y   yhat      e
## 1 4  2 -1.849  3.849
## 2 4 10 -1.849 11.849
## 3 7  4  9.948 -5.948

2. Stopping distances, in three stages

Stage 1: the estimates, by hand

Lecture 3 gave \(\hat\beta_1 = S_{xY}/S_{xx}\) and \(\hat\beta_0 = \bar Y - \hat\beta_1 \bar x\). Compute them from the summary statistics alone, with no call to lm:

b1 <- Sxy / Sxx
b0 <- mean(y) - b1 * mean(x)
c(beta1_hat = b1, beta0_hat = b0)
##  beta1_hat  beta0_hat 
##   3.932409 -17.579095

and check against the fit:

coef(fit)
## (Intercept)       speed 
##  -17.579095    3.932409

Lecture 3 also gave two other forms of the same estimator. Both agree:

r <- Sxy / sqrt(Sxx * Syy)
k <- (x - mean(x)) / Sxx                    # the weights k_i

c(from_ratio   = Sxy / Sxx,
  from_r       = r * sd(y) / sd(x),
  from_weights = sum(k * y))
##   from_ratio       from_r from_weights 
##     3.932409     3.932409     3.932409

Questions.

  1. \(\hat\beta_1 = S_{xY}/S_{xx}\), and also \(\hat\beta_1 = \widehat{\mathrm{Cov}}(x,Y)/\widehat{\mathrm{Var}}(x)\) using the \(n-1\) versions. Show why the two agree, in one line.
  2. The fitted slope is about \(3.93\). Say what that means about cars, with units, and be careful to make it a statement about a mean.
  3. The fitted intercept is about \(-17.6\). What does the literal reading claim, why is it absurd, and what has gone wrong? The speeds in this dataset run from 4 to 25 mph.
  4. \(\hat\beta_1 = \sum_i k_i Y_i\) with \(k_i = (x_i - \bar x)/S_{xx}\). The \(k_i\) do not involve \(Y\) at all. What does that buy us when we come to compute \(E(\hat\beta_1)\) and \(\mathrm{Var}(\hat\beta_1)\)?

Stage 2: the fitted line and its residuals

Three identities from Lecture 4. Check all three numerically:

c(sum_residuals   = sum(e),
  sum_x_times_e   = sum(x * e),
  sum_yhat        = sum(yhat),
  sum_y           = sum(y))
## sum_residuals sum_x_times_e      sum_yhat         sum_y 
##  1.110223e-14  2.682299e-13  2.149000e+03  2.149000e+03

and the fact that the line passes through the centre of the data:

c(fitted_at_xbar = b0 + b1 * mean(x), ybar = mean(y))
## fitted_at_xbar           ybar 
##          42.98          42.98

The residuals plotted against \(x\):

ggplot(data.frame(x, e), aes(x, e)) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_point(size = 1.3, colour = "grey25") +
  labs(x = "Speed (mph)", y = "Residual (ft)")

Questions.

  1. sum(e) prints as something like 1.1e-14 rather than exactly 0. Is the identity wrong? What is R telling you?
  2. \(\sum_i e_i = 0\) and \(\sum_i x_i e_i = 0\) are the two normal equations. Derive them from \(\partial Q / \partial \beta_0 = 0\) and \(\partial Q / \partial \beta_1 = 0\), where \(Q = \sum_i (Y_i - \beta_0 - \beta_1 x_i)^2\).
  3. Which of the four facts checked above are algebraic consequences of the normal equations, true for any dataset whatsoever, and which need an assumption about how the data were generated?
  4. Looking at the residual plot: the spread on the right looks wider than on the left. Which Gauss-Markov condition is that a threat to, and what would the consequence be?

Stage 3: how far the points sit from the line

SSE   <- sum(e^2)
MSE   <- SSE / (n - 2)
sigma_hat <- sqrt(MSE)

c(SSE = SSE, MSE = MSE, sigma_hat = sigma_hat)
##         SSE         MSE   sigma_hat 
## 11353.52105   236.53169    15.37959

There is a shortcut that avoids computing the residuals at all:

c(direct = SSE, shortcut = Syy - b1 * Sxy)
##   direct shortcut 
## 11353.52 11353.52

and R prints \(\hat\sigma\) under a name of its own:

summary(fit)$sigma
## [1] 15.37959

Questions.

  1. Prove the shortcut: \(\mathrm{SSE} = S_{YY} - \hat\beta_1 S_{xY}\). One line, starting from \(e_i = (Y_i - \bar Y) - \hat\beta_1 (x_i - \bar x)\).
  2. Why \(n - 2\) and not \(n\) or \(n - 1\)? Answer in terms of what was estimated before the residuals could be formed.
  3. \(\hat\sigma\) is about \(15.4\) feet. Say in one sentence what that number describes.
  4. Suppose every stopping distance in the dataset were recorded in metres instead of feet. Which of \(\hat\beta_1\), \(\hat\beta_0\), \(\hat\sigma\) and \(r\) change, and by what factor?

Full worked answers to every question above are in the companion solutions sheet, ProblemSet2-solutions.