Multiple & Polynomial Regression

Simple regression predicted weight from height alone, and height alone runs out quickly. People of the same height weigh different amounts, and a line that only reads height has to give all of them the same answer.

That is the whole motivation. Some of the difference between people of one height goes with another measured variable, and multiple regression lets several inputs contribute to one prediction.

Take one person before taking a whole dataset. With two inputs the prediction is

ŷ = β₁x₁ + β₂x₂ + α

and it is easiest to read as three separate pieces that get added up.

Below, the three pieces are laid end to end along a number line. The coefficients start at the values the fit chose for the twenty people, and every slider is yours after that. Push the age coefficient negative and watch its block turn round and walk the prediction backwards.

Each input contributes its value multiplied by its coefficient, and the contributions are added. That sentence is the entire model, and everything else on this page is either a way of seeing it or a way of fitting it.

Once the pieces are familiar, group them. Put the inputs in one vector and the coefficients in another,

x = (x₁, x₂)        β = (β₁, β₂)

and the weighted sum collapses to a dot product.

Nothing was added and nothing was hidden. The dot product is compact notation for the weighted sum already calculated by hand, and the linear algebra primer called it the shape of every linear prediction for exactly this reason. With twenty inputs the expression stays three symbols long.

Take the height coefficient first. Hold age at one value, raise height by one centimetre, and read the prediction again.

So a coefficient describes how the prediction changes with one input while the other inputs in the model stay where they are. Any two records that differ in one input only will differ in prediction by exactly that input’s coefficient.

With one input the prediction rule was a line. Give it a second input and the rule needs another direction to change in.

Picture the height-and-weight scatter from the regression page, height across and weight up, with its fitted line. Now let age run into the depth of the page. Every person moves back or forward according to their age, and the line, which had no opinion about depth, extends into a sheet that can tilt along the new direction too. That sheet is the plane in the box at the top of the page, and pressing “face the height axis” up there looks along the age direction so the room flattens back into the familiar scatter.

Adding a second input gives the prediction rule one more direction it can change in. That is all a plane is.

A plane is hard to read as a whole and easy to read one direction at a time. The two flat panels under the room are cuts through it.

Fix an age and the plane seen against height is a straight line, drawn in indigo on the wireframe, and its slope is the height coefficient, 0.651 kg per centimetre. Fix a height and the plane seen against age is another straight line, in amber, with slope 0.154 kg per year. Slide either fixed value and the cut moves across the plane while its slope does not change, because a plane tilts by the same amount everywhere.

That is section 4’s “holding the other input fixed” made literal. A coefficient is the slope of the plane along one axis, and the cut is what you see when you stand at a fixed value of the other axis and look.

The people do not sit on the plane, and that part is familiar. Each person has a measured weight, the plane produces a predicted weight at their height and age, and the vertical gap between them is the residual, exactly as on the regression page.

eᵢ = yᵢ − ŷᵢ

Switch on the residual drops in the box above to see them, green where the person sits above the plane and rose where they sit below it. The residual sum of squares still squares and adds those drops, 214.2 for the twenty people, and the fit is still the plane that makes that total as small as it can be. R² reads 0.865, the share of the weight variation that height and age explain together.

Multiple regression changes the shape of the prediction surface and nothing about how it is judged. Residuals, RSS and the least-squares rule carry over untouched.

It is tempting to imagine the fit working one input at a time, fitting height, locking it, then fitting age on what is left. It does not, and the twenty people show why it cannot.

Changing one coefficient changes what the others need to do, so the fit chooses all of them at once. RSS is a bowl over every coefficient together, the same bowl the regression page drew over a slope and an intercept, with one more direction per input, and the least-squares answer is its one lowest point. Each coefficient is chosen in the context of every other input in the model, which is why adding or removing an input can move every coefficient, not only the one that came or went.

Now the thrown ball. Its height rises and then falls, and a straight line cannot rise and then fall, however it is placed. A model with transformed input columns can fit a curve, though, and it does not need a second measurement to do it. A new input can be calculated from an old one.

Start with time t, and manufacture t² from it.

A model can receive transformed versions of an existing input as extra inputs, and the fit neither knows nor cares that one column was made from another. Polynomial regression is nothing more than this. It builds the power columns and hands the widened table to the machinery of Parts 1 and 2.

The two-input model on those two columns is

ĥ = β₁·t + β₂·t² + α

and it can be read two ways. As multiple regression it has two inputs, t and t², two coefficients and an intercept, and the prediction is still a weighted sum. As a function of time it is curved, since a nonzero coefficient on t² makes a parabola.

The two drawings show the same fit. On the left, time is plotted against height, producing a parabola. On the right, time and squared time are separate coordinates, and the fitted prediction is a plane.

squared-time coordinate = t² = t × t

Actual inputs cannot move freely around that plane: squared time is determined by time. The highlighted path follows that constraint. Reading the plane along this path gives the curve on the left.

The table under the pictures fits the data twice, once as a degree-2 curve in t and once as a plane over two inputs that happen to be called t and t², and the two sets of numbers agree to every digit. They are one fit. The curve in the original input is a plane in the constructed features, which is the page’s central claim.

This deserves its own moment, because “linear” is doing something specific. A model counts as linear regression when its coefficients enter linearly,

ŷ = β₁φ₁(x) + β₂φ₂(x) + α

The feature functions may themselves be nonlinear. Here, one keeps the original time and the other squares it.

φ₁(t) = t
φ₂(t) = t²

The prediction is nonlinear in time but linear in the coefficients being fitted. That second property is what allows the same linear least-squares fitting method to work.

Below, the five basic shapes 1, t, t², t³ and t⁴ are drawn unscaled along the top. Each slider scales one of them, the scaled shapes are dashed in their own colours, and the solid line is their sum. The sliders start on the coefficients recovered from the ideal throw.

Every curve this widget can draw is assembled by adding scaled copies of fixed shapes, and the only things a fit chooses are the amounts. Polynomial regression produces a curved prediction while remaining a linear model in the parameters it learns, and that is why the straight line’s solver handles it without modification.

Physics writes the height of a thrown ball as

h(t) = h₀ + v₀·t − ½·g·t²

with h₀ the starting height, v₀ the vertical launch speed and g the acceleration due to gravity. Line that up against the regression and every term has a partner.

The simulated throw starts at ground level with an upward speed of twenty metres per second.

height h(t) = 20t − 4.9t²
gravity magnitude g = −2 × coefficient on t²
                    = −2 × (−4.9) = 9.8 m/s²

The coefficient of squared time is minus half the gravitational acceleration. Its value is not the acceleration itself.

Take the noiseless case first. Five exact readings of the throw at t = 0, 1, 2, 3 and 4 are the “ideal case” button below. Press it, set the degree to 2, and read the chips.

Do not jump from a straight line to a degree-9 curve. The useful comparison is the three smallest models on the same five ideal points.

The straight line barely improves on the mean, because the best line through a symmetric arc is almost flat. The parabola explains everything. The difference between a useless model and a perfect one was one manufactured column, and the lesson is that a model can fail because the inputs it was given cannot express the pattern in the data.

The degree-1 model is not just worse by the numbers. Its misses have a shape. Plot each residual against time and the line’s residuals are negative at both ends and positive in the middle, which is the arc it could not draw, left over in the errors.

On the ideal throw the parabola’s residuals are zero. On the noisy throw they are a scatter with no obvious curve, which is what residuals look like when the model has taken the structure and left the noise. A visible pattern in the residuals says the model has left systematic structure unexplained, and that is worth checking on any fit, not only polynomial ones.

This has a name. A model too rigid to represent the pattern in the data is underfitting. It shows large systematic residuals, a poor training score and poor predictions, and on the ball the straight line is the textbook case.

Each degree adds one power of t as a new column, and each power buys the curve the possibility of one more turn.

The word is permits. A higher degree allows a more complex shape and does not require it. A degree-4 curve may turn three times, and the fitted coefficients decide how many of those turns actually appear. On the ideal throw a degree-4 fit puts 0.00 on t³ and 0.00 on t⁴ and draws the same parabola as degree 2, because the extra turns would not help. Section 11’s sliders are the place to feel this. Every shape is available and the amounts decide which powers contribute to the curve up.

Real readings of a thrown ball carry error, a wobble of the hand or the ruler that says nothing about the next throw. The noise slider in section 13’s box adds exactly that. It takes fifteen exact readings of the same 20t − 4.9t² throw and shifts each by its own draw, scaled to the amount you choose, and the “show the true throw” box draws the arc the readings came from.

What you see with noise on is three things at once. There is the true height, dashed. There is the noise, the vertical gap from the dashed curve to each dot. And there is the observed height, the dot itself, which is all the model is ever shown. The fit sees the measured values and not the hidden relationship that produced them, and every question about overfitting comes from that fact.

Now sweep the degree upward on the fifteen noisy readings, refitting at each step and scoring on the same fifteen.

The bars can only rise or hold still, and there is a formal reason. The degree-3 model contains every degree-2 model, since it can set the coefficient on t³ to zero and reproduce the degree-2 fit exactly, and the same holds one degree up from any degree. So the best degree-3 fit is at least as good on the training rows as the best degree-2 fit, and it may use the new column to do better. Training R² for nested models can only rise or stay unchanged as terms are added.

Degree 2 reads 0.996 here and degree 9 reads 0.997, with 10 coefficients doing the work of 3. Every step past 2 bought almost nothing on the training rows, and the next Part is about what it cost.

The argument against the high-degree curve is usually that it looks implausible, and a stronger argument is available. Split the readings. Fit on some of them and keep the rest hidden, and let the hidden ones judge.

Below, four of the fifteen noisy readings are hidden from the fit and from you. Pick a degree, decide where you think the four missing dots fall, and then reveal them.

A model should be judged, at least in part, by observations it did not use to fit itself. The training score measures how well the curve fits where it was told to fit. The held-out score measures whether what it learned travels.

Trace both scores across every degree and the picture arrives in one chart.

Training R² climbs from 0.222 at degree 1 to 0.999 at degree 9 and never falls. The held-out score rises to 0.989 at degree 2, holds near there through degree 5, drops to 0.693 at degree 6, and then collapses, −104 at degree 8 and −5479 at degree 9. A curve that scores 0.999 on the eleven readings it was fitted to is wrong by street widths on the four it was not.

That is overfitting. Flexibility improves the reproduction of the training sample while harming predictions on new observations. The signs are very low training error, large and unstable coefficients, sensitivity to individual measurements and residuals on the training data that are tiny while errors elsewhere grow. At degree 9 the coefficient on t⁴ is −229.5 and on t⁵ is +169.5, two huge terms cancelling each other almost everywhere and not quite cancelling between the points.

Section 13 scored 1.000 and was right to. Section 20 is heading toward 1.000 and is wrong. The number is the same and the situations are not.

A score becomes concerning because of how it was achieved and whether the model generalises, and not merely because it is high. Fourteen noisy points and a degree-13 curve will score exactly 1 by threading every one of them, and that says nothing about the fifteenth.

Set the three fits from section 19 beside each other in your head, degree 1, degree 2 and degree 9, all on the same eleven training points and judged by the same four hidden ones. The first is too stiff to follow the arc and misses both groups. The second follows the arc and scores well on both. The third follows the eleven exactly and misses the four badly.

Too little flexibility leaves structure unfitted. Too much reacts to the particular sample it was handed, noise included. A useful model sits between those, and the way to find it is performance on data outside the fit, which is the question the evaluation page takes up properly. None of the three is the simple model or the advanced one. They differ in flexibility, and that is the only axis this comparison is about.

Everything so far can be written in one line if the data is arranged as a matrix. Start from the table anyone would write down for the ideal throw, and rearrange it in stages.

The result is a matrix X with one row per observation and one column per input, a coefficient vector β, and their product, which is the vector of predictions.

ŷ = Xβ

The leading column of ones is the intercept’s seat. If α is to live inside β, X has to carry a column that multiplies it, and multiplying a column of ones by α adds the same α to every prediction. Leave that column out and the normal equations below solve for a model with no intercept at all.

Choose one row of X in the third stage above, say row 2, and it is the vector (1, 1, 1). Dot it with (α, β₁, β₂) and you get

ĥ₂ = α + β₁·1 + β₂·1 = 0 + 20 − 4.9 = 15.1

This is section 2’s weighted sum for one observation. Row 4 contains the intercept feature, time and squared time.

row 4 prediction = (1, 3, 9) · (0, 20, −4.9)
                 = 0 + 60 − 44.1
                 = 15.9

Matrix multiplication performs this dot product for every row at once. The result is the whole column of predictions.

The derivation is the regression page’s, in vector form. The residuals are a vector, RSS is that vector dotted with itself, and the bowl’s floor is where every partial derivative is zero at once.

Keeping track of shapes is most of the work of reading this, so here they are for n observations and p columns of X, where p counts the column of ones.

The normal equations describe the condition the least-squares solution satisfies. They do not say to compute (XᵀX)⁻¹Xᵀy by forming the inverse, and in practice nobody should. A working implementation solves the system directly, and numerical libraries usually go further and avoid forming XᵀX at all, using a QR decomposition or a singular value decomposition of X, which are better behaved when columns are nearly dependent. A unique solution also needs enough independent information in the columns, which is the next section.

Suppose one column holds t and another holds 2t. The second says nothing the first did not, and the fit has no way to decide how much of the effect belongs to each.

Every position of the slider is a different pair of coefficients and the identical set of predictions, so the residuals are identical too and RSS gives the fit no reason to prefer any of them. The bowl has a flat valley floor instead of a single lowest point, XᵀX has no inverse, and the normal equations have infinitely many solutions.

The name for this is collinearity, and for an exact copy like 2t it is perfect multicollinearity. A coefficient cannot be uniquely identified when its column can be reconstructed from the others. The fit behind this page refuses such a table outright rather than picking one of the infinitely many answers, and the nearly-dependent case, where two columns are almost but not quite copies, is one of the things ridge regression exists to handle.

High-degree polynomials behave dramatically outside the interval they were fitted on. Section 13’s box has a range control for this. Load the full throw, set the degree to 9, and switch from the observed range to a little beyond and then far beyond. The training points stay where they are, and the curve leaves the frame.

Even the parabola’s answers are wrong, since a real ball stops at the ground rather than continuing below it, and the degree-9 curve is wrong by a kilometre within two seconds of the last reading. A curve that fits well within the measured interval may be unreliable the moment it leaves it, and nothing in the training score can warn you.

A quieter problem comes from the columns themselves. The powers of t live on very different scales.

Wide differences in scale make XᵀX badly conditioned, so its solution swings on the last digits of the data. Centring and scaling the input before raising it to powers is the usual remedy, and the feature scaling page takes that up in full.

This page ends on two questions it has raised and not answered, and they are different questions.

How do we decide which degree predicts new data best? Section 19 held back four readings once. Doing that systematically, with every reading taking a turn as the judge, is cross-validation, and choosing a degree by it is what the held-out evaluation and grid search pages do.

How do we restrain a flexible model without throwing away the inputs it needs? Keep the degree-9 columns and penalise large coefficients, so the t⁴ and t⁵ terms cannot reach −229 and +169 in the first place. That is ridge and lasso, and it also answers section 26’s nearly-dependent columns.