These are the answers to every question in Problem Set 3, written out in full. Read it alongside the problem set rather than instead of it, and attempt each question before you look here.
The \(t\) value is the estimate divided by its standard error, and in simple linear regression the \(F\) statistic is its square (Lecture 6): \[t = \frac{10.7296}{0.3148} = 34.08, \qquad F = t^2 = 1162.\]
Each additional minute of eruption is associated with an increase of about 10.73 minutes in the mean waiting time until the next eruption.
\(H_0: \beta_1 = 10\) against \(H_a: \beta_1 \neq 10\). The null value is not zero, so the \(t\) value printed by R does not answer this question: \[t^* = \frac{\hat\beta_1 - 10}{SE(\hat\beta_1)} = \frac{10.7296 - 10}{0.3148} = 2.318 \quad \text{on } 270 \text{ degrees of freedom.}\]
The conclusion depends on the level: there is evidence against the ranger claim but it is not overwhelming (\(p \approx 0.021\)).
\[10.7296 \pm 1.969 \times 0.3148 = \big(10.110,\ 11.349\big).\]
The value 10 lies outside the interval, which is exactly the 5% rejection in (c): a two-sided test at level \(\alpha\) rejects \(\beta_1 = b\) if and only if \(b\) is outside the \((1-\alpha)\) confidence interval.
\(r = +\sqrt{0.8115} = 0.901\), positive because the slope is positive. The regression accounts for about 81% of the variability in waiting times.
From \(SE(\hat\beta_1) = \hat\sigma / \sqrt{S_{xx}}\): \[S_{xx} = \left(\frac{5.914}{0.3148}\right)^2 = 352.93, \qquad s_x^2 = \frac{S_{xx}}{n - 1} = \frac{352.93}{271} = 1.302.\]
var(faithful$eruptions)## [1] 1.302728The small difference comes from the rounding in the printed output.
\(\bar x = 384.3/20 = 19.215\) and \(\bar Y = 2726/20 = 136.30\). \[\hat\beta_1 = \frac{4769.9}{538.5} = 8.8578, \qquad \hat\beta_0 = 136.30 - 8.8578 \times 19.215 = -33.90.\] Fitted line: \(\hat Y = -33.90 + 8.86\, x\).
Each additional degree of maximum temperature is associated with about 8.86 more iced drinks sold, on average.
Literally, \(\hat\beta_0\) says that on a day with a maximum of 0 °C, the shop sells on average -33.9 iced drinks, a negative number. The data cannot support this reading: every day observed had a maximum between 12.0 and 28.4 °C, so \(x = 0\) is far outside the range of the data.
\[SSE = 50082.2 - \frac{4769.9^2}{538.5} = 7831.60, \qquad \hat\sigma^2 = \frac{SSE}{n-2} = 435.09, \qquad \hat\sigma = 20.86.\]
\(\hat Y_{n+1} = -33.902 + 8.8578 \times 25 = 187.54\), and \[SE(\hat Y_{n+1}) = 20.859 \sqrt{\frac{1}{20} + \frac{(25 - 19.215)^2}{538.5}} = 6.985,\] so the interval is \(187.54 \pm 2.101 \times 6.985 = \big(172.87,\ 202.22\big)\).
\[s_{\{pred\}} = 20.859 \sqrt{1 + \frac{1}{20} + \frac{(25 - 19.215)^2}{538.5}} = 21.997,\] so the interval is \(187.54 \pm 2.101 \times 21.997 = \big(141.33,\ 233.76\big)\), about 3.1 times wider.
The interval in (d) is about the average over all 25 °C days, which the data pin down well. This one is about one particular day, which carries its own day-to-day variability, the extra \(1\) under the square root. More data would shrink the first interval towards a point, but never the second below about \(\pm 2\sigma\).
The point prediction is \(-33.902 + 8.8578 \times 38 = 302.69\) drinks. The number itself can be computed, and so can a prediction interval, but 38 °C is almost ten degrees above the hottest day observed. The formulas assume that the straight line still holds out there, and nothing in the data can check that: sales may level off when the shop runs out of stock or staff, or when people stay indoors. The \((x_{n+1} - \bar x)^2\) term widens the interval for sampling uncertainty only, and does not protect against the model itself failing.
| Source | Df | Sum of squares | Mean square | \(F\) |
|---|---|---|---|---|
| Regression | 1 | 7581.8 | 7581.8 | 419.36 |
| Error | 29 | 524.3 | 18.08 | |
| Total | 30 | 8106.1 |
The model uses 1 degree of freedom, and the total is \(1 + 29 = 30 = n - 1\), so there are \(n = 31\) trees. \(SSE = SST - SSR\), and each mean square is its sum of squares over its degrees of freedom.
\[R^2 = \frac{7581.8}{8106.1} = 0.9353, \qquad r = +\sqrt{0.9353} = 0.967,\] positive because \(\hat\beta_1 = 5.066 > 0\). Girth accounts for about 94% of the variability in volume.
\(F^* = 419.36 > F(0.99;\,1,\,29) = 7.60\), so we reject \(H_0\) at the 1% level: there is strong evidence of a linear relationship between girth and volume.
In simple linear regression \(F^* = (t^*)^2\), so \[|t^*| = \sqrt{419.36} = 20.48, \qquad SE(\hat\beta_1) = \frac{\hat\beta_1}{|t^*|} = \frac{5.066}{20.48} = 0.2474.\]
\(\hat\sigma = \sqrt{MSE} = \sqrt{18.08} = 4.252\) cubic feet.
fit <- lm(Volume ~ Girth, data = trees)
anova(fit)## Analysis of Variance Table
##
## Response: Volume
## Df Sum Sq Mean Sq F value Pr(>F)
## Girth 1 7581.8 7581.8 419.36 < 2.2e-16 ***
## Residuals 29 524.3 18.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(fit)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.943459 3.365145 -10.97827 7.621449e-12
## Girth 5.065856 0.247377 20.47829 8.644334e-19Since \(\hat Y_i - \bar Y = \hat\beta_0 + \hat\beta_1 x_i - \bar Y\), \[\sum_i (\hat Y_i - \bar Y)\, e_i = (\hat\beta_0 - \bar Y) \sum_i e_i + \hat\beta_1 \sum_i x_i e_i = 0 + 0 = 0,\] by the two normal equations. Then, writing \(Y_i - \bar Y = e_i + (\hat Y_i - \bar Y)\) and squaring, \[SST = \sum_i e_i^2 + \sum_i (\hat Y_i - \bar Y)^2 + 2 \sum_i (\hat Y_i - \bar Y)\, e_i = SSE + SSR.\]
\(S_{xx}\) is a fixed number, so \[E(SSR) = S_{xx}\, E(\hat\beta_1^2) = S_{xx}\left[\mathrm{Var}(\hat\beta_1) + E(\hat\beta_1)^2\right] = S_{xx}\left[\frac{\sigma^2}{S_{xx}} + \beta_1^2\right] = \sigma^2 + \beta_1^2 S_{xx}.\]
A simulation check, with \(\sigma = 2\) and \(\beta_1 = 0.3\):
set.seed(1)
x <- runif(25, 0, 10); Sxx <- sum((x - mean(x))^2)
ssr <- replicate(5000, {
y <- 1 + 0.3 * x + rnorm(25, 0, 2)
sum((fitted(lm(y ~ x)) - mean(y))^2)
})
c(simulated = mean(ssr), formula = 2^2 + 0.3^2 * Sxx)## simulated formula
## 22.83404 22.72222\(MSR = SSR/1\), so \(E(MSR) = \sigma^2 + \beta_1^2 S_{xx}\), while \(E(MSE) = \sigma^2\).
So a large \(F^*\) is what we expect to see when \(\beta_1 \neq 0\). The term \(\beta_1^2 S_{xx}\) also shows the role of the design: the more spread out the \(x\) values, the larger \(S_{xx}\), and the easier a given slope is to detect. This is the same fact as \(SE(\hat\beta_1) = \sigma / \sqrt{S_{xx}}\) getting smaller.
At \(x_{n+1} = \bar x\) the \((x_{n+1} - \bar x)^2\) term vanishes, so the half-widths are \[t(0.975;\, n-2)\; \hat\sigma \sqrt{\tfrac{1}{n}} \qquad\text{and}\qquad t(0.975;\, n-2)\; \hat\sigma \sqrt{1 + \tfrac{1}{n}}.\] As \(n \to \infty\), \(\hat\sigma \to \sigma\) and \(t(0.975;\,n-2) \to 1.96\). The first half-width tends to \(0\), the second to \(1.96\,\sigma\).
The confidence interval targets a fixed number, the mean response, which more data estimates ever more precisely. The prediction interval must also cover the new observation’s own error, whose standard deviation \(\sigma\) does not shrink with \(n\).
set.seed(67)
n <- 20
x <- runif(n, 1, 10)
beta0 <- 2; beta1 <- 0.5
covers <- function(make_errors, copies = 1, reps = 2000) {
mean(replicate(reps, {
y <- beta0 + beta1 * x + make_errors(x)
xx <- rep(x, copies); yy <- rep(y, copies)
ci <- confint(lm(yy ~ xx), "xx", level = 0.95)
ci[1] <= beta1 && beta1 <= ci[2]
}))
}
normal <- function(x) rnorm(length(x), 0, 1)
funnel <- function(x) rnorm(length(x), 0, 0.25 * x)
skewed <- function(x) rexp(length(x), 1) - 1
cov5 <- c(model_true = covers(normal),
copies_3 = covers(normal, copies = 3),
copies_10 = covers(normal, copies = 10),
funnel = covers(funnel),
skewed = covers(skewed))
round(cov5, 3)## model_true copies_3 copies_10 funnel skewed
## 0.951 0.712 0.424 0.919 0.955The coverage is 0.951. It estimates the probability that the procedure produces an interval containing \(\beta_1\): the long-run proportion of datasets, drawn from this model, whose interval covers the true value. The 95% describes the procedure, not any one interval.
It is not exactly 0.95 because 2000 datasets is a finite sample. The standard error of a simulated proportion near 0.95 is about \(\sqrt{0.95 \times 0.05 / 2000} \approx 0.005\), so anything between roughly 0.94 and 0.96 is consistent with the true 95%.
With three copies the coverage falls to 0.712, and with ten to 0.424. Independence fails: the copies are not new information, but the formulas treat 60 or 200 observations as independent. The estimate \(\hat\beta_1\) is exactly the one from the 20 original points, but the standard error shrinks and the degrees of freedom grow, so the intervals become far too narrow. This is the xkcd cartoon of Lecture 5, in numbers.
The coverage drops to 0.919. Equal variance fails: \(\mathrm{Var}(\hat\beta_1) = \sigma^2 / S_{xx}\) was derived with one common \(\sigma^2\), so the standard error R reports is wrong. Here it is too small, and the intervals miss more often than they should. The estimate itself is still unbiased: only its uncertainty is misjudged.
The coverage is 0.955, essentially 95%. Normality fails, yet the interval still works. \(\hat\beta_1 = \sum k_i Y_i\) averages many errors, so it is close to normal even when each error is not.
Comparing the three: the interval for \(\beta_1\) depends most on independence, then on equal variance, and least on normality. This matches the table in Lecture 7. The exception is the prediction interval, which is about a single new error and so relies on normality much more.