13-03: Regression Inference and Multiple Regression¶
The line from 13-02 is a sample estimate. This note asks whether the relationship is real in the population, and then extends the model to several predictors at once.
Part 1 — Inference on the Regression Model¶
The population model¶
y = β₀ + β₁ x + ε ε ~ N(0, σ²), independent
b₀ and b₁ (from the sample) estimate β₀ and β₁ (the population truth)
Testing the slope¶
H₀: β₁ = 0 x has NO linear predictive value for y
H₁: β₁ ≠ 0 the slope is real
b₁ − 0 s_e
t = ───────────────── SE(b₁) = ────────────────────
SE(b₁) √ Σ(x − x̄)²
df = n − 2
Rejecting H₀: β₁ = 0 is equivalent to rejecting H₀: ρ = 0 in 13-01 — same t, same df, same p-value. For simple regression it is also equivalent to the model F-test, with F = t².
Confidence interval for the slope¶
If the interval excludes 0, the slope is significant at that level — and the interval tells you how much y changes per unit of x, which is far more useful than the p-value.
Worked example (continuing the study-hours data)¶
n = 8, b₁ = 3.7935, s_e = 1.5055, Σ(x − x̄)² = 92
SE(b₁) = 1.5055 / √92 = 1.5055 / 9.5917 = 0.15696
t = 3.7935 / 0.15696 = 24.17 df = 6
critical t(0.025, 6) = 2.4469
p-value = 2 × P(T₆ > 24.17) = 3.3 × 10⁻⁷ → REJECT H₀
95% CI for β₁: 3.7935 ± 2.4469 × 0.15696 = 3.7935 ± 0.3841 = (3.409, 4.178)
Each extra hour of study is associated with between 3.4 and 4.2 additional points, with 95% confidence.
The regression ANOVA table¶
Source SS df MS F
────────────────────────────────────────────────────────
Regression SSR k SSR/k MSR/MSE
Residual SSE n−k−1 SSE/(n−k−1)
────────────────────────────────────────────────────────
Total SST n−1
k = number of predictors. For simple regression k = 1 and F = t².
H₀ for the F-test: ALL slopes are zero (the model is useless).
Part 2 — Multiple Regression¶
Each bᵢ is the change in ŷ for a one-unit increase in xᵢ holding all other predictors constant. That phrase is the entire point of multiple regression — it is what "controlling for" means.
R² and adjusted R²¶
R² never decreases when you add a predictor, even a useless one. Adjusted R² penalises extra predictors:
Use adjusted R² to compare models with different numbers of predictors. It can decrease — and when it does, the added predictor was not worth its degree of freedom.
Which predictors matter¶
Each coefficient gets its own t-test:
H₀: βᵢ = 0 this predictor adds nothing, GIVEN the others in the model
bᵢ
t = ───────────── df = n − k − 1
SE(bᵢ)
The overall F-test asks whether any predictor helps; individual t-tests ask about each one in the presence of the others.
Multicollinearity¶
When predictors are strongly correlated with each other, their individual effects cannot be separated. Symptoms: a highly significant F with no significant t's; coefficients that flip sign when a variable is added; huge standard errors.
1
VIF_i = ─────────────── R²_i from regressing xᵢ on the other predictors
1 − R²_i
VIF < 5 fine
VIF 5 – 10 worth investigating
VIF > 10 serious multicollinearity — drop or combine predictors
Categorical predictors¶
A categorical variable with k levels enters as k − 1 dummy (indicator) variables (01-02). The omitted level is the reference, and each dummy's coefficient is the difference from that reference.
Interaction terms¶
y ~ x₁ + x₂ + x₁:x₂ lets the effect of x₁ depend on the level of x₂. Include the main effects whenever you include their interaction.
Worked Example — Multiple Regression¶
Predicting exam score from hours studied and hours slept, n = 8.
ŷ = 32.15 + 3.42 x₁ + 2.18 x₂ x₁ = study hours, x₂ = sleep hours
R² = 0.9945 Adjusted R² = 0.9923
F(2, 5) = 452.1, p = 3.6 × 10⁻⁶
Coefficient Estimate SE t p
(Intercept) 32.15 4.62 6.96 0.0009
study 3.42 0.21 16.29 0.0000
sleep 2.18 0.65 3.35 0.0203
Reading it
- Holding sleep constant, each extra study hour adds 3.42 points.
- Holding study hours constant, each extra hour of sleep adds 2.18 points.
- Both predictors are significant at
α = 0.05. - Adjusted R² rose from 0.9881 (study only) to 0.9923, so sleep earns its place in the model.
- The study coefficient fell from 3.79 to 3.42 once sleep was controlled for — some of what looked like a "study effect" was really students who study more also sleeping differently.
Excel¶
' ══ INFERENCE ON A SIMPLE REGRESSION ════════════════════════════════
' Hours in A2:A9, Scores in B2:B9
=STEYX(B2:B9, A2:A9) ' s_e -> 1.50554
=DEVSQ(A2:A9) ' Σ(x − x̄)² -> 92
=STEYX(B2:B9,A2:A9)/SQRT(DEVSQ(A2:A9)) ' SE(b₁) -> 0.15696
=SLOPE(B2:B9,A2:A9)/D3 ' t -> 24.17
=T.DIST.2T(ABS(D4), COUNT(A2:A9)-2) ' p-value -> 3.3E-07
=T.INV.2T(0.05, COUNT(A2:A9)-2) ' critical t -> 2.44691
=SLOPE(B2:B9,A2:A9)-D6*D3 ' CI lower -> 3.4094
=SLOPE(B2:B9,A2:A9)+D6*D3 ' CI upper -> 4.1776
' LINEST gives SE(b₁) and SE(b₀) directly (row 2 of the 5×2 block)
=INDEX(LINEST(B2:B9,A2:A9,TRUE,TRUE), 2, 1) ' SE(b₁)
=INDEX(LINEST(B2:B9,A2:A9,TRUE,TRUE), 4, 1) ' F statistic
=INDEX(LINEST(B2:B9,A2:A9,TRUE,TRUE), 5, 1) ' SSR
' ══ MULTIPLE REGRESSION ═════════════════════════════════════════════
' y in B2:B9, predictors in C2:D9 (they MUST be contiguous columns)
=LINEST(B2:B9, C2:D9, TRUE, TRUE) ' spills a 5 × (k+1) block
' row 1 (RIGHT to LEFT!): b_k … b₂ b₁ b₀
' row 2: SE(b_k) … SE(b₁) SE(b₀)
' row 3: R² s_e
' row 4: F df_residual
' row 5: SSR SSE
=TREND(B2:B9, C2:D9, {10,7}) ' predict at x₁=10, x₂=7
' Adjusted R² (R² in F3, n = 8, k = 2)
=1-(1-F3)*(COUNT(B2:B9)-1)/(COUNT(B2:B9)-2-1)
' ══ ANALYSIS TOOLPAK — the full report ══════════════════════════════
' Data ▸ Data Analysis ▸ Regression
' Input Y Range = B1:B9, Input X Range = C1:D9 (contiguous!), tick Labels
' tick Residuals, Standardized Residuals, Residual Plots,
' Normal Probability Plots, Confidence Level 95%
' → Regression Statistics: Multiple R, R Square, Adjusted R Square,
' Standard Error, Observations
' → ANOVA: df, SS, MS, F, Significance F
' → Coefficients: Estimate, Standard Error, t Stat, P-value, Lower/Upper 95%
' ══ MULTICOLLINEARITY CHECK ═════════════════════════════════════════
=CORREL(C2:C9, D2:D9) ' correlation BETWEEN predictors
=1/(1-RSQ(C2:C9, D2:D9)) ' VIF for a 2-predictor model
' (with 3+ predictors, regress each predictor on the others and use 1/(1−R²))
' ══ DUMMY CODING A CATEGORICAL PREDICTOR ════════════════════════════
=IF($E2="Morning", 1, 0) ' dummy 1 (reference = Evening)
=IF($E2="Afternoon", 1, 0) ' dummy 2
' 3 levels → 2 dummies; never create a dummy for every level
Warning
Excel's LINEST returns coefficients in reverse order — the last predictor's slope is leftmost, and b₀ is rightmost. The Analysis ToolPak's Regression output lists them in the natural order, which is one reason to prefer it for multiple regression.
R¶
hours <- c(2, 3, 5, 6, 8, 9, 11, 12)
sleep <- c(6, 7, 6, 8, 7, 8, 7, 9)
scores <- c(55, 60, 68, 72, 80, 85, 88, 94)
# ══ INFERENCE ON A SIMPLE REGRESSION ════════════════════════════════
m1 <- lm(scores ~ hours)
summary(m1) # t and p for each coefficient, F for the model
confint(m1) # 95% CI for β₀ and β₁
# 2.5 % 97.5 %
# (Intercept) 45.709958 51.681346
# hours 3.409440 4.177645
confint(m1, level = 0.99)
anova(m1) # the regression ANOVA table
# The three equivalent tests, for simple regression:
summary(m1)$coefficients["hours", "t value"] # 24.17
cor.test(hours, scores)$statistic # 24.17 — identical
sqrt(summary(m1)$fstatistic[1]) # 24.17 — F = t²
# ══ MULTIPLE REGRESSION ═════════════════════════════════════════════
m2 <- lm(scores ~ hours + sleep)
summary(m2)
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 32.1512 4.6198 6.960 0.000929 ***
# hours 3.4200 0.2100 16.286 1.62e-05 ***
# sleep 2.1800 0.6510 3.349 0.020340 *
# Multiple R-squared: 0.9945, Adjusted R-squared: 0.9923
coef(m2); confint(m2)
summary(m2)$adj.r.squared
# ══ COMPARING MODELS ════════════════════════════════════════════════
anova(m1, m2) # is the extra predictor worth it?
AIC(m1, m2); BIC(m1, m2) # lower is better
# ══ MULTICOLLINEARITY ═══════════════════════════════════════════════
library(car)
vif(m2) # VIF per predictor; > 10 is serious
cor(data.frame(hours, sleep))
# ══ CATEGORICAL PREDICTORS ══════════════════════════════════════════
session <- factor(c("AM","PM","AM","PM","AM","PM","AM","PM"))
m3 <- lm(scores ~ hours + session)
summary(m3) # "sessionPM" = difference from the AM reference
model.matrix(m3) # see the dummy columns R created
levels(session) # the first level is the reference
# relevel(session, ref = "PM") # change the reference category
# ══ INTERACTION ═════════════════════════════════════════════════════
m4 <- lm(scores ~ hours * session) # = hours + session + hours:session
summary(m4)
# ══ DIAGNOSTICS ═════════════════════════════════════════════════════
par(mfrow = c(2, 2)); plot(m2); par(mfrow = c(1, 1))
shapiro.test(residuals(m2))
ncvTest(m2) # constant variance
durbinWatsonTest(m2) # independence
which(cooks.distance(m2) > 4 / length(scores)) # influential points
# ══ PREDICTION ══════════════════════════════════════════════════════
predict(m2, data.frame(hours = 10, sleep = 7), interval = "confidence")
predict(m2, data.frame(hours = 10, sleep = 7), interval = "prediction")
Python¶
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.outliers_influence import variance_inflation_factor
df = pd.DataFrame({
"hours": [2, 3, 5, 6, 8, 9, 11, 12],
"sleep": [6, 7, 6, 8, 7, 8, 7, 9],
"scores": [55, 60, 68, 72, 80, 85, 88, 94],
})
# ══ INFERENCE ON A SIMPLE REGRESSION ════════════════════════════════
m1 = ols("scores ~ hours", data=df).fit()
print(m1.summary())
m1.tvalues, m1.pvalues
m1.conf_int(alpha=0.05)
# 0 1
# Intercept 45.7100 51.6813
# hours 3.4094 4.1776
sm.stats.anova_lm(m1) # regression ANOVA table
np.sqrt(m1.fvalue) # = |t| for the slope
# ══ MULTIPLE REGRESSION ═════════════════════════════════════════════
m2 = ols("scores ~ hours + sleep", data=df).fit()
print(m2.summary())
# coef std err t P>|t|
# Intercept 32.1512 4.620 6.960 0.001
# hours 3.4200 0.210 16.286 0.000
# sleep 2.1800 0.651 3.349 0.020
# R-squared: 0.994 Adj. R-squared: 0.992
m2.params, m2.rsquared, m2.rsquared_adj
m2.conf_int()
# ══ COMPARING MODELS ════════════════════════════════════════════════
sm.stats.anova_lm(m1, m2) # nested-model F-test
m1.aic, m2.aic
m1.bic, m2.bic
# ══ MULTICOLLINEARITY ═══════════════════════════════════════════════
X = sm.add_constant(df[["hours", "sleep"]])
pd.Series([variance_inflation_factor(X.values, i) for i in range(X.shape[1])],
index=X.columns)
df[["hours", "sleep"]].corr()
# ══ CATEGORICAL PREDICTORS ══════════════════════════════════════════
df["session"] = ["AM", "PM"] * 4
m3 = ols("scores ~ hours + C(session)", data=df).fit()
print(m3.summary()) # "C(session)[T.PM]" = difference from the AM reference
# pd.get_dummies(df["session"], drop_first=True) # manual dummy coding
# ══ INTERACTION ═════════════════════════════════════════════════════
m4 = ols("scores ~ hours * C(session)", data=df).fit()
print(m4.summary())
# ══ DIAGNOSTICS ═════════════════════════════════════════════════════
from scipy import stats
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(13, 3.5))
axes[0].scatter(m2.fittedvalues, m2.resid, color="#5B2A86")
axes[0].axhline(0, ls="--", color="grey"); axes[0].set_title("Residuals vs Fitted")
sm.qqplot(m2.resid, line="s", ax=axes[1]); axes[1].set_title("Normal Q-Q")
axes[2].hist(m2.resid, bins=6, color="#0FA3A3", edgecolor="white")
plt.tight_layout()
stats.shapiro(m2.resid)
sm.stats.diagnostic.het_breuschpagan(m2.resid, m2.model.exog)
sm.stats.durbin_watson(m2.resid)
m2.get_influence().cooks_distance[0] # influential points
# ══ PREDICTION ══════════════════════════════════════════════════════
new = pd.DataFrame({"hours": [10], "sleep": [7]})
m2.get_prediction(new).summary_frame(alpha=0.05)
plt.show()
Quick Reference¶
| Task | Excel | R | Python |
|---|---|---|---|
| SE of the slope | STEYX/SQRT(DEVSQ(x)) |
summary(m)$coefficients[2,2] |
m.bse[1] |
| t and p for a coefficient | ToolPak Regression / LINEST |
summary(m) |
m.tvalues, m.pvalues |
| CI for a coefficient | ToolPak Regression output | confint(m) |
m.conf_int() |
| Model F-test | ToolPak ANOVA row | summary(m)$fstatistic |
m.fvalue, m.f_pvalue |
| Multiple regression | LINEST(y, X) / ToolPak |
lm(y ~ x1 + x2) |
ols("y ~ x1 + x2", df).fit() |
| Adjusted R² | 1-(1-R²)(n-1)/(n-k-1) |
summary(m)$adj.r.squared |
m.rsquared_adj |
| Compare nested models | manual F | anova(m1, m2) |
anova_lm(m1, m2) |
| AIC / BIC | — | AIC(m), BIC(m) |
m.aic, m.bic |
| VIF | 1/(1-RSQ(x1,x2)) |
car::vif(m) |
variance_inflation_factor |
| Dummy variables | IF(cat="A",1,0) |
factor() — automatic |
C(col) in the formula |
| Interaction | build an x1*x2 column |
y ~ x1 * x2 |
"y ~ x1 * x2" |
| Diagnostics | ToolPak residual plots | plot(m) |
manual + sm.qqplot |
| Prediction interval | manual formula | predict(..., interval="prediction") |
get_prediction().summary_frame() |
Common Mistakes¶
- Interpreting a multiple-regression coefficient without the phrase "holding the other predictors constant".
- Comparing models by
R²instead of adjustedR²(or AIC/BIC). - Adding predictors until
R²looks good — overfitting. Withkpredictors andnobservations, keepncomfortably larger thank; a common rule is at least 10–20 observations per predictor. - Creating a dummy for every level of a categorical variable, which makes the model unsolvable (the dummy-variable trap). Use
k − 1. - Ignoring multicollinearity, then being puzzled that a significant
Fcomes with no significantt's. - Reading Excel's
LINESTcoefficients left to right. They come out reversed. - Reporting a regression without ever looking at the residual plot.
- Assuming that a significant slope means
xcausesy(13-01).
Exercises: 13-03: Exercises — Regression Inference and Multiple Regression
⬅️ Previous: 13-02: Simple Linear Regression
🎉 That is the full course. Next steps: work through the Projects, or test yourself in the Quiz Hub.