Multiple Linear Regression | Machine Learning | Codeayan Courses

Four columns, and every number shifts

Multiple linear regression is the same model with more inputs: instead of one slope you fit one per column, and predictions come from adding up all the contributions plus an intercept. The arithmetic is a small step from linear regression with a single feature. The interpretation is not.

Here is the result that catches people. Fit rent against bedrooms alone and each bedroom is worth 7,438 rupees. Add floor area to the model and the same bedroom coefficient drops to 2,366. Nothing about the flats changed. The question the coefficient answers changed.

That shift is the subject of this chapter, along with what happens when two columns carry the same information and why your training score can never get worse when you add a column.

Solving for all the coefficients at once

With k features the model is y = b + w₁x₁ + w₂x₂ + ... + wₖxₖ. Each w is the change in the prediction for a one-unit change in that column, holding the others fixed. The closed-form solution is the same least-squares problem written in matrix form, using the operations from linear algebra for machine learning.

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(14)
N = 900
area = rng.normal(1100, 260, N).clip(400, 2400)
bedrooms = np.clip((area / 420 + rng.normal(0, 0.5, N)).round(), 1, 5)
age = rng.integers(0, 35, N).astype(float)
metro_km = rng.gamma(2.0, 1.4, N).clip(0.2, 12)

rent = (6000 + 21.5 * area + 2400 * bedrooms - 180 * age - 1500 * metro_km
        + rng.normal(0, 3000, N))

flats = pd.DataFrame({"area_sqft": area, "bedrooms": bedrooms,
                      "age_years": age, "metro_km": metro_km})

X_with_intercept = np.hstack([np.ones((N, 1)), flats.to_numpy()])
theta = np.linalg.lstsq(X_with_intercept, rent, rcond=None)[0]
print(theta.round(2))
# [ 4941.22    22.66  2366.36  -178.17 -1521.02]

model = LinearRegression().fit(flats, rent)
print(round(float(model.intercept_), 2), model.coef_.round(2))
# 4941.22 [   22.66  2366.36  -178.17 -1521.02]

Identical, because LinearRegression runs the same least-squares solve. Read the coefficients in order: 22.66 rupees per square foot, 2,366 per bedroom, minus 178 per year of age, minus 1,521 per kilometre from the nearest metro station.

Use np.linalg.lstsq rather than inverting anything by hand. It handles the case where columns are linearly dependent, which np.linalg.solve cannot.

What multiple linear regression does to your coefficients

solo = LinearRegression().fit(flats[["bedrooms"]], rent)
print("bedrooms alone:", round(float(solo.coef_[0]), 1))          # 7437.6
print("bedrooms with area:", round(float(model.coef_[1]), 1))     # 2366.4
print("correlation:", round(float(np.corrcoef(area, bedrooms)[0, 1]), 3))  # 0.746

Bigger flats have more bedrooms, with a correlation of 0.746. On its own, the bedroom coefficient absorbs the value of all that extra floor area. Once area is in the model, the bedroom coefficient answers a narrower question: what is an extra bedroom worth between two flats of the same size, same age and same distance from the metro? That is 2,366 rupees, and it is the number a developer deciding how to partition a fixed floor plate actually needs.

Neither value is wrong. They answer different questions, and the phrase “holding the others constant” is the whole difference. Quote a coefficient without naming the rest of the model and you have said nothing checkable.

One consequence worth keeping: adding or removing a column changes every other coefficient. There is no stable “effect of bedrooms” independent of specification.

When two columns say the same thing

Multicollinearity is what happens when one feature can be predicted from the others. The extreme case is easy to construct.

flats2 = flats.copy()
flats2["area_sqm"] = flats2["area_sqft"] * 0.0929      # the same column, rescaled

for seed in [0, 1, 2]:
    idx = np.random.default_rng(seed).choice(N, N, replace=True)
    m = LinearRegression().fit(flats2.iloc[idx], rent[idx])
    print(f"seed {seed}: area_sqft {m.coef_[0]:7.2f}   area_sqm {m.coef_[4]:7.2f}")
# seed 0: area_sqft   22.18   area_sqm    2.06
# seed 1: area_sqft   22.95   area_sqm    2.13
# seed 2: area_sqft   22.17   area_sqm    2.06

scikit-learn does not explode, which is worth stating because textbooks often imply it will. The least-squares solver returns a minimum-norm solution, so it splits the area effect across both columns rather than producing wild numbers.

The damage is to interpretation. Neither 22.18 nor 2.06 is the value of a square foot. The real effect is the two combined, and no single coefficient in that table means anything on its own. Had you used np.linalg.solve on the normal equations instead, the matrix would have been singular and the call would have failed outright.

The standard diagnostic is the variance inflation factor: regress each column on all the others and see how predictable it is.

def vif(frame, col):
    others = [c for c in frame.columns if c != col]
    r2 = LinearRegression().fit(frame[others], frame[col]).score(frame[others], frame[col])
    return 1 / (1 - r2) if r2 < 1 else np.inf

for c in flats2.columns:
    print(f"{c:11s} {vif(flats2, c):.1f}")
# area_sqft   inf
# bedrooms    2.3
# age_years   1.0
# metro_km    1.0
# area_sqm    inf

Infinite for the duplicated pair, around 1 for the independent columns. A rough working rule: above 5 is worth investigating, above 10 means the coefficients are not individually interpretable.

The fix is usually to drop one of the pair. When the columns are genuinely correlated rather than duplicated, and you still need them all, a penalty on coefficient size stabilises the fit, which is what regularised regression is for. Note that multicollinearity hurts interpretation far more than prediction: a model with two copies of the same column still predicts fine.

Why your training score never falls

noise = pd.DataFrame(rng.normal(0, 1, (N, 10)),
                     columns=[f"junk_{i}" for i in range(10)])
wider = pd.concat([flats.reset_index(drop=True), noise], axis=1)

print(round(LinearRegression().fit(flats, rent).score(flats, rent), 5))    # 0.87892
print(round(LinearRegression().fit(wider, rent).score(wider, rent), 5))    # 0.88065

Ten columns of pure random noise raised the training score. They always do, because least squares can only improve or tie when given another column to use, and with enough junk columns it will fit the noise exactly.

Which is why a training score is not evidence of anything. Judge on held-out data, and if you need an in-sample number, use an adjusted version that charges for each column added. Categorical inputs go through encoding first, and note that one-hot encoding creates exactly the collinearity discussed above unless you drop one level. Solver details are documented under numpy.linalg.lstsq.

Frequently Asked Questions

What is the difference between simple and multiple linear regression?

Simple regression uses one input and fits one slope. Multiple regression uses several and fits one coefficient each, with every coefficient meaning “the effect of this column holding the others fixed”. That conditioning is why a coefficient changes value when you add or remove other columns.

What is multicollinearity and why does it matter?

It occurs when one feature can be predicted from the others, which splits a shared effect across several coefficients unpredictably. Prediction quality barely suffers; interpretation breaks completely. Check with variance inflation factors and treat anything above 10 as not individually interpretable.

Why does adding features always increase R squared?

Least squares can set a useless column’s coefficient to zero, so adding one can never make the training fit worse, and random noise usually improves it slightly. Ten junk columns raised training R squared here. Use held-out data or an adjusted measure that penalises column count.

Key Takeaways

State the full model whenever you quote a coefficient, because each one means “holding the other columns fixed” and changes when the specification changes.
Compute variance inflation factors before interpreting anything, and treat values above 10 as a sign the coefficients cannot be read individually.
Expect scikit-learn to return a stable minimum-norm answer under collinearity rather than failing, which means the problem is silent.
Drop one column from any duplicated pair, and reach for a coefficient penalty when correlated columns must all stay in the model.
Ignore training R squared entirely when comparing models of different sizes, since ten columns of noise raised it here.