STAC67: Problem Set 5 solutions
Worked answers to every question
STAC67, Fall 2026
Week 6
About this document
These are the answers to every question in Problem Set 5, 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. Least squares with matrices, by hand
(a) \(\boldsymbol{X}'\boldsymbol{X}\) and \(\boldsymbol{X}'\boldsymbol{Y}\)
\[\boldsymbol{Y} = \begin{bmatrix} 2 \\ 4 \\ 7 \\ 7 \\ 13 \end{bmatrix}, \qquad \boldsymbol{X} = \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 4 \\ 1 & 5 \\ 1 & 8 \end{bmatrix},\] \[\boldsymbol{X}'\boldsymbol{X} = \begin{bmatrix} n & \sum x_i \\ \sum x_i & \sum x_i^2 \end{bmatrix} = \begin{bmatrix} 5 & 20 \\ 20 & 110 \end{bmatrix}, \qquad \boldsymbol{X}'\boldsymbol{Y} = \begin{bmatrix} \sum Y_i \\ \sum x_i Y_i \end{bmatrix} = \begin{bmatrix} 33 \\ 177 \end{bmatrix}.\]
(b) The inverse
\(\det(\boldsymbol{X}'\boldsymbol{X}) = 5 \times 110 - 20^2 = 150\), so \[(\boldsymbol{X}'\boldsymbol{X})^{-1} = \frac{1}{150}\begin{bmatrix} 110 & -20 \\ -20 & 5 \end{bmatrix} = \begin{bmatrix} 0.7333 & -0.1333 \\ -0.1333 & 0.0333 \end{bmatrix}.\] With \(\bar x = 4\), \(\sum x_i^2/n = 22\) and \(S_{xx} = 110 - 5 \times 4^2 = 30\), the formula gives \(\frac{1}{30}\begin{bmatrix} 22 & -4 \\ -4 & 1 \end{bmatrix}\), the same matrix (multiply numerator and denominator by 5).
(c) The estimates
\[\hat{\boldsymbol{\beta}} = \frac{1}{150}\begin{bmatrix} 110 \times 33 - 20 \times 177 \\ -20 \times 33 + 5 \times 177 \end{bmatrix} = \frac{1}{150}\begin{bmatrix} 90 \\ 225 \end{bmatrix} = \begin{bmatrix} 0.6 \\ 1.5 \end{bmatrix}.\] Check: \(S_{xY} = \sum x_i Y_i - n\bar x\bar Y = 177 - 5 \times 4 \times 6.6 = 45\), so \(\hat\beta_1 = 45/30 = 1.5\) and \(\hat\beta_0 = 6.6 - 1.5 \times 4 = 0.6\). The fitted line is \(\hat Y = 0.6 + 1.5\,x\).
(d) Fitted values, residuals and \(MSE\)
Each row of \(\boldsymbol{X}\hat{\boldsymbol{\beta}}\) is \(0.6 + 1.5\,x_i\): \[\hat{\boldsymbol{Y}} = \boldsymbol{X}\hat{\boldsymbol{\beta}} = \begin{bmatrix} 2.1 \\ 3.6 \\ 6.6 \\ 8.1 \\ 12.6 \end{bmatrix}, \qquad \boldsymbol{e} = \boldsymbol{Y} - \hat{\boldsymbol{Y}} = \begin{bmatrix} -0.1 \\ 0.4 \\ 0.4 \\ -1.1 \\ 0.4 \end{bmatrix}.\] The residuals add up to zero, as they must. Then \[SSE = \boldsymbol{e}'\boldsymbol{e} = 0.01 + 0.16 + 0.16 + 1.21 + 0.16 = 1.70, \qquad MSE = \frac{SSE}{n - 2} = \frac{1.70}{3} = 0.567.\]
(e) The determinant and equal \(x\) values
\[\det(\boldsymbol{X}'\boldsymbol{X}) = n\sum_i x_i^2 - \Big(\sum_i x_i\Big)^2 = n\Big(\sum_i x_i^2 - n\bar x^2\Big) = n\,S_{xx}.\] Here \(5 \times 30 = 150\). If all the \(x_i\) are equal, \(S_{xx} = 0\), so \(\det(\boldsymbol{X}'\boldsymbol{X}) = 0\): \(\boldsymbol{X}'\boldsymbol{X}\) is singular, and the normal equations have no unique solution.
A singular \(\boldsymbol{X}'\boldsymbol{X}\) means there is a linear relation between the columns of \(\boldsymbol{X}\). A square matrix is singular exactly when it sends some nonzero vector to \(\boldsymbol{0}\), so there is a \(\boldsymbol{v} = (v_0, v_1)' \neq \boldsymbol{0}\) with \(\boldsymbol{X}'\boldsymbol{X}\boldsymbol{v} = \boldsymbol{0}\). Multiplying on the left by \(\boldsymbol{v}'\), \[0 = \boldsymbol{v}'\boldsymbol{X}'\boldsymbol{X}\boldsymbol{v} = (\boldsymbol{X}\boldsymbol{v})'(\boldsymbol{X}\boldsymbol{v}) = \sum_{i=1}^n (v_0 + v_1 x_i)^2 .\] A sum of squares is zero only if every term is zero, so \(\boldsymbol{X}\boldsymbol{v} = \boldsymbol{0}\): \[v_0 \begin{bmatrix} 1 \\ \vdots \\ 1 \end{bmatrix} + v_1 \begin{bmatrix} x_1 \\ \vdots \\ x_n \end{bmatrix} = \boldsymbol{0}, \quad \text{with } v_0, v_1 \text{ not both zero.}\] This is a linear relation between the two columns of \(\boldsymbol{X}\): \(\boldsymbol{X}\) is not of full column rank (its rank is 1, not 2). Here \(v_1 \neq 0\), since \(v_1 = 0\) would leave \(v_0 \boldsymbol{1} = \boldsymbol{0}\) and \(v_0 = 0\). Hence \(x_i = -v_0/v_1\) for every \(i\): the second column of \(\boldsymbol{X}\) is a multiple of the first, which is exactly the case of equal \(x_i\). Conversely, if \(\boldsymbol{X}\boldsymbol{v} = \boldsymbol{0}\) for some \(\boldsymbol{v} \neq \boldsymbol{0}\), then \(\boldsymbol{X}'\boldsymbol{X}\boldsymbol{v} = \boldsymbol{0}\) and \(\boldsymbol{X}'\boldsymbol{X}\) is singular. The same argument holds for a design matrix with any number of columns: \(\boldsymbol{X}'\boldsymbol{X}\) is invertible exactly when \(\boldsymbol{X}\) has full column rank, that is, when no column of \(\boldsymbol{X}\) is a linear combination of the others (Exercise 3).
In terms of the data: if every observation was taken at the same \(x\), the data say nothing about how the mean of \(Y\) changes when \(x\) changes. Every line through the point \((\bar x, \bar Y)\) (which any least-square regression line should do) fits equally well, whatever its slope.
X <- cbind(1, c(1, 2, 4, 5, 8)); Y <- c(2, 4, 7, 7, 13)
t(X) %*% X; t(X) %*% Y## [,1] [,2]
## [1,] 5 20
## [2,] 20 110## [,1]
## [1,] 33
## [2,] 177beta_hat <- solve(t(X) %*% X) %*% t(X) %*% Y; beta_hat## [,1]
## [1,] 0.6
## [2,] 1.5e <- Y - X %*% beta_hat; t(e) %*% e## [,1]
## [1,] 1.7Exercise 2. Scottish hill races
(a) The fitted equation and its coefficients
library(MASS)
fit <- lm(time ~ dist + climb, data = hills)
coef(fit)## (Intercept) dist climb
## -8.99203896 6.21795571 0.01104791\[\widehat{\text{time}} = -8.99 + 6.22\,\text{dist} + 0.0110\,\text{climb}.\]
- Distance. Between two races with the same climb, one of them a mile longer, the record time is on average 6.2 minutes longer for the longer race.
- Climb. Between two races of the same distance, the record time is on average 0.0110 minutes longer per extra foot of climb, that is, about 11 minutes per 1000 feet.
- Intercept. It would be the mean record time of a race with no distance and no climb. No such race exists: the races run from 2 to 28 miles and climb from 300 to 7500 feet, so this point is outside the scope of the model, and a negative time is impossible anyway. The intercept only positions the (fitted regression) plane.
(b) The slope of climb, alone and with distance
coef(lm(time ~ climb, data = hills))## (Intercept) climb
## 12.69917062 0.02488635cor(hills$dist, hills$climb)## [1] 0.6523461Alone, the slope of climb is 0.0249, about 25 minutes
per 1000 feet: more than twice the 11 minutes of (a).
The two slopes answer different questions. Distance and climb are
positively correlated (\(r = 0.65\)):
races that climb more are also longer (in terms of distance). In the
simple regression, comparing two races that differ by 1000 feet of climb
also compares, on average, a longer race with a shorter one, so the
slope of climb collects part of the effect of distance. In
the multiple regression the comparison is made at a given distance, and
only the effect of the climb remains.
(c) Three simple regressions
e_time <- resid(lm(time ~ dist, data = hills)) # time, with distance removed
e_climb <- resid(lm(climb ~ dist, data = hills)) # climb, with distance removed
coef(lm(e_time ~ e_climb))[2]## e_climb
## 0.01104791The slope is 0.01105, the coefficient of climb in the
multiple regression. This is what “holding distance fixed” means: \(\hat\beta_{\text{climb}}\) relates the part
of the time that distance does not explain to the part of the climb that
distance does not explain (informally, because “a bit of”
distance is “contained” in climb already).
(d) The same numbers from the matrices
X <- cbind(1, hills$dist, hills$climb); Y <- hills$time
dim(X)## [1] 35 3\(\boldsymbol{X}\) has \(n = 35\) rows, one per race, and \(p = 3\) columns: a column of 1s and one column per predictor.
n <- nrow(X); p <- ncol(X)
XtX_inv <- solve(t(X) %*% X)
beta_hat <- XtX_inv %*% t(X) %*% Y
e <- Y - X %*% beta_hat
MSE <- sum(e^2) / (n - p)
s2_beta <- MSE * XtX_inv
cbind(manual = drop(beta_hat), lm = coef(fit))## manual lm
## (Intercept) -8.99203896 -8.99203896
## dist 6.21795571 6.21795571
## climb 0.01104791 0.01104791c(manual = sqrt(MSE), lm = sigma(fit))## manual lm
## 14.67552 14.67552s2_beta## [,1] [,2] [,3]
## [1,] 18.513523217 -1.2606623775 -1.580489e-03
## [2,] -1.260662378 0.3613787787 -8.042704e-04
## [3,] -0.001580489 -0.0008042704 4.206156e-06vcov(fit)## (Intercept) dist climb
## (Intercept) 18.513523217 -1.2606623775 -1.580489e-03
## dist -1.260662378 0.3613787787 -8.042704e-04
## climb -0.001580489 -0.0008042704 4.206156e-06All three agree. \(SSE = 6891.9\) on
\(n - p = 32\) degrees of freedom, so
\(MSE = 215.4\) and \(\sqrt{MSE} = 14.68\) minutes, the “Residual
standard error” of summary(fit).
The standard error of the coefficient of climb is the
square root of the third diagonal entry of \(\boldsymbol{s^2}\{\hat{\boldsymbol{\beta}}\}\):
\[s\{\hat\beta_{\text{climb}}\} = \sqrt{4.206
\times 10^{-6}} = 0.00205,\] that is, about 2.1 minutes per 1000
feet, for an estimate of 11.0.
(e) The residuals
d <- data.frame(fitted = fitted(fit), resid = resid(fit), race = rownames(hills))
ggplot(d, aes(fitted, resid)) +
geom_hline(yintercept = 0, colour = "grey60") +
geom_point(colour = "grey25") +
geom_text(data = subset(d, resid > 25), aes(label = race), hjust = -0.1, size = 3) +
xlim(0, 200) +
labs(x = "Fitted value (minutes)", y = "Residual (minutes)")
Most residuals lie between \(-16\) and \(10\) minutes. Two races stand far above the rest.
- Knock Hill, with a residual of 65 minutes. It is a
race of 3 miles with 350 feet of climb, and a record of 78.65 minutes,
when the model predicts 13.5. The other races of at most 3 miles have
records between 16 and 19 minutes. This looks like an error in the data
rather than a slow race: the help page of
hillsreports that a later article “corrects” the time to 18.65 minutes. - Bens of Jura, with a residual of 31 minutes. It has the largest climb of all the races (7500 feet), and its record is half an hour longer than the regression plane predicts.
Without Knock Hill (removing it from the dataset and reproducing the
regression analysis), \(\sqrt{MSE}\)
drops from 14.7 to 8.8 minutes, and the coefficient of
climb moves from 11.0 to 11.9 minutes per 1000 feet. How to
detect such observations, and what to do with them, comes later in the
course.
Among the other races, the residuals are mostly negative for fitted values between 50 and 100 minutes: a hint that a plane may not describe the mean record time completely and a more complex model might give a better fit.
Exercise 3. Can the coefficients be estimated?
(a) The five models
1. Celsius and Fahrenheit: no. The columns of \(\boldsymbol{X}\) are the column of 1s, \(C\) and \(F\), and \[F = 32 \times \boldsymbol{1} + 1.8 \times C.\] Neither temperature column is a multiple of the other, but with the column of 1s the three columns are linearly dependent. In terms of the data: once the temperature in Celsius is known, the temperature in Fahrenheit brings nothing new, and “the effect of \(F\) with \(C\) held fixed” has no meaning, since \(F\) cannot change when \(C\) is fixed.
set.seed(67)
tempC <- round(runif(20, 5, 30), 1); tempF <- 32 + 1.8 * tempC
y <- 200 + 12 * tempC + rnorm(20, 0, 25)
coef(lm(y ~ tempC + tempF))## (Intercept) tempC tempF
## 201.91804 11.68326 NAR notices the dependence, drops the last column involved and reports
NA for it.
2. Sand, silt and clay: no. For every sample, \[\text{sand} + \text{silt} + \text{clay} = 100 \times \boldsymbol{1},\] so the column of 1s is a linear combination of the three predictors. Knowing two of the percentages gives the third: one of the three has to be left out of the model.
3. \(x\) and \(x^2\): yes. \(x^2\) is a function of \(x\), but not a linear one, and only linear relations between columns make \(\boldsymbol{X}'\boldsymbol{X}\) singular. A relation \(a + b\,x_i + c\,x_i^2 = 0\) holding for every observation would give the polynomial \(a + bx + cx^2\) three different roots or more, which is impossible unless \(a = b = c = 0\). This is why the condition on the data matters: if \(x\) took only two values, \(u\) and \(v\), then \(x_i^2 = (u + v)\,x_i - uv\) for every observation, and the model could not be fitted.
4. Two predictors and two observations: no. \(\boldsymbol{X}\) has two rows and \(p = 3\) columns. Three columns with two entries each are always linearly dependent, so the \(3 \times 3\) matrix \(\boldsymbol{X}'\boldsymbol{X}\) is singular. In terms of the data: the two observations are two points in space, and infinitely many planes pass through both of them, each with \(SSE = 0\). Least squares cannot choose between them. In general \(n \geq p\) is needed, and \(n > p\) to estimate \(\sigma^2\), since \(MSE = SSE/(n - p)\).
The figure below shows two observations, \((x_1, x_2, Y) = (1, 2, 3)\) and \((3, 1, 6)\), from two angles. The blue plane \(\hat Y = 5.25 + 0.75\,x_1 - 1.5\,x_2\) and the orange plane \(\hat Y = -2.25 + 2.25\,x_1 + 1.5\,x_2\) both pass through the two points, so both give \(SSE = 0\). In fact, for every number \(t\), the plane \[\hat Y = (1.5 - 2.5\,t) + (1.5 + 0.5\,t)\,x_1 + t\,x_2\] passes through both points (\(t = -1.5\) gives the blue plane, \(t = 1.5\) the orange one). Each of these planes contains the black line through the two points: turning a plane around this line gives infinitely many least squares fits.
5. Correlation 0.999: yes. A correlation of 0.999 is not 1: no exact linear relation holds, \(\boldsymbol{X}'\boldsymbol{X}\) can be inverted and the estimates are unique. But \(\boldsymbol{X}'\boldsymbol{X}\) is close to a singular matrix, so its inverse has very large entries, and with them the variances in \(\sigma^2(\boldsymbol{X}'\boldsymbol{X})^{-1}\). The data hardly separate the effects of the two predictors.
set.seed(67); m <- 50
x1 <- rnorm(m); eps <- rnorm(m)
x2 <- x1 + rnorm(m, 0, 0.045) # nearly a copy of x1
x3 <- rnorm(m) # unrelated to x1
round(c(cor(x1, x2), cor(x1, x3)), 3)## [1] 0.999 -0.099round(summary(lm(1 + x1 + x2 + eps ~ x1 + x2))$coefficients[, 1:2], 2)## Estimate Std. Error
## (Intercept) 0.99 0.15
## x1 -0.11 3.18
## x2 2.12 3.19round(summary(lm(1 + x1 + x3 + eps ~ x1 + x3))$coefficients[, 1:2], 2)## Estimate Std. Error
## (Intercept) 0.99 0.15
## x1 0.99 0.15
## x3 0.88 0.14In both simulated models the two true slopes are equal to 1. With nearly identical predictors, their standard errors are about twenty times larger, and the two estimates are far from 1. Their sum is close to 2 all the same: the data give the combined effect of the two predictors, not how it splits between them. This is called multicollinearity, and it comes back later in the course.
(b) Which temperature column to remove
For the fitted values, it does not matter. Since \(C = (F - 32)/1.8\), \[\beta_0 + \beta_1\,C = \Big(\beta_0 - \frac{32}{1.8}\,\beta_1\Big) + \frac{\beta_1}{1.8}\,F.\] Every line in \(C\) is a line in \(F\), and conversely: the two models describe the same set of possible mean functions, so least squares selects the same fitted values, the same residuals and the same \(MSE\) in both.
For the coefficients, it does. The slope in Fahrenheit is the slope in Celsius divided by 1.8, because a degree Fahrenheit is a smaller step. The intercepts differ too: one is the mean response at 0 °C, the other at 0 °F.
fit_C <- lm(y ~ tempC); fit_F <- lm(y ~ tempF)
coef(fit_C); coef(fit_F)## (Intercept) tempC
## 201.91804 11.68326## (Intercept) tempF
## -5.784281 6.490698all.equal(fitted(fit_C), fitted(fit_F))## [1] TRUEc(coef(fit_C)[2] / 1.8, coef(fit_C)[1] - 32 * coef(fit_C)[2] / 1.8)## tempC (Intercept)
## 6.490698 -5.784281