Linear and Logistic Regression

Article · 11:00

The models to try first: fast, interpretable, hard to overfit when regularised, and often within a point or two of anything more elaborate. Also the only ones whose decisions you can explain in a sentence.

This lesson has 10 runnable examples. Press Try yourself! under any of them to open it in the editor.

Linear and logistic regression are the models to try first and the ones to beat. They fit in seconds, their coefficients are readable, and a problem where they do badly tells you something before you have spent a day on anything larger.

Linear regression

The model is a weighted sum: prediction = intercept + w1 * x1 + w2 * x2 + ..., with the weights chosen to minimise the squared error.

PY Example 1 30 lines

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, r2_score
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(3)
n = 2000

area = rng.normal(95, 25, n).clip(30)
bedrooms = np.clip((area / 30 + rng.normal(0, 0.5, n)).round(), 1, 6)
age = rng.gamma(2, 12, n).clip(0, 120)

# A genuinely linear relationship, plus noise.
price = (42_000 + 2_150 * area + 9_500 * bedrooms - 380 * age
         + rng.normal(0, 22_000, n))

frame = pd.DataFrame({'area': area, 'bedrooms': bedrooms, 'age': age, 'price': price})
X_tr, X_te, y_tr, y_te = train_test_split(
    frame[['area', 'bedrooms', 'age']], frame['price'], test_size=0.3, random_state=0)

model = LinearRegression().fit(X_tr, y_tr)
predicted = model.predict(X_te)

print(pd.DataFrame({'feature': X_tr.columns, 'coefficient': model.coef_.round(0),
                    'true': [2150, 9500, -380]}).to_string(index=False))
print(f'intercept {model.intercept_:,.0f} (true 42,000)')
print()
print(f'R2  {r2_score(y_te, predicted):.3f}')
print(f'MAE {mean_absolute_error(y_te, predicted):,.0f}')

Open this example in the editor and change it. Nothing you type is saved.

Each coefficient reads as: holding the other features constant, one more unit of this feature changes the prediction by this much. That sentence is the whole reason to start here - no other model family gives it to you for free.

R2 is the share of the variance the model accounts for. It is comparable within one dataset and meaningless across datasets, because it depends on how much variance there was to begin with.

The assumptions, and what breaks them

PY Example 2 24 lines

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

rng = np.random.default_rng(7)
n = 1500
x = rng.uniform(1, 20, n)

shapes = {
    'linear': 3 * x + rng.normal(0, 3, n),
    'monotonic curve': 0.4 * x ** 2 + rng.normal(0, 3, n),
    'u-shaped': 0.5 * (x - 10) ** 2 + rng.normal(0, 3, n),
    'periodic': 20 * np.sin(x) + rng.normal(0, 3, n),
    'step at the middle': np.where(x > 10, 40, 10) + rng.normal(0, 3, n),
}
for name, y in shapes.items():
    model = LinearRegression().fit(x.reshape(-1, 1), y)
    print(f'{name:<20} R2 {r2_score(y, model.predict(x.reshape(-1, 1))):.3f}')
print()
print('A monotonic curve is approximated passably by a straight line;')
print('a U-shape and a periodic signal are not approximated at all,')
print('because no single slope describes them. The fix is usually a')
print('transformed feature, not a different algorithm.')

Open this example in the editor and change it. Nothing you type is saved.

The fix for a curve is often a feature, not a model:

PY Example 3 17 lines

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score

rng = np.random.default_rng(7)
n = 1500
x = rng.uniform(1, 20, n)
y = 0.5 * (x - 10) ** 2 + rng.normal(0, 3, n)       # the U-shape

plain = LinearRegression().fit(x.reshape(-1, 1), y)
with_square = LinearRegression().fit(np.column_stack([x, x ** 2]), y)

print(f'x only        R2 {r2_score(y, plain.predict(x.reshape(-1, 1))):.3f}')
print(f'x and x**2    R2 {r2_score(y, with_square.predict(np.column_stack([x, x ** 2]))):.3f}')
print()
print('"Linear" means linear in the PARAMETERS, not in the inputs.')
print('A squared term, a log, an interaction - all still linear models.')

Open this example in the editor and change it. Nothing you type is saved.

Two more assumptions worth checking:

Constant error variance. If the errors grow with the prediction, the model is overconfident at the top end. The usual fix is to model log(y) instead of y.

PY Example 4 20 lines

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_percentage_error
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(11)
n = 2000
x = rng.uniform(1, 10, n)
# Error proportional to the level - very common with money.
y = 1000 * np.exp(0.4 * x) * rng.lognormal(0, 0.25, n)

X_tr, X_te, y_tr, y_te = train_test_split(x.reshape(-1, 1), y, test_size=0.3, random_state=0)

direct = LinearRegression().fit(X_tr, y_tr)
logged = LinearRegression().fit(X_tr, np.log(y_tr))

print(f'on y:        MAPE {mean_absolute_percentage_error(y_te, direct.predict(X_te)):.1%}')
print(f'on log(y):   MAPE '
      f'{mean_absolute_percentage_error(y_te, np.exp(logged.predict(X_te))):.1%}')

Open this example in the editor and change it. Nothing you type is saved.

No extreme multicollinearity. Two features that carry the same information make the coefficients unstable, even though the predictions stay fine.

PY Example 5 19 lines

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

rng = np.random.default_rng(13)
n = 500
a = rng.normal(0, 1, n)
b = a + rng.normal(0, 0.02, n)           # almost identical to a
y = 3 * a + rng.normal(0, 0.5, n)

for seed in range(3):
    sample = rng.choice(n, n, replace=True)
    model = LinearRegression().fit(np.column_stack([a[sample], b[sample]]), y[sample])
    print(f'sample {seed}: coefficients {model.coef_.round(2)}  '
          f'(their sum {model.coef_.sum():.2f})')
print()
print('The individual coefficients swing wildly between samples; their')
print('SUM is stable. Do not interpret a coefficient when two features')
print('are nearly the same thing.')

Open this example in the editor and change it. Nothing you type is saved.

Logistic regression

For a yes/no label, the same weighted sum is passed through a function that squashes it into a probability.

PY Example 6 32 lines

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, accuracy_score
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

rng = np.random.default_rng(17)
n = 3000

tenure = rng.integers(1, 60, n)
tickets = rng.poisson(1.6, n)
spend = rng.gamma(3, 25, n)

logit = -0.5 - 0.045 * tenure + 0.42 * tickets - 0.011 * spend
probability = 1 / (1 + np.exp(-logit))
churned = (rng.random(n) < probability).astype(int)

X = np.column_stack([tenure, tickets, spend])
X_tr, X_te, y_tr, y_te = train_test_split(X, churned, test_size=0.3,
                                          random_state=0, stratify=churned)

model = Pipeline([('scale', StandardScaler()),
                  ('model', LogisticRegression(max_iter=1000))]).fit(X_tr, y_tr)

probabilities = model.predict_proba(X_te)[:, 1]
print(f'positive rate {churned.mean():.1%}')
print(f'accuracy {accuracy_score(y_te, model.predict(X_te)):.3f}')
print(f'ROC AUC  {roc_auc_score(y_te, probabilities):.3f}')
print()
print('first five probabilities:', probabilities[:5].round(3))

Open this example in the editor and change it. Nothing you type is saved.

The output is a probability, which is more useful than the class. predict() applies a 0.5 threshold, which is a choice - and rarely the right one, as the evaluation lesson shows.

Coefficients in logistic regression are log-odds, and exponentiating them gives an odds ratio:

PY Example 7 30 lines

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

rng = np.random.default_rng(19)
n = 4000
tenure = rng.integers(1, 60, n)
tickets = rng.poisson(1.6, n)
spend = rng.gamma(3, 25, n)
logit = -0.5 - 0.045 * tenure + 0.42 * tickets - 0.011 * spend
churned = (rng.random(n) < 1 / (1 + np.exp(-logit))).astype(int)

# Unscaled, so the coefficients are in the features' own units.
model = LogisticRegression(max_iter=2000).fit(
    np.column_stack([tenure, tickets, spend]), churned)

table = pd.DataFrame({
    'feature': ['tenure_months', 'tickets_90d', 'monthly_spend'],
    'coefficient': model.coef_[0].round(4),
    'true': [-0.045, 0.42, -0.011],
})
table['odds_ratio'] = np.exp(model.coef_[0]).round(3)
print(table.to_string(index=False))
print()
print('Read the third row as: each extra ticket multiplies the odds of')
print(f'churn by {np.exp(model.coef_[0][1]):.2f} - about a '
      f'{(np.exp(model.coef_[0][1]) - 1) * 100:.0f}% increase in the odds,')
print('holding tenure and spend constant.')

Open this example in the editor and change it. Nothing you type is saved.

Why start here

PY Example 8 33 lines

import numpy as np
import pandas as pd
import time
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(23)
n = 4000
X = rng.normal(size=(n, 8))
# A LINEAR relationship, which is the common case in tabular business data.
logit = 1.2 * X[:, 0] - 0.8 * X[:, 1] + 0.5 * X[:, 2] + rng.normal(0, 0.8, n)
y = (logit > 0).astype(int)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)

models = {
    'logistic regression': Pipeline([('s', StandardScaler()),
                                     ('m', LogisticRegression(max_iter=1000))]),
    'random forest': RandomForestClassifier(n_estimators=300, random_state=0),
    'gradient boosting': GradientBoostingClassifier(random_state=0),
}
for name, model in models.items():
    start = time.perf_counter()
    model.fit(X_tr, y_tr)
    seconds = time.perf_counter() - start
    auc = roc_auc_score(y_te, model.predict_proba(X_te)[:, 1])
    print(f'{name:<22} AUC {auc:.3f}   fitted in {seconds * 1000:7.0f}ms')
print()
print('On a linear problem the simplest model wins and is far faster.')
print('That is not a coincidence - it is the right inductive bias.')

Open this example in the editor and change it. Nothing you type is saved.

Four practical reasons it earns the first attempt:

It is fast, so you learn whether the problem is learnable in minutes. It is readable, so a coefficient with the wrong sign reveals a data problem immediately. It extrapolates, unlike tree models, which matters when inputs move outside the training range. It is a baseline that is hard to argue with. If the complicated model does not beat it, ship the simple one.

Regularisation: ridge and lasso

PY Example 9 37 lines

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.metrics import mean_absolute_error
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

rng = np.random.default_rng(29)
n, p = 300, 60

X = rng.normal(size=(n, p))
true_weights = np.zeros(p)
true_weights[:5] = [3.0, -2.0, 1.5, -1.0, 0.8]      # only 5 of 60 matter
y = X @ true_weights + rng.normal(0, 1.0, n)

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0)

rows = []
for name, estimator in [('linear', LinearRegression()),
                        ('ridge (a=1)', Ridge(alpha=1.0)),
                        ('ridge (a=10)', Ridge(alpha=10.0)),
                        ('lasso (a=0.1)', Lasso(alpha=0.1, max_iter=5000)),
                        ('lasso (a=0.3)', Lasso(alpha=0.3, max_iter=5000))]:
    model = Pipeline([('scale', StandardScaler()), ('model', estimator)]).fit(X_tr, y_tr)
    coefficients = model.named_steps['model'].coef_
    rows.append({
        'model': name,
        'test MAE': round(mean_absolute_error(y_te, model.predict(X_te)), 3),
        'nonzero coefs': int((np.abs(coefficients) > 1e-6).sum()),
        'found the 5': int((np.abs(coefficients[:5]) > 1e-6).sum()),
    })
print(pd.DataFrame(rows).to_string(index=False))
print()
print('With 60 features and 210 training rows, plain linear regression')
print('fits the noise. Lasso sets most coefficients to exactly zero and')
print('keeps the five that matter.')

Open this example in the editor and change it. Nothing you type is saved.

Ridge shrinks every coefficient towards zero without reaching it, which stabilises them when features are correlated. Lasso drives some to exactly zero, which performs feature selection. ElasticNet does both.

Both require scaled features - the penalty is on the coefficient size, so a feature measured in pounds is penalised a thousand times more than the same feature in thousands of pounds.

What they cannot do

PY Example 10 31 lines

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(31)
n = 3000
X = rng.normal(size=(n, 2))

# The label depends on an INTERACTION: positive when the signs agree.
y = ((X[:, 0] > 0) == (X[:, 1] > 0)).astype(int)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0)

plain = LogisticRegression(max_iter=1000).fit(X_tr, y_tr)
print(f'logistic, raw features:   AUC '
      f'{roc_auc_score(y_te, plain.predict_proba(X_te)[:, 1]):.3f}')

# Give it the interaction explicitly and it recovers.
with_interaction = LogisticRegression(max_iter=1000).fit(
    np.column_stack([X_tr, X_tr[:, 0] * X_tr[:, 1]]), y_tr)
print(f'logistic, with x0 * x1:   AUC '
      f'{roc_auc_score(y_te, with_interaction.predict_proba(np.column_stack([X_te, X_te[:, 0] * X_te[:, 1]]))[:, 1]):.3f}')

forest = RandomForestClassifier(n_estimators=200, random_state=0).fit(X_tr, y_tr)
print(f'random forest:            AUC '
      f'{roc_auc_score(y_te, forest.predict_proba(X_te)[:, 1]):.3f}')
print()
print('A linear model cannot find an interaction unless you give it')
print('one. A tree finds it automatically - which is the main reason')
print('trees are the next thing to try.')

Open this example in the editor and change it. Nothing you type is saved.

The three limitations, in order of how often they bite:

No interactions unless you write them. x1 * x2 must be an explicit feature.
No non-linearity unless you transform. log(x) and x**2 must be explicit.
Sensitive to outliers, because squared error punishes a single large residual enormously. HuberRegressor is the robust alternative.

When it goes wrong

Symptom	Cause
A coefficient has the wrong sign	Collinearity, or a confounder in the features
Coefficients swing between samples	Two features carry the same information
Good R2, obviously curved residuals	Linear fit on a non-linear relationship
The solver warns about convergence	Unscaled features; add a `StandardScaler`
Lasso zeroed a feature you believe in	It is correlated with another one it kept
Errors grow with the prediction	Model `log(y)` instead of `y`
The model misses an obvious interaction	Add the product term, or use a tree
Regularisation hurt badly	Features are unscaled; the penalty is uneven

A check you can run

Fit a linear or logistic regression before the model you intended to build, on the same split, and write both scores down.

If the gap is small, ship the simple one: it is faster, explainable, and it extrapolates. If the gap is large, you have learned something specific - the relationship is non-linear or has interactions - and that tells you which features to engineer rather than which library to install.