Feature Engineering · Notebook 3 companion articleFeature Transformation for Numeric Variables

Start hereIntroduction#

The first two notebooks fixed broken things: a customer ID that shouldn't have been a feature, a blank cell hiding as text, a missing age. This one is different. Every number here is correct. Nothing is mistyped or missing. Yet some of these correct numbers will make a model do worse than if it ignored the data completely.

That's because a model sees your numbers through its own assumptions, and those assumptions are rarely stated. K-nearest neighbours assumes distance means the same thing in every column. A linear model assumes a straight line fits the relationship. Squaring a column assumes it has no extreme values, because squaring makes them even more extreme. When these assumptions are broken, nothing errors. The model just gets worse.

Notebook 2 already gave two previews of this. Its KNN imputer had to scale the columns first, or Fare would have controlled every distance. And log1p brought Fare's skew down from 4.79 to 0.39 without deleting a row. This notebook explains both properly.

It closes the gap between what your numbers look like and what your model needs, using four kinds of transformation:

• scaling, so no column dominates just because it uses bigger numbers
• non-linear transforms, to fix columns that are heavily skewed
• binning, which turns a number into ranges, and when that's worth doing
• polynomial features, which let a model follow a curve

The two datasets#

Telco comes back from the first two notebooks for the classification work. This notebook also adds California Housing, a 20,640-row census dataset that comes with scikit-learn. It gives us a regression target and, more importantly, columns that are badly skewed, which Telco's rarely are.

Two things to know about it. Each row is a census block group, a neighborhood of roughly a few hundred to a few thousand people, so columns like AveRooms are averages, not one house. And some columns were capped when the data was made: median income stops at 15 (in tens of thousands of dollars), house age at 52 years, and the target, median house value in hundreds of thousands, at 5, or $500,000.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression, Ridge
from sklearn.model_selection import KFold, StratifiedKFold, cross_val_score, train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import (KBinsDiscretizer, MinMaxScaler, PolynomialFeatures,
                                   PowerTransformer, RobustScaler, StandardScaler)

plt.style.use('seaborn-v0_8-whitegrid')
pd.set_option('display.max_columns', 30)
pd.set_option('display.width', 170)

RANDOM_STATE = 42

# Telco, cleaned exactly as notebook 2 left it
telco = pd.read_csv('datasets/telco_churn.csv')
telco['TotalCharges'] = pd.to_numeric(telco['TotalCharges'], errors='coerce')
telco = telco[telco['tenure'] > 0].copy()
churned = (telco['Churn'] == 'Yes').astype(int)

# California housing: 20,640 census block groups, target is median house value
ca_raw = fetch_california_housing()
housing = pd.DataFrame(ca_raw.data, columns=ca_raw.feature_names)
house_value = pd.Series(ca_raw.target, name='house_value')

print(f'Telco    : {telco.shape[0]:,} rows, target = churn (classification)')
print(f'Housing  : {housing.shape[0]:,} rows, target = median house value (regression)')
print()
print('Housing columns, sorted by how skewed they are:')
print(housing.skew().sort_values(ascending=False).round(2).to_string())

Telco    : 7,032 rows, target = churn (classification)
Housing  : 20,640 rows, target = median house value (regression)

Housing columns, sorted by how skewed they are:
AveOccup      97.64
AveBedrms     31.32
AveRooms      20.70
Population     4.94
MedInc         1.65
Latitude       0.47
HouseAge       0.06
Longitude     -0.30

That skew column sums up the whole notebook. HouseAge, at 0.06, is almost perfectly balanced and needs nothing. AveOccup, at 97.6, is so lopsided that a handful of extreme blocks control its standard deviation, as section 2 shows. Same table, same census, and two columns that need completely different treatment. Applying one transformation to every numeric column out of habit is exactly the mistake to avoid.

Section 1Why scaling matters#

A real, measurable way that useful features get ignored

The usual explanation is that "some algorithms are sensitive to scale." That's true but not very helpful, because it doesn't tell you when it matters or how much. Telco's numeric columns show the problem clearly enough to measure.

Distance-based algorithms, especially k-nearest neighbours, decide which rows are similar by measuring how far apart they are. Euclidean distance takes the difference in each column, squares it, and adds them up. Because of the squaring, a column with big numbers adds a lot to the total and a column with small numbers adds very little. That happens no matter which column is actually more useful.

# Defined once here and reused for the rest of the notebook
telco_numeric = telco[['tenure', 'MonthlyCharges', 'TotalCharges']]

spread = pd.DataFrame({
    'min': telco_numeric.min(),
    'max': telco_numeric.max(),
    'range': telco_numeric.max() - telco_numeric.min(),
    'std': telco_numeric.std(),
    'variance': telco_numeric.var(),
})
spread['share_of_total_variance'] = (spread['variance'] / spread['variance'].sum())

print(spread.round(4).to_string())

                  min      max   range        std      variance  share_of_total_variance
tenure           1.00    72.00    71.0    24.5453  6.024698e+02                   0.0001
MonthlyCharges  18.25   118.75   100.5    30.0860  9.051658e+02                   0.0002
TotalCharges    18.80  8684.80  8666.0  2266.7714  5.138252e+06                   0.9997

TotalCharges makes up 99.97% of the combined variance across the three columns. So in a distance calculation, it decides the answer almost by itself. tenure and MonthlyCharges each hold one or two ten-thousandths of the total, so they barely affect which neighbours get picked. (Variance is the right thing to measure here: on average, the squared gap between two random rows in a column is exactly twice its variance.)

To see what this costs, the notebook adds four yes/no features that we know predict churn: month-to-month contract, fibre internet, paying by electronic check, and no tech support. Each is a plain 0 or 1.

features = pd.DataFrame({
    'tenure': telco['tenure'],
    'MonthlyCharges': telco['MonthlyCharges'],
    'TotalCharges': telco['TotalCharges'],
    'is_month_to_month': (telco['Contract'] == 'Month-to-month').astype(int),
    'has_fiber': (telco['InternetService'] == 'Fiber optic').astype(int),
    'electronic_check': (telco['PaymentMethod'] == 'Electronic check').astype(int),
    'no_tech_support': (telco['TechSupport'] == 'No').astype(int),
})

variance_share = features.var() / features.var().sum()
print('Share of total variance held by each feature:')
print(variance_share.round(6).to_string())

Share of total variance held by each feature:
tenure               0.000117
MonthlyCharges       0.000176
TotalCharges         0.999706
is_month_to_month    0.000000
has_fiber            0.000000
electronic_check     0.000000
no_tech_support      0.000000

All four flags round to a variance share of 0.000000, at six decimal places. A column that's only ever 0 or 1 is tiny next to a column that runs into the thousands.

That's the actual math. Two customers 8,100 dollars apart on TotalCharges add 65.6 million to the squared distance. Two customers on opposite sides of the month-to-month flag add just 1. You could give this model the most predictive yes/no feature in the dataset, and it would barely change a prediction.

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)

def churn_auc(model, X=features):
    return cross_val_score(model, X, churned, cv=cv, scoring='roc_auc').mean()


knn = lambda: KNeighborsClassifier(n_neighbors=15)
flags = features[['is_month_to_month', 'has_fiber', 'electronic_check', 'no_tech_support']]

print('UNSCALED, does adding the four flags help?')
print(f'  3 numeric columns only   : AUC {churn_auc(knn(), telco_numeric):.4f}')
print(f'  + the 4 binary flags     : AUC {churn_auc(knn(), features):.4f}')
print()
print('SCALED, the same seven columns:')
print(f'  StandardScaler + KNN     : AUC {churn_auc(make_pipeline(StandardScaler(), knn()), features):.4f}')
print()
print('Are those flags even informative? Test them on their own, scaled:')
print(f'  the 4 flags ALONE        : AUC {churn_auc(make_pipeline(StandardScaler(), knn()), flags):.4f}')

UNSCALED, does adding the four flags help?
  3 numeric columns only   : AUC 0.7846
  + the 4 binary flags     : AUC 0.7848

SCALED, the same seven columns:
  StandardScaler + KNN     : AUC 0.8272

Are those flags even informative? Test them on their own, scaled:
  the 4 flags ALONE        : AUC 0.8022

Read those four numbers together. On their own, with scaling, the four flags score an AUC of 0.8022, higher than the 0.7846 the three numeric columns score unscaled. Add the same flags to an unscaled KNN, and AUC goes from 0.7846 to 0.7848, which is just noise. Scale the same seven columns and AUC jumps to 0.8272. The data didn't change between those last two runs. Only the units did.

Section 2The three scalers, and what each does to an outlier#

All three put columns on a similar scale, but they define "similar" differently

Three formulas, each making columns comparable in a different way:

StandardScaler subtracts the mean and divides by the standard deviation, so every column ends up with a mean of 0 and a standard deviation of 1.

z = (x − μ) / σ

MinMaxScaler squeezes a column into a fixed range, 0 to 1 by default.

x' = (x − x_min) / (x_max − x_min)

RobustScaler subtracts the median and divides by the interquartile range.

x' = (x − median) / (Q₃ − Q₁)

The formulas already tell you how each one handles extreme values. Notebook 2 showed that the mean and standard deviation get pulled around by a few huge values, when three Fare values moved the z-score cutoff by 27 pounds. StandardScaler uses the mean and standard deviation. MinMaxScaler uses the two most extreme values, which makes it the most fragile. RobustScaler uses the median and quartiles, which barely move no matter how extreme a few outliers are.

California Housing's AveOccup column (average people per household in each census block) is a perfect test. The median is a normal 2.82 people, but the maximum is 1,243.33, most likely a place like a dormitory or institution, where many people are counted but very few households.

occupancy = housing[['AveOccup']]
print(f'AveOccup: min {occupancy.values.min():.2f}, median {np.median(occupancy):.2f}, max {occupancy.values.max():.2f}')
print()

rows = []
for name, scaler in [('StandardScaler', StandardScaler()),
                     ('MinMaxScaler', MinMaxScaler()),
                     ('RobustScaler', RobustScaler())]:
    scaled = scaler.fit_transform(occupancy).ravel()
    q1, q3 = np.percentile(scaled, [25, 75])
    lo, hi = np.percentile(scaled, [0.5, 99.5])
    rows.append({
        'scaler': name,
        'min': scaled.min(),
        'median': np.median(scaled),
        'max': scaled.max(),
        'IQR_width': q3 - q1,
        'middle_99pct_spans': f'[{lo:.3f}, {hi:.3f}]',
    })

print(pd.DataFrame(rows).round(4).to_string(index=False))

AveOccup: min 0.69, median 2.82, max 1243.33

        scaler     min  median       max  IQR_width middle_99pct_spans
StandardScaler -0.2290 -0.0243  119.4191     0.0821    [-0.156, 0.294]
  MinMaxScaler  0.0000  0.0017    1.0000     0.0007     [0.001, 0.004]
  RobustScaler -2.4936  0.0000 1455.1161     1.0000    [-1.608, 3.876]

Focus on the IQR_width column. It shows how much room the normal neighborhoods have after scaling.

MinMaxScaler fails badly here. The middle 99% of neighborhoods gets squeezed into the range 0.001 to 0.004. The middle half of neighborhoods differ by almost a whole person per household, yet they now fit in a range just 0.0007 wide. To a distance-based model like KNN, they look almost identical. One far-away outlier sets the top of the scale for everyone.

StandardScaler does okay, but the outlier still hurts. The middle 99% spans about 0.45 units, which is tighter than it should be, because the extreme value inflates the standard deviation everyone is divided by.

RobustScaler isn't affected. Its interquartile width is 1.0 by design, and the middle 99% spans about 5.5 units, which is plenty of room. The outlier still sits far out at 1,455, but it doesn't squash everyone else.

It's easy to assume scaling changes the shape of the data. Let's check: does it change the shape, or only the units?

fig, axes = plt.subplots(1, 4, figsize=(17, 3.8))

income = housing[['MedInc']]
panels = [('Original', income.values.ravel()),
          ('StandardScaler', StandardScaler().fit_transform(income).ravel()),
          ('MinMaxScaler', MinMaxScaler().fit_transform(income).ravel()),
          ('RobustScaler', RobustScaler().fit_transform(income).ravel())]

for ax, (title, values) in zip(axes, panels):
    ax.hist(values, bins=60, color='#3b76af', edgecolor='white')
    ax.set(title=f'{title}\nskew {pd.Series(values).skew():.3f}', xlabel='value')

axes[0].set_ylabel('block groups')
plt.tight_layout()
plt.show()

Every chart shows a skew of 1.647, the same as the original. Scaling only shifts and stretches the data evenly, and neither of those changes its shape. A long right tail stays exactly the same after all three scalers.

This is the most useful fact in this section, and the most often confused: scaling changes units, not shape. A skewed column stays just as skewed after scaling. Fixing shape is a different job, covered in section 5.

Section 3Which algorithms actually need scaling#

Many projects scale every numeric column by habit. Seeing where it does nothing shows what it's really for

Trees are the easiest case. A decision tree asks questions like "is TotalCharges above 1,400?" If you multiply every value by 1,000, subtract the mean or squeeze it into 0 to 1, the threshold number changes, but the same rows still land on each side. Rescaling can only relabel the threshold. It can't change which rows are above or below it.

forest = RandomForestClassifier(n_estimators=120, random_state=RANDOM_STATE, n_jobs=1)

forest_raw = churn_auc(forest)
forest_scaled = churn_auc(make_pipeline(StandardScaler(), forest))

print(f'RandomForest on raw features : AUC {forest_raw:.4f}')
print(f'RandomForest + StandardScaler: AUC {forest_scaled:.4f}')
print(f'Difference                   : {abs(forest_scaled - forest_raw):.4f}')

RandomForest on raw features : AUC 0.8012
RandomForest + StandardScaler: AUC 0.8020
Difference                   : 0.0008

A difference of 0.0008 is just noise. In perfect math the two numbers would be identical, because a tree only compares values to thresholds and scaling keeps the order the same. The tiny gap comes from computer rounding: rescaled values give slightly different split points, and a tie between two equally good splits can go a different way. With a different random seed, the gap could flip direction.

So for tree models, scaling does nothing. But it also costs nothing, so it's fine to include a scaler in your pipeline even if today's model is a tree. Tomorrow's model might not be.

raw_model = LogisticRegression(max_iter=5000).fit(telco_numeric, churned)
scaled_model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000)).fit(telco_numeric, churned)

raw_auc = cross_val_score(LogisticRegression(max_iter=5000), telco_numeric, churned, cv=cv, scoring='roc_auc').mean()
scaled_auc = cross_val_score(make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000)),
                             telco_numeric, churned, cv=cv, scoring='roc_auc').mean()

print(f'LogisticRegression raw    : AUC {raw_auc:.4f}   converged in {raw_model.n_iter_[0]} iterations')
print(f'LogisticRegression scaled : AUC {scaled_auc:.4f}   converged in {scaled_model.named_steps["logisticregression"].n_iter_[0]} iterations')
print()
print('Coefficients:')
print(pd.DataFrame({
    'raw': raw_model.coef_[0],
    'scaled': scaled_model.named_steps['logisticregression'].coef_[0],
}, index=telco_numeric.columns).round(5).to_string())

LogisticRegression raw    : AUC 0.8087   converged in 54 iterations
LogisticRegression scaled : AUC 0.8087   converged in 11 iterations

Coefficients:
                    raw   scaled
tenure         -0.06711 -1.61603
MonthlyCharges  0.03020  0.91280
TotalCharges    0.00015  0.29866

Here the AUC doesn't change at all: 0.8087 either way on these three columns. That doesn't mean scaling was pointless, because two other things changed that a single score won't show.

The model needed 11 steps to converge instead of 54. When one column's scale is much bigger than the others, the optimizer behind logistic regression has to zigzag through a long, narrow valley to find the answer. Scaling makes that valley round, so it gets there faster.

On three columns and 7,000 rows, a few extra steps don't matter. On a wide dataset with a slower solver, the model can hit its step limit before it finishes. Scikit-learn then shows a ConvergenceWarning, and the model you get is only partly trained.

The coefficients became comparable. In the unscaled model, TotalCharges has a coefficient of 0.00015 and tenure has −0.06711. It's tempting to think the first is less important because it's smaller. That's wrong. One is "per dollar of total charges" and the other is "per month of tenure," so you can't compare them directly. After scaling, the coefficients become −1.616, 0.913 and 0.299, all in the same unit (one standard deviation), so their sizes can now be compared.

One caution: tenure and TotalCharges rise together, so the model splits credit between them in ways that are hard to read. Comparable units don't make coefficients a clean measure of importance.

There's a third reason, and it explains why the AUC didn't move. Scikit-learn's LogisticRegression actually applies a small penalty by default, but with 7,000 rows it's too weak to matter here. Ridge and Lasso regression penalize large coefficients. A feature measured in small units needs a large coefficient to have any effect, so the penalty punishes it just because of its units, not because it's less useful.

If you skip scaling before a regularized linear model, you're not testing which features matter. You're testing which ones happen to use small numbers. Notebook 8 covers this in more detail.

Plain linear and logistic regression sit in between. Without a penalty, scaling doesn't change their predictions, but it still speeds up training and makes coefficients comparable, as you just saw.

Scaling costs almost nothing. So when a pipeline might serve more than one type of model, or you don't yet know which model will win, adding a scaler by default is a safe habit.

Section 4Scalers are transformers, so the leakage rule applies again#

The third notebook in a row where the same rule comes up

Notebook 1 set this rule for a scaler, and notebook 2 applied it again to an imputer's median. A scaler learns a mean, standard deviation, minimum, maximum or quartiles from the data it sees, and all of those must come from training rows only. California Housing's outliers make the cost of breaking this rule easy to see.

train, test = train_test_split(housing, test_size=0.2, random_state=RANDOM_STATE)

leak_check = pd.DataFrame({
    'full_mean': housing[['MedInc', 'AveOccup']].mean(),
    'train_mean': train[['MedInc', 'AveOccup']].mean(),
    'full_std': housing[['MedInc', 'AveOccup']].std(),
    'train_std': train[['MedInc', 'AveOccup']].std(),
}).round(4)
leak_check['std_differs_by'] = (leak_check['full_std'] - leak_check['train_std']).abs().round(4)

print(leak_check.to_string())

          full_mean  train_mean  full_std  train_std  std_differs_by
MedInc       3.8707      3.8808    1.8998     1.9043          0.0045
AveOccup     3.0707      3.0970   10.3860    11.5787          1.1927

MedInc barely changes: a full-data mean of 3.8707 against a training mean of 3.8808. That difference hardly matters. AveOccup is different. Its standard deviation goes from 10.386 on the full data to 11.579 on the training rows, about an 11% difference.

The reason is the outliers from section 2. Where a few extreme values land in a random split has a big effect on any statistic calculated from them. So the columns that most need careful scaling are also the ones where fitting on the wrong rows changes the numbers most.

The effect on a final score is usually small, and notebook 9 measures it. The reason to follow the rule anyway is that the same habit protects you from leaks that are big, which later notebooks show.

The fix is the same as in notebook 1: put the scaler inside a pipeline and let it fit on training data only.

leak_free = make_pipeline(StandardScaler(), Ridge())
leak_free.fit(train[['MedInc', 'AveOccup']], house_value.loc[train.index])

stored = leak_free.named_steps['standardscaler']

print(pd.DataFrame({
    'mean_stored_in_pipeline': stored.mean_,
    'train_mean': train[['MedInc', 'AveOccup']].mean().values,
    'full_data_mean': housing[['MedInc', 'AveOccup']].mean().values,
}, index=['MedInc', 'AveOccup']).round(4).to_string())

print()
print(f'The scaler learned from {len(train):,} training rows.')
print(f'It will transform the {len(test):,} test rows using those numbers, without refitting.')

          mean_stored_in_pipeline  train_mean  full_data_mean
MedInc                     3.8808      3.8808          3.8707
AveOccup                   3.0970      3.0970          3.0707

The scaler learned from 16,512 training rows.
It will transform the 4,128 test rows using those numbers, without refitting.

The mean_stored_in_pipeline column matches train_mean exactly and is different from full_data_mean. That's the whole point in one table. The fitted scaler keeps the training statistics and uses them, unchanged, on the 4,128 test rows. It never looks at the test data to decide what "normal" is.

Section 5Non-linear transformations for skewed data#

Scaling moves a distribution. It can't change its shape. This section covers the tools that can

Skew matters for some models, not all. A common myth is that linear models need normally distributed inputs. They don't. The real problems are these: in a strongly right-skewed column, a few extreme rows pull hard on the fitted line, and the relationship with the target is often closer to a straight line after a log. A linear model has to fit one slope that covers both the tail and the normal values, and it ends up describing neither well.

Tree models don't care about skew, because a log keeps the order of values, and trees only use order.

Three transformations, from gentlest to strongest:

• Square root compresses gently. It works with zero but not with negative values.
• Log compresses more strongly and is the usual first choice for money, counts and durations. It doesn't work at zero, so in practice people use log1p, which calculates log(1 + x) instead of log(x).
• Reciprocal (1/x) is so strong it's risky, and it reverses the order of all values. It's rarely what you want.

income_values = housing['MedInc'].values

candidates = {
    'original': income_values,
    'sqrt': np.sqrt(income_values),
    'log1p': np.log1p(income_values),
    'reciprocal': 1 / income_values,
}

print('Median income, skew after each transformation:')
for name, values in candidates.items():
    print(f'  {name:12s} skew = {pd.Series(values).skew():7.3f}')

Median income, skew after each transformation:
  original     skew =   1.647
  sqrt         skew =   0.689
  log1p        skew =   0.226
  reciprocal   skew =   2.448

Log brings MedInc's skew from 1.647 down to 0.226, close enough to balanced. Square root gets partway, to 0.689. The reciprocal is the interesting failure, so let's look at how it fails.

fig, axes = plt.subplots(1, 4, figsize=(17, 3.8))

for ax, (name, values) in zip(axes, candidates.items()):
    ax.hist(values, bins=60, color='#3b76af', edgecolor='white')
    ax.set(title=f'{name}\nskew {pd.Series(values).skew():.3f}', xlabel='value')

axes[0].set_ylabel('block groups')
plt.tight_layout()
plt.show()

The reciprocal made skew worse, from 1.647 to 2.448, because squeezing the tail that hard creates a new tail on the other side. It also flips the order of every value: the highest income becomes the smallest number. That reverses the direction of every relationship the column is part of. If you read a coefficient on the transformed column and forgot this, you'd reach the opposite conclusion.

Box-Cox and Yeo-Johnson#

Trying square root, then log, and comparing by hand raises a question: why guess, when you could search for the transformation that makes a column most balanced? That's what power transforms do. They use a single number called lambda and search for the lambda that makes the column look as close to a normal distribution as possible.

The Box-Cox transform is defined like this:

x^(λ) = (x^λ − 1) / λ   if λ ≠ 0,     log(x)   if λ = 0

Different lambda values give transformations you already know: 1 leaves the data almost unchanged, 0.5 acts like a square root, 0 acts like a log and −1 acts like a reciprocal. Box-Cox picks the best lambda from the data instead of you choosing one.

income_2d = housing[['MedInc']]

box_cox = PowerTransformer(method='box-cox').fit(income_2d)
yeo = PowerTransformer(method='yeo-johnson').fit(income_2d)

bc_values = box_cox.transform(income_2d).ravel()
yj_values = yeo.transform(income_2d).ravel()

print(f'Original skew          : {housing["MedInc"].skew():7.3f}')
print(f'Best manual (log1p)    : {pd.Series(np.log1p(income_values)).skew():7.3f}')
print()
print(f'Box-Cox     lambda = {box_cox.lambdas_[0]:7.4f}  ->  skew {pd.Series(bc_values).skew():7.3f}')
print(f'Yeo-Johnson lambda = {yeo.lambdas_[0]:7.4f}  ->  skew {pd.Series(yj_values).skew():7.3f}')

Original skew          :   1.647

Best manual (log1p)    :   0.226

Box-Cox     lambda =  0.0909  ->  skew   0.005
Yeo-Johnson lambda = -0.1985  ->  skew  -0.003

Box-Cox picked a lambda of 0.0909, very close to zero. So on its own, it found something very close to a log transform. It reached a skew of 0.005, compared with 0.226 for the hand-picked log, because it searches all possible values instead of trying a few named options.

One detail that matters later: scikit-learn's PowerTransformer also standardizes its output to mean 0 and standard deviation 1 by default. So a pipeline that uses it doesn't need a separate scaler.

The formula has a limit. Both log(x) and fractional powers need values above zero, so Box-Cox does too. Longitude, which is negative everywhere in California, tests this right away.

longitude = housing[['Longitude']]
print(f'Longitude range: {longitude.values.min():.2f} to {longitude.values.max():.2f}')
print()

try:
    PowerTransformer(method='box-cox').fit(longitude)
    print('Box-Cox succeeded')
except ValueError as err:
    print(f'Box-Cox raises  : ValueError: {err}')

yeo_lon = PowerTransformer(method='yeo-johnson').fit(longitude)
print(f'Yeo-Johnson     : fits fine, lambda = {yeo_lon.lambdas_[0]:.4f}')

Longitude range: -124.35 to -114.31

Box-Cox raises  : ValueError: The Box-Cox transformation can only be applied to strictly positive data

Yeo-Johnson     : fits fine, lambda = 8.0189

That error isn't a bug to work around. The method is correctly telling you it doesn't work on this column, which is better than quietly returning something wrong. Yeo-Johnson was made to remove this limit. It uses a modified formula that also accepts zero and negative values, which is why it's the default in scikit-learn's PowerTransformer.

But notice the lambda it picked: 8.0. That's an extreme reshaping, applied to a column with a skew of only −0.30 that never needed it. Running without an error isn't the same as being the right choice.

So far we've measured skew, which is useful but not the real goal. Does fixing skew actually help a model?

housing_features = ['MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup']
skewed = ['MedInc', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup']   # HouseAge is already symmetric

raw_features = housing[housing_features]
logged_features = raw_features.copy()
logged_features[skewed] = np.log1p(logged_features[skewed])

cv_reg = KFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)

def housing_r2(model, X):
    return cross_val_score(model, X, house_value, cv=cv_reg, scoring='r2').mean()


print(f"raw features                 R2 = {housing_r2(make_pipeline(StandardScaler(), Ridge()), raw_features):.4f}")
print(f"log1p on the 5 skewed cols   R2 = {housing_r2(make_pipeline(StandardScaler(), Ridge()), logged_features):.4f}")
print(f"Yeo-Johnson on all columns   R2 = {housing_r2(make_pipeline(PowerTransformer(method='yeo-johnson'), Ridge()), raw_features):.4f}")

raw features                 R2 = 0.5348
log1p on the 5 skewed cols   R2 = 0.5699

Yeo-Johnson on all columns   R2 = 0.5676

R-squared goes from 0.5348 to 0.5699. The transformed features explain about three and a half more percentage points of the variation in house prices, just from changing the shape of five columns. Yeo-Johnson applied to all six columns automatically scores 0.5676, about the same as the hand-picked logs, without anyone deciding which columns needed it. (These models leave out Latitude and Longitude on purpose. Location matters a lot for house prices, but not in a straight line, so it needs different tools.)

To be clear, Yeo-Johnson matched the careful manual work but didn't beat it. It also transformed HouseAge, which has a skew of 0.06 and never needed changing. And its pipeline had no separate scaler, because PowerTransformer standardizes its output by default.

Section 6Binning, and when it costs more than it gives#

An honest comparison, including the result most tutorials don't show

Binning turns a continuous number into ranges, like ages 0 to 17, 18 to 34, 35 to 54 and so on. The appeal is real. Each bin becomes its own column with its own coefficient, so a simple linear model can give each range a different effect. It doesn't need one straight line to fit them all. That's the same curve problem notebook 1 raised.

Three ways to choose where bins start and end:

• Equal width splits the range into same-size intervals, using pd.cut. It's easy to understand but gives very uneven bins when the data is skewed.
• Equal frequency splits so each bin has the same number of rows, using pd.qcut. The counts are balanced, but the edges land on numbers with no real meaning.
• Domain-based uses edges that mean something to the people using the model, like a contract's renewal dates.

Notebook 1 found that churn against tenure is strongly curved, and fixed it with log_tenure and tenure_squared. Binning is another way to handle the same curve, so we can compare the two on the same data.

tenure = telco['tenure']

equal_width = pd.cut(tenure, bins=5)
equal_freq = pd.qcut(tenure, q=5)

comparison = pd.DataFrame({
    'equal_width_counts': equal_width.value_counts().sort_index().values,
    'equal_width_edges': [str(i) for i in equal_width.value_counts().sort_index().index],
    'equal_freq_counts': equal_freq.value_counts().sort_index().values,
    'equal_freq_edges': [str(i) for i in equal_freq.value_counts().sort_index().index],
})

print('Same column, same number of bins, very different memberships:')
print(comparison.to_string(index=False))

Same column, same number of bins, very different memberships:
 equal_width_counts equal_width_edges  equal_freq_counts equal_freq_edges
               2459     (0.929, 15.2]               1470     (0.999, 6.0]
               1099      (15.2, 29.4]               1397      (6.0, 20.0]
                917      (29.4, 43.6]               1408     (20.0, 40.0]
                947      (43.6, 57.8]               1350     (40.0, 60.8]
               1610      (57.8, 72.0]               1407     (60.8, 72.0]

Equal width creates bins with between 917 and 2,459 customers, almost three times as many in the biggest as the smallest. Equal frequency keeps every bin close to 1,400. Where the rows fall matters more than the counts. Equal width puts all 2,459 customers from 1 to 15 months into one bin, right where churn risk changes fastest. Equal frequency spends two of its five bins on that same stretch.

other_numeric = telco[['MonthlyCharges', 'TotalCharges']]   # tenure deliberately left out, we re-add it below

def logistic_auc(X):
    return cross_val_score(make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000)),
                           X, churned, cv=cv, scoring='roc_auc').mean()


results = {}
results['raw tenure'] = logistic_auc(telco_numeric)
results['log + square (notebook 1)'] = logistic_auc(
    telco_numeric.assign(log_tenure=np.log1p(tenure), tenure_squared=tenure ** 2))

for strategy in ['uniform', 'quantile']:
    for n_bins in [5, 10]:
        kwargs = dict(n_bins=n_bins, encode='onehot-dense', strategy=strategy)
        if strategy == 'quantile':
            kwargs['quantile_method'] = 'averaged_inverted_cdf'
        binned = pd.DataFrame(KBinsDiscretizer(**kwargs).fit_transform(telco[['tenure']]),
                              index=telco.index).add_prefix('bin_')
        results[f'binned: {strategy}, {n_bins} bins'] = logistic_auc(
            pd.concat([other_numeric, binned], axis=1))

domain_bins = pd.cut(tenure, [0, 6, 12, 24, 48, 72], labels=['0-6m', '6-12m', '1-2y', '2-4y', '4y+'])
results['binned: domain, 5 bins'] = logistic_auc(
    pd.concat([other_numeric, pd.get_dummies(domain_bins, prefix='ten').astype(int)], axis=1))

scores = pd.Series(results).sort_values(ascending=False)
print('How to represent tenure, ranked by ROC AUC:')
print(scores.round(4).to_string())

How to represent tenure, ranked by ROC AUC:
log + square (notebook 1)    0.8159
binned: quantile, 10 bins    0.8155
binned: quantile, 5 bins     0.8118
binned: domain, 5 bins       0.8097
raw tenure                   0.8087
binned: uniform, 10 bins     0.8077
binned: uniform, 5 bins      0.8038

Most introductions to binning skip this result. Equal-width binning did worse than doing nothing. Five equal bins scored an AUC of 0.8038, against 0.8087 for raw tenure. We did extra work, added four columns, and ended up worse off.

The smooth transformation from notebook 1 wins, at 0.8159. Ten quantile bins come close at 0.8155, but need ten columns where the smooth version used three (tenure, its log and its square). The reason is that binning throws information away on purpose.

A customer at 3 months and one at 13 months can land in the same equal-width bin and look identical to the model, even though in notebook 1's data, churn is 56.8% in the first three months and 28.7% for customers between one and two years. Equal-frequency binning does better because it puts more edges where there's more data, but it still loses the order within each bin.

Section 7Polynomial and interaction features#

Changing the order of two steps turns the worst result in this notebook into the best one

Notebook 1 made the case for this: a linear model can only draw straight lines, and multiplying two columns lets it draw a curve. We did that by hand for one feature. PolynomialFeatures does the same thing automatically, for every column at once.

For two inputs a and b at degree 2, it creates a, b, a², ab and b². The ab term is the interaction: it captures how the effect of one variable can depend on the other. The squared terms let the model bend.

tiny = pd.DataFrame({'income': [3.0, 5.0], 'rooms': [4.0, 6.0]})

poly_demo = PolynomialFeatures(degree=2, include_bias=False)
expanded = poly_demo.fit_transform(tiny)

print('Input:')
print(tiny.to_string(index=False))
print()
print('After PolynomialFeatures(degree=2):')
print(pd.DataFrame(expanded, columns=poly_demo.get_feature_names_out()).to_string(index=False))
print()
print('2 columns became', expanded.shape[1])

Input:
 income  rooms
    3.0    4.0
    5.0    6.0

After PolynomialFeatures(degree=2):
 income  rooms  income^2  income rooms  rooms^2
    3.0    4.0       9.0          12.0     16.0
    5.0    6.0      25.0          30.0     36.0

2 columns became 5

Two columns became five. On the raw housing data, the next cell gives the most dramatic result in the notebook, and it's a bad one.

# Apply it to the housing data exactly as it comes
print(f"{'baseline: 6 raw features':38s} R2 = {housing_r2(make_pipeline(StandardScaler(), Ridge()), raw_features):8.4f}")

poly2 = PolynomialFeatures(degree=2, include_bias=False)
n_cols = poly2.fit_transform(raw_features.head()).shape[1]
naive_r2 = housing_r2(make_pipeline(poly2, StandardScaler(), Ridge()), raw_features)

print(f"{'degree 2 on raw features':38s} R2 = {naive_r2:8.4f}   ({n_cols} columns)")

baseline: 6 raw features               R2 =   0.5348
degree 2 on raw features               R2 = -11.1644   (27 columns)

A negative R-squared means the model does worse than just predicting the average house price for every block. This isn't a small disappointment. It's a total failure, and it's what happens when a standard technique is applied to data that wasn't checked first.

print('What squaring does to columns that already had extreme values:')
print()
print(f'{"column":12s} {"max":>12s} {"max squared":>18s}')
print('-' * 44)
for col in housing_features:
    top = raw_features[col].max()
    print(f'{col:12s} {top:>12.2f} {top ** 2:>18.2f}')

What squaring does to columns that already had extreme values:

column                max        max squared
--------------------------------------------
MedInc              15.00             225.00
HouseAge            52.00            2704.00
AveRooms           141.91           20138.19
AveBedrms           34.07            1160.54
Population       35682.00      1273205124.00
AveOccup          1243.33         1545877.78

Population goes up to 35,682, so its square reaches 1.27 billion. AveOccup squared reaches 1.5 million. Scaling afterward doesn't fix this. A scaler divides by the standard deviation, but that standard deviation comes from the same column, which is dominated by those huge squared values.

It gets worse at test time. When one of these extreme blocks lands in a test fold, the model multiplies its huge squared value by a coefficient learned mostly from normal rows. The prediction can be wildly off, and a few predictions like that are enough to drag the average R-squared far below zero.

The outlier stays just as extreme, and everything else gets squashed toward zero. It's the same MinMaxScaler failure from section 2, arriving another way.

# Could stronger regularisation rescue it?
for alpha in [1.0, 100.0, 10_000.0]:
    score = housing_r2(make_pipeline(PolynomialFeatures(2, include_bias=False),
                                     StandardScaler(), Ridge(alpha=alpha)), raw_features)
    print(f'degree 2, Ridge alpha = {alpha:<10g} R2 = {score:8.4f}')

degree 2, Ridge alpha = 1          R2 = -11.1644
degree 2, Ridge alpha = 100        R2 =   0.4491
degree 2, Ridge alpha = 10000      R2 =   0.4507

Turning up regularization helps a bit, recovering to 0.4507 with a Ridge alpha of 10,000. But it never gets back to 0.5348, which the six original features scored with no tuning. Regularization can't fix features that were already broken before the model saw them. The fix has to come earlier: deal with the skew first, using the tools from section 5.

print(f"{'raw features':40s} R2 = {housing_r2(make_pipeline(StandardScaler(), Ridge()), raw_features):.4f}")
print(f"{'log1p only':40s} R2 = {housing_r2(make_pipeline(StandardScaler(), Ridge()), logged_features):.4f}")
print()

for degree, interaction_only in [(2, True), (2, False), (3, False)]:
    poly = PolynomialFeatures(degree=degree, include_bias=False, interaction_only=interaction_only)
    n_cols = poly.fit_transform(logged_features.head()).shape[1]
    score = housing_r2(make_pipeline(poly, StandardScaler(), Ridge()), logged_features)
    label = f'log1p + degree {degree}' + (', interactions only' if interaction_only else '')
    print(f'{label:40s} R2 = {score:.4f}   ({n_cols} columns)')

raw features                             R2 = 0.5348
log1p only                               R2 = 0.5699

log1p + degree 2, interactions only      R2 = 0.6122   (21 columns)
log1p + degree 2                         R2 = 0.6505   (27 columns)

log1p + degree 3                         R2 = 0.6294   (83 columns)

Same columns, same model, same random seed. Squaring first, then scaling, gave an R-squared of −11.16. Taking the log first, then squaring, then scaling, gave 0.6505, the biggest gain in this notebook. Only the order changed. So here's the rule: fix skew and extreme values before creating polynomial features, because squaring an outlier makes it much worse.

Two more results from that table stand out.

Degree 3 is worse than degree 2. It scores 0.6294 with 83 columns, against 0.6505 with 27 for degree 2. More columns started to hurt, because 83 closely related columns give the model enough freedom to fit random noise in the training data. Treat the degree as a setting to test, not something to keep increasing.

Interaction terms alone aren't enough here. Using only interactions, with no squared terms, scores 0.6122 with 21 columns, short of the full 0.6505. So the squared terms add real value, which suggests that some relationships here, such as income and house value, are curved.

The biggest cost is the number of columns. Degree 2 on 6 columns gives 27. Degree 3 on the same 6 gives 83. Degree 2 on 100 columns would give 5,150, because the count grows roughly with the square of the number of columns. So PolynomialFeatures works best on small, carefully chosen sets of columns. Once you have dozens or hundreds, it's more practical to build specific combinations you have a reason to expect.

# Hand-built features: fewer, interpretable, and each one is a hypothesis you can defend
#   rooms_per_person : crowding, low values suggest cramped housing
#   bedroom_ratio    : what share of rooms are bedrooms, a proxy for layout and unit type
#   households       : how many dwellings the block group contains
ratios = pd.DataFrame({
    'rooms_per_person': housing['AveRooms'] / housing['AveOccup'],
    'bedroom_ratio': housing['AveBedrms'] / housing['AveRooms'],
    'households': housing['Population'] / housing['AveOccup'],
})
domain = pd.concat([logged_features, np.log1p(ratios)], axis=1)   # log them too, same reasoning as before

print(f"{'log1p features only':34s} R2 = {housing_r2(make_pipeline(StandardScaler(), Ridge()), logged_features):.4f}   ({logged_features.shape[1]} columns)")
print(f"{'+ 3 hand-built ratios':34s} R2 = {housing_r2(make_pipeline(StandardScaler(), Ridge()), domain):.4f}   ({domain.shape[1]} columns)")
print()
print('Compare against the 27-column polynomial expansion above.')

log1p features only                R2 = 0.5699   (6 columns)

+ 3 hand-built ratios              R2 = 0.5963   (9 columns)

Compare against the 27-column polynomial expansion above.

Here are three hand-built ratios, each with a clear meaning: crowding, the share of rooms that are bedrooms, and household count. They're compared against the 27-column automatic expansion. The automatic version still wins on R-squared, and it's important to say that plainly.

What the hand-built version gives you instead is that every column has a clear meaning, the model is small enough for a person to check, and nothing breaks unexpectedly when a new column appears next quarter.

In practice, the two approaches work well together. Hand-built features capture what you already know. Polynomial features search for patterns you didn't think of. Start with the first, and use the second once you've run out of your own ideas.