One predictor is rarely enough. The world is multivariate: college GPA depends on high-school GPA and study time and sleep and a dozen other things. Multiple regression (often called MRC, for Multiple Regression/Correlation) is simply the line from the last chapter grown a dimension:
\[y = b_1 x_1 + b_2 x_2 + \cdots + b_k x_k + a\]
The logic does not change one bit. We still minimize the sum of squared residuals, and we still quantify our uncertainty about every coefficient. What does change, and what trips up almost everyone, is what a coefficient now means when the predictors are tangled up with each other. That tangle is the story of this chapter.
To see the interesting behavior, we need predictors that are correlated with each other. The cleanest way to manufacture that is to build them from a shared source, then add noise - so we know the truth going in.
set.seed(1936) # Fisher's "The Design of Experiments"
base <- rnorm(100)
y <- base + rnorm(100)
x1 <- base + rnorm(100)
x2 <- base + rnorm(100)
dat <- data.frame(y, x1, x2)
# x1 and x2 are correlated because they share 'base'
round(cor(dat), 2) y x1 x2
y 1.00 0.51 0.56
x1 0.51 1.00 0.57
x2 0.56 0.57 1.00CORRELATIONS /VARIABLES=y x1 x2.using Statistics
cor(Matrix(df))Now fit the two-predictor model:
lm3 <- lm(y ~ x1 + x2, data = dat)
coef(lm3)(Intercept) x1 x2
-0.06525968 0.27246178 0.38880552 REGRESSION /STATISTICS COEFF R ANOVA /DEPENDENT y /METHOD=ENTER x1 x2.using GLM
coeftable(lm(@formula(y ~ x1 + x2), df))Here is the puzzle that drives the chapter. If you regress y on x1 alone, you get one slope. In the joint model above you get a different slope. Why? Because in the joint model, b1 is the effect of x1 with x2 held constant - it has been partialled. R did not lie to you in either model; the two slopes answer two different questions.
When predictors overlap, there are three distinct quantities, and sloppy language treats them as one. Picture the classic Venn diagram from Cohen, Cohen, West, and Aiken (Aiken 2002): circle \(Y\) overlaps circles \(X_1\) and \(X_2\), and \(X_1\) and \(X_2\) overlap each other. Label the piece of \(Y\) that only \(X_1\) explains as area \(a\), and the piece \(Y\) has left over (unexplained) as area \(e\).
lm() gives you.All three are built from the same three correlations, but they are not the same number. Here they are, by formula:
r_y1 <- cor(x1, y); r_y2 <- cor(x2, y); r_12 <- cor(x1, x2)
# semipartial correlation of x1 (unique slice = area a)
sr1 <- (r_y1 - r_y2 * r_12) / sqrt(1 - r_12^2)
# partial correlation of x1 (a / (a + e))
pr1 <- (r_y1 - r_y2 * r_12) / sqrt((1 - r_y2^2) * (1 - r_12^2))
# standardized regression coefficient (beta) for x1
beta1 <- (r_y1 - r_y2 * r_12) / (1 - r_12^2)
round(c(semipartial = sr1, partial = pr1, beta = beta1), 3)semipartial partial beta
0.234 0.282 0.284 Notice they share the same numerator and differ only in the denominator - which is exactly why they get muddled. Let us confirm that \(\beta_1\) is really the standardized slope R reports:
lm3z <- lm(scale(y) ~ scale(x1) + scale(x2), data = dat)
round(coef(lm3z)[["scale(x1)"]], 3) # compare to beta1 above[1] 0.284* REGRESSION prints the standardized Beta by default.
REGRESSION /STATISTICS COEFF R ANOVA /DEPENDENT y /METHOD=ENTER x1 x2.using GLM, Statistics
df.zy = (df.y .- mean(df.y)) ./ std(df.y)
df.zx1 = (df.x1 .- mean(df.x1)) ./ std(df.x1)
df.zx2 = (df.x2 .- mean(df.x2)) ./ std(df.x2)
coeftable(lm(@formula(zy ~ zx1 + zx2), df))Same number. So the coefficient lm() reports is a partial effect: it controls for the other predictors. That is the single most important sentence in multiple regression. When someone says “controlling for X,” this is the machinery they are invoking.
Three questions, three numbers. Do not let a textbook (or a reviewer) tell you they are interchangeable.
Statistics is about variance, and \(R^2\) is the cleanest variance story we have. \(R\) is just the correlation between what we observed and what the model predicts; \(R^2\) is the ratio of predicted-to-total variance:
R <- cor(y, predict(lm3))
R2 <- var(predict(lm3)) / var(y)
round(c(R = R, R2 = R2, R2_from_summary = summary(lm3)$r.squared), 3) R R2 R2_from_summary
0.607 0.368 0.368 That is it: \(R^2 = \frac{\text{var}(\hat{y})}{\text{var}(y)}\), the fraction of Y’s variance the model reproduces. But \(R^2\) has a flaw: add any predictor, even pure noise, and it never goes down. So we penalize it for complexity to get the adjusted (or “shrunken”) \(R^2\):
\[R^2_{adj} = 1 - (1 - R^2)\,\frac{n - 1}{n - k - 1}\]
n <- length(y); k <- 2
R2_adj <- 1 - (1 - R2) * (n - 1) / (n - k - 1)
round(c(R2_adj = R2_adj, adj_from_summary = summary(lm3)$adj.r.squared), 3) R2_adj adj_from_summary
0.355 0.355 The adjustment is why you should never judge a model by raw \(R^2\) alone. Throw in enough predictors and you can make \(R^2\) look impressive while the model is worthless.
As always, we do not trust the machine until we can reproduce it. The standard error of a coefficient is:
\[SE_{b_1} = \frac{s_y}{s_{x_1}} \cdot \sqrt{\text{VIF}_1} \cdot \sqrt{\frac{1 - R^2}{n - k - 1}}\]
where the variance inflation factor \(\text{VIF}_1 = \frac{1}{1 - R^2_1}\), and \(R^2_1\) comes from regressing \(x_1\) on all the other predictors. Watch:
R2_1 <- summary(lm(x1 ~ x2, data = dat))$r.squared # x1 explained by x2
VIF1 <- 1 / (1 - R2_1)
SEb1 <- (sd(y) / sd(x1)) * sqrt(VIF1) * sqrt((1 - R2) / (n - k - 1))
b1 <- coef(lm3)[["x1"]]
t1 <- b1 / SEb1
p1 <- 2 * pt(-abs(t1), df = n - k - 1) # df = n - k - 1, NOT n - 1
round(c(SE = SEb1, t = t1, p = p1), 4) SE t p
0.0941 2.8941 0.0047 summary(lm3)$coefficients["x1", ] # compare Estimate Std. Error t value Pr(>|t|)
0.272461785 0.094143554 2.894109837 0.004696413 The by-hand values match R’s table. Two things to burn in:
For the whole model, the omnibus test is an \(F\) built from \(R^2\):
\[F = \frac{R^2 / k}{(1 - R^2)/(n - k - 1)} \quad \text{with } df_1 = k,\ df_2 = n - k - 1\]
Fmodel <- (R2 / k) / ((1 - R2) / (n - k - 1))
pF <- pf(Fmodel, df1 = k, df2 = n - k - 1, lower.tail = FALSE)
round(c(F = Fmodel, p = pF), 5) F p
28.2781 0.0000 summary(lm3)$fstatistic # compare F and its df value numdf dendf
28.2781 2.0000 97.0000 Let us watch the VIF at work directly. As two predictors become more correlated, the VIF explodes:
curve(1 / (1 - x^2), from = 0, to = 0.99, n = 500,
xlab = "correlation between X1 and X2", ylab = "VIF",
main = "Collinearity inflates the variance of your estimates")
abline(h = 10, col = "red", lty = 2) # a common (rough) rule-of-thumb threshold
Past a correlation of about 0.95 the curve rises very steeply. A predictor that is 95% predictable from the others carries almost no unique information, so the model cannot pin down its slope. The red line marks a rough \(\text{VIF} = 10\) rule of thumb, but treat it as a warning sign, not a hard rule.
The VIF told us that collinearity inflates our uncertainty. Here is the deeper reason, and it is a matter of bookkeeping: when predictors overlap, the variance they jointly explain cannot be cleanly assigned to one or the other. How you divide it up is a choice, and there are three main conventions, the types of sums of squares.
Type I (sequential). Each predictor is credited with the variance it explains after the predictors entered before it. It is order-dependent, and the pieces add up to the model total. Watch the credit shift when we swap the order:
set.seed(42); n <- 200
x1 <- rnorm(n); x2 <- 0.7 * x1 + rnorm(n, sd = sqrt(1 - 0.49)) # x1 and x2 correlate ~ .7
y <- 0.5 * x1 + 0.4 * x2 + rnorm(n)
anova(lm(y ~ x1 + x2))[1:2, "Sum Sq"] # x1 first, then x2[1] 92.77167 22.07744anova(lm(y ~ x2 + x1))[1:2, "Sum Sq"] # x2 first, then x1 — same data, different split[1] 100.0556 14.7935* Type I (sequential) sums of squares; DESIGN order sets the entry order.
UNIANOVA y WITH x1 x2 /METHOD=SSTYPE(1) /DESIGN=x1 x2.
UNIANOVA y WITH x1 x2 /METHOD=SSTYPE(1) /DESIGN=x2 x1.using GLM, StatsModels
# R's anova() gives Type I (sequential) sums of squares, so the split
# depends on term order; swap the formula terms to change the credit.
anova(lm(@formula(y ~ x1 + x2), df))
anova(lm(@formula(y ~ x2 + x1), df))Same data, same model, but x1 is credited with far more variance when it enters first. Whoever goes first gets the shared variance. Sequential sums of squares reward order, which is appropriate when theory dictates an order and misleading when it does not.
Type III (unique, or “each last”). Each predictor is credited only with the variance that is uniquely its own, the part no other predictor can account for. This is order-free: every predictor is evaluated as if it had entered last. It is exactly the squared semipartial correlation from earlier in this chapter, and you can read it off as the drop in fit when you remove that one predictor from the full model:
full <- lm(y ~ x1 + x2)
c(x1_unique = anova(lm(y ~ x2), full)[["Sum of Sq"]][2], # what x1 adds, controlling x2
x2_unique = anova(lm(y ~ x1), full)[["Sum of Sq"]][2]) # what x2 adds, controlling x1x1_unique x2_unique
14.79350 22.07744 * Type III (unique / "each last") sums of squares; SPSS's default.
UNIANOVA y WITH x1 x2 /METHOD=SSTYPE(3) /DESIGN=x1 x2.using GLM, StatsModels
# Type III (unique) sums of squares: each term controlling for all others.
# SPSS defaults to Type III; R's anova() gives Type I, so build it from
# the full model down (drop one term at a time), as the R tab does.
anova(lm(@formula(y ~ x1 + x2), df))These do not depend on order, but they do not add up to the model total:
c(unique_sum = anova(lm(y ~ x2), full)[["Sum of Sq"]][2] +
anova(lm(y ~ x1), full)[["Sum of Sq"]][2],
model_SS = sum(anova(full)[1:2, "Sum Sq"]))unique_sum model_SS
36.87093 114.84910 The gap is the shared variance - the part both predictors could explain, which Type III refuses to award to either. With collinear predictors that shared piece can be large (here it is most of the model), and that is the collinearity problem stated exactly: much of the explained variance cannot be attributed to one predictor over the other.
Type II sits between the two: like Type III it credits each term controlling for the others, but it respects marginality (a main effect is tested controlling for other main effects, but not for the interactions that contain it). For a model without interactions, Type II and Type III agree.
There is no universally correct type; there is the type that answers your question.
The point to carry away is this: whenever predictors are correlated, some of the explained variance is shared and unattributable, and different conventions handle it differently. Reporting a single “variance explained by X” without saying how you partitioned it is incomplete. Collinearity is not only a nuisance to diagnose with a VIF; it is a real limit on what the data can say about any one predictor, and how you partition the variance is part of the answer.
There is no single fix; there are trade-offs. In rough order of how often they are the right call:
rowMeans(scale(cbind(x1, x2))). Often the most defensible move when the predictors are indicators of one underlying construct.composite <- rowMeans(data.frame(x1z = scale(x1), x2z = scale(x2)))
round(summary(lm(y ~ composite))$coefficients, 3) Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.076 0.074 -1.032 0.303
composite 0.828 0.080 10.288 0.000One predictor, no collinearity, and often a cleaner story than two fighting over the same variance.
Let’s leave the simulated world and put all of this to work on the real GPA dataset from the last chapter. There we predicted college GPA from high-school GPA alone. Now we add the other things an admissions office actually has: SAT total and the quality of the applicant’s letters of recommendation.
gpa <- read.csv("data/gpa.csv")
fit <- lm(CollegeGPA ~ HSGPA + SATtotal + QLOR, data = gpa)
round(summary(fit)$coefficients, 4) Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1533 0.3229 -0.4746 0.6362
HSGPA 0.3764 0.1143 3.2938 0.0014
SATtotal 0.0012 0.0003 4.0464 0.0001
QLOR 0.0227 0.0510 0.4449 0.6574c(R2 = summary(fit)$r.squared, adj_R2 = summary(fit)$adj.r.squared) R2 adj_R2
0.3997354 0.3809771 REGRESSION /STATISTICS COEFF R ANOVA
/DEPENDENT CollegeGPA /METHOD=ENTER HSGPA SATtotal QLOR.using GLM
coeftable(lm(@formula(CollegeGPA ~ HSGPA + SATtotal + QLOR), df))Every coefficient here is a partial effect - the pull of that predictor with the other two held constant. Notice HSGPA’s slope shrank from its bivariate value: some of what looked like “high-school GPA” was really SAT and letter quality traveling alongside it. That is partialling doing its job.
Now check for collinearity by hand, exactly as we defined it - each VIF is \(1/(1-R^2_j)\) where \(R^2_j\) regresses that predictor on the others:
preds <- c("HSGPA", "SATtotal", "QLOR")
vif <- sapply(preds, function(p) {
others <- setdiff(preds, p)
r2 <- summary(lm(as.formula(paste(p, "~", paste(others, collapse = "+"))), data = gpa))$r.squared
1 / (1 - r2)
})
round(vif, 2) HSGPA SATtotal QLOR
1.94 1.24 1.66 All three VIFs sit near 1 - our predictors carry mostly distinct information, so the model can pin down each slope cleanly. (Recall the HSGPA–QLOR correlation was the highest at 0.63; that is real but nowhere near the danger zone.) This is the whole workflow in miniature: fit, read the partial effects, and check that collinearity has not quietly wrecked your standard errors.
Simulate three predictors where x3 is nearly a copy of x1 (e.g., x3 <- x1 + rnorm(100, 0, 0.2)), plus an outcome y that truly depends on x1 and x2.
y ~ x1 + x2 + x3 and report the coefficients, their standard errors, and each VIF (compute the VIFs by hand with 1/(1 - R2_j)).x3 and refit. What happened to the standard error and \(t\) for x1? Explain why in terms of the VIF sitting inside the SE.x2 in the full model, and confirm its square equals the drop in \(R^2\) when you remove x2.If you can explain part 2 to a skeptical colleague, you understand multicollinearity better than the software does.
We have been treating every predictor as continuous. But some of the most important predictors are categorical - treatment vs. control, cat vs. dog, before vs. after. It turns out ANOVA is not a separate technique at all; it is just multiple regression with cleverly coded categorical predictors. That unification - the General Linear Model - is the subject of the next chapter, ANOVA and the GLM.