About this document

These are the answers to every question in Problem Set 4, written out in full. Read it alongside the problem set rather than instead of it, and attempt each question before you look here.

Exercise 1. Reading residual plots

(a) What each panel shows

  • A, a funnel. The residuals spread out as the fitted value grows: equal variance (E) is in doubt.
  • B, a curve. Positive residuals at both ends, negative in the middle: the mean of \(Y\) is not a straight line in \(x\), so linearity (L) fails.
  • C, no pattern. An even, structureless band around zero: nothing is in doubt. This is what a residual plot should look like.
  • D, one point far above the rest. An outlier. It may signal heavy-tailed errors (N), or simply a recording error: check the data before anything else.

(b) What is at risk

  • A (E): the estimates are still unbiased, but the standard errors, intervals and tests are wrong, and prediction intervals are too wide where the spread is small and too narrow where it is large.
  • B (L): the estimates themselves are misleading. Here the fitted slope is 0.078 and the correlation only 0.07, even though \(Y\) depends strongly on \(x\): \(r\) measures only linear association, and a straight line is the wrong summary, whatever its standard errors.
  • D (outlier): a single point inflates \(\hat\sigma\), which widens every interval, and can pull the fitted line towards it.

(c) Checking independence

Plot the residuals against time, or the order in which the observations were collected. Independent errors scatter without pattern; dependent errors show runs: long stretches of consecutive residuals with the same sign, or a slow wave. A plot of each residual against the previous one (\(e_i\) against \(e_{i-1}\)) makes it sharper: independence gives a round cloud, positive dependence a rising band.

Here is what this gives for the COVID-19 data, with the fit of Lecture 7. The rows of COVID-19.csv are not stored in date order, so the days must be sorted first; one row has no date and is left out of these plots. A few days are missing from the file, so three of the “consecutive” pairs are two or three days apart, which does not change the picture.

covid <- read.csv("COVID-19.csv")
fitC <- lm(Cases ~ Tests, data = covid)
rc <- data.frame(date = as.Date(covid$Date, format = "%m/%d/%Y"), e = resid(fitC))
rc <- rc[!is.na(rc$date), ]
rc <- rc[order(rc$date), ]                      # consecutive days, in order
lagged <- data.frame(previous = head(rc$e, -1), current = tail(rc$e, -1))

p_time <- ggplot(rc, aes(date, e)) +
  geom_hline(yintercept = 0, colour = "grey60") + geom_line(colour = "grey70") +
  geom_point(size = 1) + labs(x = "Date (2020)", y = "Residual")
p_lag <- ggplot(lagged, aes(previous, current)) +
  geom_hline(yintercept = 0, colour = "grey60") + geom_vline(xintercept = 0, colour = "grey60") +
  geom_point(size = 1) +
  labs(x = expression(paste("Previous day's residual, ", e[i - 1])),
       y = expression(paste("Residual, ", e[i])))
gridExtra::grid.arrange(p_time, p_lag, nrow = 1)

cor(lagged$previous, lagged$current)
## [1] 0.8754165

Against time (left), the residuals move in long waves: they stay negative for weeks, then positive for weeks. Against the previous day’s residual (right), the points lie along a rising band rather than in a round cloud: a day whose residual is positive is usually followed by another positive one. The correlation between the residuals of consecutive days is 0.88, far from the 0 that independent errors would give. Neither pattern can be seen when the residuals are plotted against the fitted values or against the number of tests.

Exercise 2. A very high \(R^2\) for the wrong model

fitL <- lm(height ~ age, data = Loblolly)
summary(fitL)$r.squared
## [1] 0.9799281
coef(fitL)
## (Intercept)         age 
##   -1.312396    2.590523

(a) \(R^2\)

\(R^2 = 0.9799\): age accounts for 98.0% of the variation in height. On its own, a value this close to 1 would usually be read as an excellent fit. The fitted slope says that the trees gain about 2.59 feet per year on average.

(b) The residual plot

rL <- data.frame(age = Loblolly$age, fitted = fitted(fitL), resid = resid(fitL))
ggplot(rL, aes(fitted, resid)) +
  geom_hline(yintercept = 0, colour = "grey60") + geom_point(size = 1) +
  labs(x = "Fitted height (ft)", y = "Residual (ft)")

round(tapply(rL$resid, rL$age, mean), 2)      # mean residual at each age
##     3     5    10    15    20    25 
## -2.22 -1.44  2.85  3.00  0.97 -3.16

With six ages, the points form six vertical stripes. Their centres trace an inverted U: the residuals are negative on average at the youngest and the oldest ages (-2.2 ft at 3 years, -3.2 ft at 25 years) and positive in between (3.0 ft at 15 years). The mean height is not a straight-line function of age: linearity fails. The stripes also widen from left to right, from a spread of about 0.4 ft at 3 years to 2.3 ft at 25 years: older, taller trees vary more, so equal variance is in doubt as well.

(c) Lowess against the band

ggplot(Loblolly, aes(age, height)) +
  geom_point(size = 1) +
  geom_smooth(method = "lm") +
  geom_smooth(method = "loess", colour = "orange", fill = "orange") +
  labs(x = "Age (years)", y = "Height (ft)")

Both bands are there, but they are only 1 to 2 feet wide on an axis that runs from 0 to 64 feet, so they are hard to see. To see them, subtract the fitted line from everything: the line becomes the horizontal axis, with its band around zero, and the lowess curve becomes its distance from the line, with its own band.

grid <- data.frame(age = seq(3, 25, by = 0.25))
pl <- predict(fitL, grid, se.fit = TRUE)
po <- predict(loess(height ~ age, data = Loblolly), grid, se = TRUE)
gap <- data.frame(age = grid$age,
                  line_lwr = -qt(0.975, pl$df) * pl$se.fit,
                  line_upr =  qt(0.975, pl$df) * pl$se.fit,
                  lo       = po$fit - pl$fit,
                  lo_lwr   = po$fit - pl$fit - qt(0.975, po$df) * po$se.fit,
                  lo_upr   = po$fit - pl$fit + qt(0.975, po$df) * po$se.fit)
ggplot(gap, aes(age)) +
  geom_ribbon(aes(ymin = line_lwr, ymax = line_upr), fill = "#3366FF", alpha = 0.3) +
  geom_hline(yintercept = 0, colour = "#3366FF") +
  geom_ribbon(aes(ymin = lo_lwr, ymax = lo_upr), fill = "orange", alpha = 0.4) +
  geom_line(aes(y = lo), colour = "orange", linewidth = 1) +
  labs(x = "Age (years)", y = "Distance from the line (ft)")

The lowess curve bends downward: it lies below the blue band of the line at the youngest and the oldest ages, and above it from about 10 to 20 years. The blue band shows where the true line could plausibly be if the model were right. Lowess estimates the mean height at each age without assuming a line. Leaving the band in this systematic way means that no straight line compatible with the data follows the means: the line overstates the height of young and old trees and understates it in between.

The orange band does the same job for the lowess curve: it shows where the mean height at each age could plausibly be, with no assumption on the shape. It is narrow too, and the straight line lies outside it at 6 of the 6 ages at which trees were measured. In the second plot, the line is the zero level: the orange band stays clear of zero except where the two curves cross, near 6 and 21 years. Away from those crossings, each curve lies outside the other’s band: the gap between the line and the lowess curve is much larger than the uncertainty of either.

(d) Why \(R^2\) did not warn us

\(R^2 = SSR/SST\) measures how much of the variation in \(Y\) the fitted line accounts for. Heights range from about 3 to 64 feet, and most of that spread is explained by age, so a straight line leaves only a little unexplained variation, even though it has the wrong shape. \(R^2\) does not measure whether the line is the right shape: a curved relation can be summarized well, but not correctly, by a line. The pattern in the residuals is small compared with the range of heights, and only a residual plot shows that. So, \(R^2\) is not a complete diagnostic.

(e) The same seed sources at every age

Independence. Trees from a fast-growing seed source tend to be tall at every age, so the errors of the six measurements of one source are related: knowing that a source is above the line at 10 years tells you it is probably above the line at 15 and 20 years too.

The plot in (b) cannot show this: it does not say which point belongs to which seed source. To see it, you would have to follow each source from one age to the next, for example by joining its six points with a line. With 84 heights from only 14 sources, the data contain less independent information than 84 independent trees would, and the usual standard errors are too small, as in the repeated-measurements part of Problem Set 3, Exercise 5.

Exercise 3. Stopping distances

fitC <- lm(dist ~ speed, data = cars)
rC <- data.frame(fitted = fitted(fitC), resid = resid(fitC))
p_res <- ggplot(rC, aes(fitted, resid)) +
  geom_hline(yintercept = 0, colour = "grey60") + geom_point() +
  labs(x = "Fitted distance (ft)", y = "Residual (ft)")
p_qq <- ggplot(rC, aes(sample = resid)) + stat_qq() + stat_qq_line() +
  labs(x = "Normal quantile", y = "Residual (ft)")
gridExtra::grid.arrange(p_res, p_qq, nrow = 1)

(a) The plots

The residuals spread out as the fitted distance grows, a funnel: equal variance is in doubt. There is also a slight bend (residuals mostly positive at the two ends), a hint that the mean is not quite linear. On the Q-Q plot, the largest residuals lie above the line: the upper tail is longer than a normal one (i.e., a right skew).

(b) Shapiro-Wilk

\(H_0\): the errors are normally distributed, against \(H_a\): they are not. Since \(p = 0.0215 < 0.05\), we reject \(H_0\) at the 5% level: the residuals show evidence against normality. This is the right skew seen at the top of the Q-Q plot, where a few cars took much longer to stop than the line predicts.

(c) Brown-Forsythe

The pooled standard deviation is \[s = \sqrt{\frac{22 \times 9.086^2 + 26 \times 11.334^2}{48}} = 10.364,\] and \[t_{BF} = \frac{\bar d_1 - \bar d_2}{s\sqrt{1/n_1 + 1/n_2}} = \frac{8.947 - 13.016}{10.364 \times \sqrt{1/23 + 1/27}} = -1.384.\] Since \(|t_{BF}| = 1.384 < 2.011\) (\(p = 0.173\)), we do not reject constant variance at the 5% level.

That does not mean the variance is constant. The deviations are larger in the faster group, \(\bar d_2 > \bar d_1\), in the direction of the funnel; the test simply lacks the power to confirm it with 50 cars split into two groups. Failing to reject is not evidence of equality, and the residual plot remains the main tool: here it shows the funnel clearly.

(d) Box-Cox

\(\hat\lambda = 0.43\) suggests a square root, \(\lambda = 1/2\), which lies in the 95% interval \((0.23, 0.66)\). The interval excludes \(\lambda = 1\), so the data do not support leaving dist unchanged, and it excludes \(\lambda = 0\), so a log transformation is not supported either.

(e) After the transformation

fitS <- lm(sqrt(dist) ~ speed, data = cars)
rS <- data.frame(fitted = fitted(fitS), resid = resid(fitS))
gridExtra::grid.arrange(p_res %+% rS + labs(x = "Fitted sqrt(distance)", y = "Residual"),
                        p_qq %+% rS + labs(y = "Residual"), nrow = 1)

Yes. The funnel is gone: the residuals form an even band, and the Q-Q plot is close to its line. Both tests agree: Shapiro-Wilk no longer rejects (\(p = 0.314\)), and the Brown-Forsythe statistic, \(t_{BF} = 0.068\), is close to zero, with the two groups now equally spread. A square root compresses large distances more than small ones, so it shrinks the spread of the fast cars and their long upper tail at the same time. The model for \(\sqrt{\text{dist}}\) satisfies the assumptions much better.