Regression lab
Open up lm() and see what’s inside
Everything on this page runs real R in your browser. Nothing to install, and you can’t break anything: reload the page to start fresh.
How the code boxes work. Edit the code, then press Run Code (or Ctrl/Cmd + Enter). In exercises, replace each ______ with your answer and press Submit Answer to have it checked. Stuck? Use the Hint button first, and the Solution button only after a real attempt.
Warm-up: meet the data
A small made-up dataset of 50 workers: hourly wage, years of education and years of experience. Run the cell to look at it.
we built the data with β_educ = 0.08, so we know the truth
Find the coefficient on educ in the output. It should be near 0.08: each extra year of education is associated with roughly 8% higher wages, holding experience fixed. Now let’s rebuild that number from scratch.
Exercise 1: OLS by the formula
The lesson gave you \hat\beta = (X^\top X)^{-1} X^\top y. In R, t(X) is the transpose, %*% is matrix multiplication, and solve(A) is the inverse of A. The matrix X and vector y are already built for you.
Read the formula from left to right: the inverse of t(X) %*% X, times t(X), times y. Each piece needs %*% between it and the next.
X <- cbind(1, d$educ, d$exper)
y <- log(d$wage)
beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y
beta_hatIn practice solve(t(X) %*% X, t(X) %*% y) is faster and numerically safer, because it avoids computing an inverse explicitly.
Exercise 2: the right angle
The lesson claimed residuals are always at a right angle to every regressor, so X^\top \hat\varepsilon = 0. Compute the residuals and check the claim.
A residual is the gap between what happened and what the line predicts: y minus y_hat.
y_hat <- X %*% beta_hat
e <- y - y_hat
t(X) %*% eExercise 3: R² by hand
R^2 = 1 - \text{RSS}/\text{TSS}. Compute it and compare with what summary() reports.
sum((y - mean(y))^2)
e <- residuals(fit)
RSS <- sum(e^2)
TSS <- sum((y - mean(y))^2)
1 - RSS / TSSExercise 4: heteroskedasticity and robust standard errors
The slides’ example is y_i = 1 + 3x_i + x_i\,\varepsilon_i: the bigger |x_i|, the noisier y_i. The usual standard-error formula assumes the noise is the same everywhere, so here it’s wrong. Run the cell to see the fan shape.
The heteroskedasticity-robust (HC0, or “sandwich”) variance estimate replaces \hat\sigma^2 (X^\top X)^{-1} with
\hat V = (X^\top X)^{-1}\Big(\sum_i X_i X_i^\top \hat\varepsilon_i^2\Big)(X^\top X)^{-1}.
The two slices of bread are the same (X^\top X)^{-1}; the filling lets each observation bring its own error size. Fill in the filling.
diag(v) puts the vector v on the diagonal of a matrix. You want the squared residuals there: e^2.
meat <- t(X) %*% diag(e^2) %*% X
V_robust <- bread %*% meat %*% bread
c(classical = summary(fit)$coefficients["x", "Std. Error"],
robust = sqrt(V_robust[2, 2]))In real work you’d use sandwich::vcovHC(fit, type = "HC0") (or HC1–HC3, which make small-sample corrections), together with lmtest::coeftest().
Does it matter? Check by simulation
A standard error is “right” if the 95% interval built from it covers the truth 95% of the time. Run this to compare the two methods over 1,000 simulated datasets. It takes a few seconds.
aim: 0.95
The classical intervals cover the truth noticeably less than 95% of the time: they claim more certainty than the data justify. The robust ones come close to 95%.
Exercise 5: collinearity inflates variance
When two regressors move together, the regression struggles to separate their effects, and standard errors blow up. The lesson’s formula is
\text{var}(\hat\beta_j) = \sigma^2\,\frac{1}{1 - R_j^2}\,\frac{1}{\sum_i (x_{ij} - \bar x_j)^2},
where R_j^2 comes from regressing x_j on the other regressors. The factor 1/(1-R_j^2) is the variance inflation factor (VIF). Compute it for x1.
1 / (1 - R2_1)
R2_1 <- summary(lm(x1 ~ x2))$r.squared
VIF <- 1 / (1 - R2_1)
VIFA common rule of thumb flags VIFs above 10. Collinearity doesn’t bias the estimates; it makes them imprecise. Often the honest conclusion is “the data can’t separate these two effects”.
Exercise 6: a t-test is a regression
Two groups, A and B. A classic two-sample t-test (with equal variances) asks whether their means differ. The lesson claimed a regression on a group dummy is the same test. Complete the t-test call so it assumes equal variances, and compare its t statistic with the regression’s.
The argument takes a logical value: TRUE or FALSE.
reg <- summary(lm(y ~ group))$coefficients
tt <- t.test(y ~ group, var.equal = TRUE)
c(regression_t = reg["groupB", "t value"], t_test_t = unname(tt$statistic))Free play
A scratchpad with the wage data loaded. Some ideas: add an interaction educ:exper; compare wage with log(wage) as the outcome and interpret the coefficient each way; plot residuals(fit) against fitted(fit) and look for patterns.
Next lesson: Maximum likelihood.