Multiple linear regression, and what controlling for another variable really does

CS1 · Actuarial Statistics · 13 min read

The problem it solves

Fit marks against hours of study and you get a slope: so many marks per extra hour. Now notice that the students who studied least were also the ones sitting the paper for the second time. Is the slope measuring the effect of study, or partly the effect of being a resitter?

This is confounding, and a single-predictor regression cannot see it. Multiple regression can, by fitting both predictors at once and reporting what each contributes with the other held constant. Almost every actuarial model is of this shape — claim frequency against age and vehicle group and region, all at the same time — and so is the generalised linear model that follows it in the syllabus.

And the headline fact, which surprises students the first time they see it: adding a second predictor changes the coefficient on the first. Not its standard error, its actual value. In the worked example the coefficient on study hours drops from 2.8452381 to 2.0698690 when past attempts joins the model. Neither number is wrong; they answer different questions.

The model, and the equations that fit it

The errors are assumed independent, with mean zero and the same variance throughout, and normal if you want tests and intervals rather than just estimates. Least squares chooses the coefficients to minimise the sum of squared residuals, and differentiating that sum gives one equation per coefficient — the normal equations. For two predictors, written in terms of corrected sums of squares and products, they are

where S12 = Σ(x1 − x̄1)(x2 − x̄2) and so on. Two linear equations in two unknowns, solved by elimination or by the 2 by 2 inverse:

The intercept comes last and for free, because the fitted plane always passes through the point of means:

A warm-up with round numbers: when the predictors do not interfere

Suppose S12 = 0 — the two predictors are uncorrelated in your data. Then the equations uncouple completely and

which are exactly the two simple regression slopes. With uncorrelated predictors, fitting them together and fitting them separately give identical answers, and nothing is confounded with anything.

That is the reference case to hold in mind, because everything difficult about multiple regression is a consequence of S12 not being zero. The cross-term is the confounding, the shifting coefficients and the collinearity, all at once.

What controlling for a variable actually does

Here is the interpretation that makes the topic click, and it can be verified arithmetically. To get β̂2:

Regress y on x1 alone and keep the residuals — the part of y that x1 cannot explain. Regress x2 on x1 alone and keep those residuals too — the part of x2 that is not just x1 in disguise. Now regress the first set of residuals on the second. The slope you get is exactly β̂2 from the full fit.

In words: a multiple regression coefficient is the effect of the part of that predictor which is unrelated to the others, on the part of the response the others have not already accounted for.

Two consequences follow immediately, and both are examinable. First, a coefficient only means something relative to the other variables in the model — change the model and you change the question. Second, if x2 is almost perfectly predictable from x1, there is almost nothing left of it after the first regression, so its coefficient is estimated from a scrap of information and its standard error explodes. That is multicollinearity, and this is why it hurts.

The ANOVA table, R squared, and the F test

The same decomposition as in analysis of variance, with "explained by the predictors" in place of "between groups":

With n observations and p predictors, the degrees of freedom are p for the regression and n − p − 1 for the residual. The residual mean square estimates the error variance:

and the F statistic MSregression/s² tests whether all the coefficients are zero together, which is a different hypothesis from any individual t test.

R squared is the share of the variation the model explains, SSregression/Syy. Its flaw is fatal if you use it to choose a model: it can never fall when you add a predictor, however irrelevant, because least squares can always set a useless coefficient near zero and do no worse. Adjusted R squared repairs this by penalising parameters:

It can fall, and when it falls the variable you added was not worth its degree of freedom.

Multicollinearity

When predictors are strongly correlated with each other, the fit is unstable in a specific way: the coefficients become badly determined even though the predictions stay fine. The standard measure is the variance inflation factor, which for two predictors is

and it multiplies the variance of a coefficient relative to the uncorrelated case, so the standard error inflates by √VIF. The usual rough signals: a VIF above 5 is worth noting, above 10 is a problem.

The symptoms to recognise in an exam answer: large standard errors, coefficients that move wildly when a variable is added or dropped, a significant F with no significant individual t, and coefficients whose sign is the opposite of what the subject matter says. The cures are to drop one of the offenders, combine them into a single sensible variable, or get data in which they are not so entangled — and note that for prediction within the range of the data, collinearity does surprisingly little damage. It is interpretation that suffers.

Where it is used after this topic

Generalised linear models. Same linear predictor, same design matrix, same idea of a partial effect, with a link function and a non-normal response added. Everything on this page transfers.

Analysis of variance. Coding groups as dummy variables turns ANOVA into a regression, and the F tests coincide.

Reserving and pricing models, where several rating factors act at once and the whole point is to separate their effects.

Where marks are lost

Interpreting a coefficient without the phrase "holding the others constant". Omit it and the sentence is simply false.

Residual degrees of freedom. It is n − p − 1 — one for each predictor and one for the intercept. With 8 observations and 2 predictors there are 5, and using 6 or 7 changes every standard error and every p-value.

Reading a high R squared as a good model. It is a statement about fit to the data in hand, not about the model being right, and it rises automatically with every variable added.

Reading the F test as a test of each coefficient. F tests them jointly. With correlated predictors the two can genuinely disagree: a significant F and two insignificant t values is the classic collinearity signature, not a contradiction.

Extrapolating. A fitted plane is a statement about the region the data occupies, and that region is not a rectangle. A prediction at a combination of predictor values that never occurs in the data — many hours of study and many past attempts, when the two go the other way — is extrapolation even if each value is individually unremarkable.

Questions students ask

Why did my first coefficient change when I added the second variable?

Because it is now answering a different question. Before, it measured the whole association between x1 and y, including whatever part of that association really belongs to x2. Now it measures only the part of the association that x2 cannot account for. The size of the change is a direct measure of how confounded the two predictors are, and it is zero exactly when S12 = 0.

How many predictors can I use?

Fewer than you have observations, as hard arithmetic: with p = n − 1 the fit passes exactly through every point, R squared is 1 and the residual degrees of freedom are zero, so nothing can be tested. Long before that limit, each predictor costs a degree of freedom and widens every interval. The practical answer is that predictors should be justified by the subject matter, not by whether they improve the fit.

Can I compare coefficients to see which predictor matters more?

Not as they stand, because they carry units: marks per hour and marks per attempt are not comparable quantities. Standardising the variables first makes the coefficients comparable in a narrow statistical sense, but "matters more" usually means something about the real world rather than about the data, and the safer answer is to quote the effect of a realistic change in each predictor instead.

A worked example

Eight candidates. For each one, the hours of revision a week, the number of previous attempts at the paper, and the mark. The data is invented for this example; it is not a record of anyone's students.

Notice the two dips, at 8 hours and at 12 hours. Both of those candidates had sat the paper before. That is the pattern the second predictor exists to capture.

The sums, both ways

The means are x̄1 = 88/8 = 11.0, x̄2 = 7/8 = 0.875 and ȳ = 514/8 = 64.25. The corrected sums of squares and products come out as

Every one of those was computed twice, once from the deviations and once from the raw totals Σx1² = 1136, Σx2² = 11, Σx1x2 = 58, Σx1y = 6132, Σx2y = 377 and Σy² = 34572. For instance S1y = 6132 − 8(11)(64.25) = 6132 − 5654 = 478, and S12 = 58 − 8(11)(0.875) = 58 − 77 = −19. The raw route is quicker under exam pressure and the deviation route is easier to check; doing both is how you find a slip before it reaches the answer.

The negative S12 is the confounding, visible before any fitting: the candidates who revised least are the ones who had sat the paper before.

The coefficients

The determinant of the system is 168(4.875) − (−19)² = 819 − 361 = 458, and then

So the fitted model is mark = 47.4803493 + 2.0698690(hours) − 6.8558952(attempts): a little over two marks for each extra hour a week, and nearly seven marks lower for each previous attempt, each holding the other constant.

The check: each coefficient as a simple regression on residuals

Take the claim in the body seriously and verify it. Regress y on hours alone; the slope is 478/168 = 2.8452381. Regress attempts on hours alone; that slope is −19/168 = −0.1130952. Form both sets of residuals, and then regress the y-residuals on the attempts-residuals:

Exactly β̂2, to every digit. Doing the same with the roles of the two predictors swapped returns 2.0698690, exactly β̂1. Two routes, two agreements, and an interpretation earned rather than asserted: "controlling for hours" means "using only the part of attempts that hours cannot predict".

Notice also what the denominator says. S22 is 4.875, but after removing what hours can explain only 2.7261905 of it is left — 56 per cent. The coefficient on attempts is estimated from just over half the information the column appeared to contain. That is the price of the correlation between the predictors, in the currency it is actually charged in.

Fitted values and residuals

The residuals add to zero, and so do their products with each predictor. Those three sums are forced to be zero by the normal equations themselves, so they are not evidence that the model is good — but if they come out non-zero, the arithmetic is wrong, which makes them the fastest check available.

The ANOVA table and R squared

Both coefficients are exact fractions over 458, so keep them that way here. Multiplying the seven-decimal versions by 478 and 72.75 magnifies their rounding into the sixth decimal of the answer, and gives 1488.1637578 instead.

Squaring and adding the eight residuals in the table above gives 59.3362445 as well, so the two routes to the residual sum of squares agree.

The 5 per cent point of an F with 2 and 5 degrees of freedom is 5.7861350, so an F of 62.7 is overwhelming; the p-value is 0.0002879. The residual standard deviation is √11.8672489 = 3.4448874 marks, which is the right way to describe the model's accuracy to someone non-technical: predictions are typically about three and a half marks out.

Check R squared the other way: 1 − 59.3362445/1547.5 = 0.9616567. And a third route, through the F statistic: (R²/2)/((1−R²)/5) = 62.7005 — the F from the table, to the four decimals that a seven-decimal R squared can support. Three agreements.

Is the second predictor worth having?

Three ways of asking, all consistent.

Its t test. The standard errors are se(β̂1) = 0.3554100 and se(β̂2) = 2.0863969, so

On 5 degrees of freedom the 2.5 per cent point is 2.5705818, so both are significant; the p-value for attempts is 0.0218096 and a 95 per cent interval for its coefficient runs from −12.2191492 to −1.4926412. Significant, and still very imprecise — eight observations buy a conclusion about the sign and little more about the size.

The extra sum of squares. Fitting hours alone leaves a residual sum of squares of 187.4761905; adding attempts brings it to 59.3362445, a reduction of 128.1399459 for one degree of freedom. Dividing by the residual mean square gives 128.1399459/11.8672489 = 10.7977803, which is exactly (−3.2859976)². The partial F test and the t test are the same test, again.

R squared. Hours alone explains 0.8788522 of the variation; both together explain 0.9616567. The adjusted figure rises as well, so the predictor has paid for its degree of freedom.

How much damage is the correlation doing?

So the variance of each coefficient is 1.79 times what it would have been with uncorrelated predictors, and the standard error is √1.7882 = 1.3372 times larger. Check that directly: with the predictors orthogonal, se(β̂1) would have been √(11.8672489/168) = 0.2657788, and 0.3554100/0.2657788 = 1.3372. The two agree, which confirms what the variance inflation factor is actually measuring.

A VIF of 1.79 is mild and nothing here needs fixing. It is worth knowing what the same arithmetic would say at r12 = 0.95: a VIF of 10.3 and standard errors more than three times inflated, which is the regime where coefficients start changing sign between samples.

One prediction, and one warning

A candidate revising 10 hours a week with one previous attempt:

with a standard error of 1.2468606 for the mean response, giving a 95 per cent interval of 58.1179868 to 64.5283014. That combination sits comfortably inside the data, so the interval is trustworthy.

Now ask for 18 hours a week with 2 previous attempts. The model answers 47.4803493 + 2.0698690(18) − 6.8558952(2) = 71.0262009 marks without hesitation, and the answer is worthless: nobody in these eight candidates combined long hours with repeated attempts, because the two go in opposite directions in this data. The fitted plane has never been tested in that corner. Extrapolation in multiple regression is not about going past the ends of each variable — it is about going outside the cloud the variables occupy together, and that is far easier to do by accident than the single-predictor case ever was.

Where this fits in CS1

This page explains one idea from CS1 Actuarial Statistics. The full CS1 syllabus, with what each part is examined on, sits on the course page, and past-paper style mock exams with full solutions are where you find out whether the idea has stuck under time pressure.

Other CS1 topics explained

• Moment generating functions, and how to get moments out of them
• The central limit theorem, and when the normal approximation is safe
• Maximum likelihood estimation, step by step
• Confidence intervals: what the 95% actually means
• Hypothesis testing and p-values, without the folklore
• Simple linear regression: least squares, fitted by hand
• Generalised linear models: the three parts, and the link function
• Bayesian estimation, conjugate priors and the credibility formula inside them
• Conditional expectation, and the variance decomposition that comes with it
• Sampling distributions: what the sample mean and sample variance actually do
• The method of moments, and when maximum likelihood is worth the extra work
• The likelihood ratio test, and why the Wald and score tests exist too
• The chi-square goodness of fit test, and how to get its degrees of freedom right
• Contingency tables: testing whether two classifications are related
• Analysis of variance: comparing three or more means at once

Stuck on a term rather than the idea? The actuarial glossary defines the notation used above. If something here is wrong or unclear, tell us — these pages get corrected.