About this document

These are the answers to every question in Problem Set 6, 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 a regression output

(a) The missing values

The model has an intercept and three slopes: \(p = 4\) parameters, and \(n - p = 21 - 4 = 17\) error degrees of freedom.

  • (a) \(s\{\hat\beta\} = \hat\beta / t^* = 0.7156 / 5.307 = 0.1348\).
  • (b) \(t^* = \hat\beta / s\{\hat\beta\} = -0.1521 / 0.1563 = -0.973\).
  • (c) \(n - p = 17\).
  • (d) From \(R^2\), as in Lecture 12: \[F^* = \frac{R^2/(p - 1)}{(1 - R^2)/(n - p)} = \frac{0.9136/3}{0.0864/17} = 59.9.\]
  • (e) \(R^2_a = 1 - \dfrac{n - 1}{n - p}\,(1 - R^2) = 1 - \dfrac{20}{17} \times 0.0864 = 0.8984\).

(b) Interval for the coefficient of Water.Temp

With \(t(0.975;\ 17) = 2.110\): \[1.2953 \pm 2.110 \times 0.3680 = (0.52,\ 2.07).\] With 95% confidence, at a given air flow and acid concentration, each extra degree of cooling-water temperature goes with an increase in mean stack loss of between 0.52 and 2.07: the plant is less efficient when its cooling water is warmer.

(c) The test for Acid.Conc.

\(H_0: \beta_3 = 0\) against \(H_a: \beta_3 \neq 0\). Here \(|t^*| = 0.97 < 2.110\) (\(p = 0.344\)): we do not reject \(H_0\).

This does not mean that acid concentration has nothing to do with the efficiency of the plant. The test is about acid concentration once air flow and water temperature are in the model: the data show no evidence that it adds anything to them. On its own, acid concentration is correlated with stack loss (\(r = 0.40\)), partly because it is correlated with air flow (\(r = 0.50\)). And not rejecting \(H_0\) is not proof that \(\beta_3 = 0\): with 21 days, a small effect could potentially go undetected.

(d) The overall \(F\) test

\(H_0: \beta_1 = \beta_2 = \beta_3 = 0\) against \(H_a\): at least one of them is not zero. The critical value is \(F(0.95;\ 3,\ 17) = 3.197\). Since \(F^* = 59.9\) is far above it, we reject \(H_0\): the three settings of the plant, taken together, are related to stack loss.

summary(lm(stack.loss ~ Air.Flow + Water.Temp + Acid.Conc., data = stackloss))
## 
## Call:
## lm(formula = stack.loss ~ Air.Flow + Water.Temp + Acid.Conc., 
##     data = stackloss)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.2377 -1.7117 -0.4551  2.3614  5.6978 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -39.9197    11.8960  -3.356  0.00375 ** 
## Air.Flow      0.7156     0.1349   5.307  5.8e-05 ***
## Water.Temp    1.2953     0.3680   3.520  0.00263 ** 
## Acid.Conc.   -0.1521     0.1563  -0.973  0.34405    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.243 on 17 degrees of freedom
## Multiple R-squared:  0.9136, Adjusted R-squared:  0.8983 
## F-statistic:  59.9 on 3 and 17 DF,  p-value: 3.016e-09

Exercise 2. Two days of operation, and one of them is unusual

(a) A typical day

fit <- lm(stack.loss ~ Air.Flow + Water.Temp + Acid.Conc., data = stackloss)
new1 <- data.frame(Air.Flow = 60, Water.Temp = 20, Acid.Conc. = 87)
predict(fit, new1, interval = "confidence")
##       fit      lwr      upr
## 1 15.6898 13.98976 17.38984
predict(fit, new1, interval = "prediction")
##       fit      lwr      upr
## 1 15.6898 8.638884 22.74072
Xh1 <- c(1, 60, 20, 87)
sqrt(drop(t(Xh1) %*% vcov(fit) %*% Xh1))                    # s{Y-hat_h}
## [1] 0.8057782
predict(fit, new1, se.fit = TRUE)$se.fit                     # the same
## [1] 0.8057782

At this point, \(\hat Y_h = 15.69\), the 95% confidence interval for the mean stack loss is \((13.99,\ 17.39)\), and the 95% prediction interval for one day is \((8.64,\ 22.74)\). Since vcov(fit) is \(MSE\,(\boldsymbol{X}'\boldsymbol{X})^{-1}\), \(\boldsymbol{X}_h'\,\)vcov(fit)\(\,\boldsymbol{X}_h\) is \(s^2\{\hat Y_h\}\), and its square root agrees with the standard error returned by predict.

(b) A day with high air flow and cool water

new2 <- data.frame(Air.Flow = 75, Water.Temp = 18, Acid.Conc. = 87)
predict(fit, new2, interval = "confidence")
##        fit    lwr      upr
## 1 23.83383 17.569 30.09866
predict(fit, new2, interval = "prediction")
##        fit      lwr      upr
## 1 23.83383 14.55626 33.11141
sapply(stackloss[, 1:3], range)
##      Air.Flow Water.Temp Acid.Conc.
## [1,]       50         17         72
## [2,]       80         27         93

Each value is within the range of the data: air flow from 50 to 80, water temperature from 17 to 27, acid concentration from 72 to 93. Yet the confidence interval, \((17.57,\ 30.10)\), is 3.7 times wider than at the typical day, and the prediction interval, \((14.56,\ 33.11)\), is also wider.

(c) Hidden extrapolation

pts <- data.frame(Air.Flow = c(60, 75), Water.Temp = c(20, 18), day = c("(a)", "(b)"))
ggplot(stackloss, aes(Air.Flow, Water.Temp)) +
  geom_point(colour = "grey30") +
  geom_point(data = pts, colour = "#C05A16", size = 3, shape = 17) +
  geom_text(data = pts, aes(label = day), colour = "#C05A16", nudge_x = 1.6) +
  labs(x = "Air flow", y = "Water temperature (°C)")

X <- model.matrix(fit)
XtX_inv <- solve(t(X) %*% X)
lev <- function(xh) drop(t(xh) %*% XtX_inv %*% xh)
c(day_a = lev(c(1, 60, 20, 87)), day_b = lev(c(1, 75, 18, 87)))
##      day_a      day_b 
## 0.06172196 0.83818175
round(range(hatvalues(fit)), 3)
## [1] 0.052 0.412

In the data, air flow and water temperature rise together (\(r = 0.78\)): every day with an air flow above 70 had cooling water at 25 °C or more. A day with an air flow of 75 and water at 18 °C is unlike any day observed, although each value taken alone is in range. This is hidden extrapolation.

The numbers say the same. For the 21 observed days, \(\boldsymbol{X}_i'(\boldsymbol{X}'\boldsymbol{X})^{-1}\boldsymbol{X}_i\) ranges from 0.052 to 0.412. At the typical day it is 0.062, among the smallest. At the unusual day it is 0.838, about twice the largest value in the data. Since \(s^2\{\hat Y_h\} = MSE\,\boldsymbol{X}_h'(\boldsymbol{X}'\boldsymbol{X})^{-1}\boldsymbol{X}_h\), the intervals of (b) are wide, and more importantly they rely on the model holding in a region where there are no data to check it. They should not be trusted.

Exercise 3. Adding predictors

(a) \(R^2\) and adjusted \(R^2\)

\(p\) \(R^2 = 1 - SSE/SST\) \(R^2_a = 1 - \frac{n - 1}{n - p}\,\frac{SSE}{SST}\)
Model A 2 0.600 \(1 - \frac{24}{23} \times 0.4 = 0.583\)
Model B 4 0.625 \(1 - \frac{24}{21} \times 0.375 = 0.571\)

\(R^2\) prefers model B; the adjusted \(R^2\) prefers model A. The two soil measurements reduce \(SSE\) by only 10, which does not make up for the two degrees of freedom they cost.

(b) The overall \(F\) test of model B

\(H_0: \beta_1 = \beta_2 = \beta_3 = 0\). With \(SSR = 400 - 150 = 250\) on 3 degrees of freedom and \(MSE = 150/21\): \[F^* = \frac{250/3}{150/21} = 11.67 > F(0.95;\ 3,\ 21) = 3.072:\] we reject \(H_0\). Model B is related to yield, which is no surprise since rainfall alone already explains 60% of the variation.

(c) Why \(SSE\) cannot increase

Every column of \(\boldsymbol{X}_A\) is also a column of \(\boldsymbol{X}_B\), so every vector \(\boldsymbol{X}_A\boldsymbol{\beta}\) of the column space of model A is also in the column space of model B: \(C(\boldsymbol{X}_A) \subset C(\boldsymbol{X}_B)\). Least squares finds the point of the column space closest to \(\boldsymbol{Y}\), and \(SSE\) is the squared distance from \(\boldsymbol{Y}\) to it. The closest point in a larger set is at least as close: \[SSE_B = \min_{\boldsymbol{v} \in C(\boldsymbol{X}_B)} \|\boldsymbol{Y} - \boldsymbol{v}\|^2 \ \le\ \min_{\boldsymbol{v} \in C(\boldsymbol{X}_A)} \|\boldsymbol{Y} - \boldsymbol{v}\|^2 = SSE_A.\] Since \(SST\) does not depend on the model, \(R^2 = 1 - SSE/SST\) can only increase (or stay the same) when predictors are added, even useless ones. This is why the adjusted \(R^2\), which charges for each parameter, is used to compare models of different sizes.

(d) Why large values of \(F^*\) count against \(H_0\)

\(MSE\) estimates \(\sigma^2\) whatever the true coefficients. \(MSR\) estimates \(\sigma^2 + \frac{1}{p - 1}\sum_i (\mu_i - \bar\mu)^2\). Under \(H_0\), every \(\mu_i = \beta_0\), so the sum is zero and both mean squares estimate \(\sigma^2\): \(F^*\) should be close to 1. When \(H_0\) is false, the means \(\mu_i\) differ, the sum is positive, and \(MSR\) tends to exceed \(MSE\). Large values of \(F^*\) are therefore what we expect under \(H_a\), not under \(H_0\), which is why the test rejects for large \(F^*\) only.