Feature Engineering · Notebook 7 companion articleFeature Construction and Domain-Driven Features
Earlier notebooks built new columns along the way, like notebook 1's service count and notebook 6's style features. This notebook makes construction the main topic and asks which new column should exist. The danger is that a new feature can look like a discovery but just be a column you already had in a new form, or worse, a copy of the answer you're trying to predict. Both failures stay silent until you check.
Start hereIntroduction#
A model can only combine features in the ways its algorithm allows. A linear model adds them up with fixed weights. A tree splits on them one at a time. Neither will create "revenue per customer" or "number of services" by itself. A tree can only approximate them with many splits. If that quantity is what drives the outcome, the model wastes effort approximating something you could have just given it. That's the case for feature construction.
It would be a short, happy notebook if building features couldn't go wrong. But it can, and in a way worse than wasted effort: a bad feature can look like a real discovery, pass a quick look at the code, and quietly rebuild a column you already had, or the answer itself.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LogisticRegression, Ridge
from sklearn.ensemble import RandomForestClassifier, HistGradientBoostingRegressor
from sklearn.model_selection import cross_val_score, StratifiedKFold, KFold, train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import roc_auc_score
from scipy import stats
plt.style.use('seaborn-v0_8-whitegrid')
pd.set_option('display.width', 200)
pd.set_option('display.max_columns', 40)
# Reused everywhere so every score in this notebook is comparable.
CV5 = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
KF5 = KFold(n_splits=5, shuffle=True, random_state=42)
print('Ready.')
Ready.
Section 1Two shapes of data#
"Construction" means two different things, depending on how your table is laid out
One row per entity. Telco has 7,032 rows and 7,032 customers. Everything about a customer is already on their row, so construction here means combining columns across the row: ratios, counts and flags.
A transaction log. The e-commerce order file, covered in section 5, has 2,000 rows but only 200 customers, because each customer appears once per order. "How much has this person spent in total" isn't a column yet. Construction here means collapsing many rows into one summary per customer, which is a very different operation.
telco = pd.read_csv('datasets/telco_churn.csv')
telco['TotalCharges'] = pd.to_numeric(telco['TotalCharges'], errors='coerce')
telco = telco[telco['tenure'] > 0].reset_index(drop=True)
y_churn = (telco['Churn'] == 'Yes').astype(int).values
orders = pd.read_csv('datasets/ecommerce_transactions.csv')
orders['order_date'] = pd.to_datetime(orders['order_date'])
print('Telco :', telco.shape, ' customers:', telco['customerID'].nunique(),
f' churn rate {y_churn.mean():.4f}')
print('Orders :', orders.shape, ' customers:', orders['customer_id'].nunique(),
' products:', orders['product_id'].nunique())
print()
print('orders per customer:')
print(orders.groupby('customer_id').size().describe()[['min', '50%', 'max']].to_string())
print()
print(orders.head(4).to_string(index=False))
Telco : (7032, 21) customers: 7032 churn rate 0.2658
Orders : (2000, 12) customers: 200 products: 50
orders per customer:
min 3.0
50% 10.0
max 19.0
order_id customer_id product_id product_category quantity unit_price discount_pct customer_age customer_income num_dependents is_returned order_date
1 103 22 Toys 3 415.08 0.03 59 52462.0 0 0 2023-12-25
2 180 39 Food 8 44.69 0.25 35 34010.0 0 0 2023-08-16
3 93 2 Books 14 241.04 0.06 32 16363.0 2 0 2023-10-20
4 15 38 Clothing 7 316.90 0.12 45 48491.0 2 0 2023-02-21Telco's customerID is unique on every row. The order log has between 3 and 19 rows per customer. That one-to-many structure is what makes aggregation possible, and, as section 6 shows, what makes it dangerous. Most real projects have both shapes. We start with the simpler one.
Section 2Arithmetic combinations, and a check people skip#
The obvious first step is to combine numeric columns. It skips a one-line check that can save you an afternoon
Telco has three numbers to combine: tenure, MonthlyCharges and TotalCharges. The average lifetime bill is TotalCharges / tenure. Whether the current bill differs from that average sounds like a natural second feature.
def churn_auc(frame, cols, model='logistic'):
# 5-fold ROC AUC on Telco. Scaler and model live in a Pipeline so nothing
# is fitted outside a fold.
if model == 'logistic':
est = Pipeline([('scale', StandardScaler()),
('clf', LogisticRegression(max_iter=3000))])
else:
est = RandomForestClassifier(n_estimators=300, min_samples_leaf=5,
random_state=42, n_jobs=1)
return cross_val_score(est, frame[cols], y_churn, cv=CV5,
scoring='roc_auc', n_jobs=1).mean()
BASE = ['tenure', 'MonthlyCharges', 'TotalCharges']
print(f'baseline, three raw numeric columns: {churn_auc(telco, BASE):.4f}')
print('(This is the same 0.8087 used as the reference point since notebook 1.)')
baseline, three raw numeric columns: 0.8087 (This is the same 0.8087 used as the reference point since notebook 1.)
telco['avg_monthly_bill'] = telco['TotalCharges'] / telco['tenure']
telco['bill_ratio'] = telco['MonthlyCharges'] / telco['avg_monthly_bill']
telco['bill_change'] = telco['MonthlyCharges'] - telco['avg_monthly_bill']
print('correlation with churn:')
for col in BASE + ['avg_monthly_bill', 'bill_ratio', 'bill_change']:
print(f' {col:20} {np.corrcoef(telco[col], y_churn)[0, 1]:+.4f}')
correlation with churn: tenure -0.3540 MonthlyCharges +0.1929 TotalCharges -0.1995 avg_monthly_bill +0.1920 bill_ratio +0.0116 bill_change +0.0022
bill_ratio has a correlation of just +0.0116 with churn, which looks like nothing. But correlation only catches straight-line relationships, so before dropping it, the notebook checks it in buckets.
buckets = pd.qcut(telco['bill_ratio'], 5)
summary = pd.DataFrame({
'churn_rate': pd.Series(y_churn).groupby(buckets, observed=True).mean().round(4),
'customers': pd.Series(y_churn).groupby(buckets, observed=True).size(),
})
print(summary.to_string())
churn_rate customers bill_ratio (0.635, 0.974] 0.2530 1407 (0.974, 0.995] 0.2070 1406 (0.995, 1.005] 0.3969 1406 (1.005, 1.029] 0.2070 1406 (1.029, 1.451] 0.2651 1407
That looks like a real finding. The middle bucket, where the current bill is almost exactly the lifetime average, churns at 39.7%, against about 21% to 27% everywhere else. That's a U-shape a single correlation number could never show. The churn difference is real, but the story it suggests is fake, and it's worth finding out why.
middle = (telco['bill_ratio'] > 0.995) & (telco['bill_ratio'] <= 1.005)
print(f'customers in the middle bucket : {middle.sum():,}')
print(f' mean tenure : {telco.loc[middle, "tenure"].mean():.1f} months')
print(f' median tenure : {telco.loc[middle, "tenure"].median():.0f} months')
print(f' share with tenure <= 2 : {(telco.loc[middle, "tenure"] <= 2).mean():.1%}')
print()
print(f'whole dataset:')
print(f' mean tenure : {telco["tenure"].mean():.1f} months')
print(f' share with tenure <= 2 : {(telco["tenure"] <= 2).mean():.1%}')
customers in the middle bucket : 1,394 mean tenure : 26.3 months median tenure : 13 months share with tenure <= 2 : 44.5% whole dataset: mean tenure : 32.4 months share with tenure <= 2 : 12.1%
44.5% of that bucket have been customers for two months or less, against 12.1% overall. The cause is math, not behavior. A customer in their first month has TotalCharges exactly equal to MonthlyCharges, so the ratio lands at exactly 1.0 automatically. The middle bucket was never "customers with a stable bill." It was "customers who just joined," and customers in their first three months churn at 56.8%, as notebook 3 noted.
The feature is mostly a disguised signal of tenure, so it shouldn't add anything to a model that already has tenure.
print(f'{"feature set":34}{"ROC AUC":>9}')
for label, cols in [
('base', BASE),
('base + avg_monthly_bill', BASE + ['avg_monthly_bill']),
('base + bill_ratio', BASE + ['bill_ratio']),
('base + bill_change', BASE + ['bill_change']),
]:
print(f'{label:34}{churn_auc(telco, cols):9.4f}')
feature set ROC AUC base 0.8087 base + avg_monthly_bill 0.8086 base + bill_ratio 0.8086 base + bill_change 0.8087
Nothing improves. Two of the three versions score slightly below the baseline, at 0.8086, the cost of giving a regularized model an extra coefficient for a column with nothing new in it. The reason applies beyond this one ratio, because of something true of all three original columns.
implied = telco['tenure'] * telco['MonthlyCharges']
resid = telco['TotalCharges'] - implied
print('Is TotalCharges just tenure x MonthlyCharges?')
print(f' correlation between them : {np.corrcoef(telco["TotalCharges"], implied)[0, 1]:.6f}')
print(f' R^2 of that identity : {1 - resid.var() / telco["TotalCharges"].var():.6f}')
print(f' median relative error : {(resid.abs() / telco["TotalCharges"]).median():.2%}')
print()
print(f' corr(TotalCharges / MonthlyCharges, tenure) = '
f'{np.corrcoef(telco["TotalCharges"] / telco["MonthlyCharges"], telco["tenure"])[0, 1]:.6f}')
Is TotalCharges just tenure x MonthlyCharges? correlation between them : 0.999560 R^2 of that identity : 0.999120 median relative error : 2.00% corr(TotalCharges / MonthlyCharges, tenure) = 0.998931
TotalCharges is tenure times MonthlyCharges, with an R-squared of 0.9991 and a median error of just 2%. Three columns that look like three pieces of information are really two. Once you know that, the earlier failure makes sense: TotalCharges / MonthlyCharges rebuilds tenure with a correlation of 0.9989, and every other ratio of these three is the same trick in a different form.
One subtlety: TotalCharges still has a use in a linear model. It's roughly the product of the other two, and section 3 shows that a linear model can't build a product by itself. So it adds a shape the model couldn't make, just not new facts.
Section 3What the model can already do for itself#
A second reason a new feature can add nothing: the model may already be able to build it
A linear model calculates a weighted sum of its inputs. If you give it x1 + x2 as a new column, you've given it something it could already make from the coefficients it fits anyway. The same goes for any difference or weighted sum. Ratios and products are different: no linear coefficients can turn x1 and x2 into x1 / x2.
Telco can't show this cleanly, since section 2 showed its columns are already redundant. California Housing, from notebook 3, has six columns that aren't built from each other, though some are related (AveRooms and AveBedrms rise together).
housing = fetch_california_housing(as_frame=True).frame
y_house = housing['MedHouseVal'].values
H_BASE = ['MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup']
h = housing[H_BASE].copy()
h['rooms_plus_bedrooms'] = h['AveRooms'] + h['AveBedrms'] # linear combination
h['rooms_minus_bedrooms'] = h['AveRooms'] - h['AveBedrms'] # linear combination
h['rooms_per_bedroom'] = h['AveRooms'] / h['AveBedrms'] # not linear
h['occupancy_x_population'] = h['AveOccup'] * h['Population'] # not linear
def house_r2(cols, model='ridge'):
if model == 'ridge':
est = Pipeline([('scale', StandardScaler()), ('m', Ridge(alpha=1.0))])
else:
est = HistGradientBoostingRegressor(max_iter=150, random_state=42)
return cross_val_score(est, h[cols], y_house, cv=KF5, scoring='r2', n_jobs=1).mean()
print('Ridge regression, R^2 to six decimals so exact ties are visible:')
print(f'{"feature set":34}{"R^2":>11}')
for label, cols in [
('base (6 columns)', H_BASE),
('+ rooms + bedrooms (sum)', H_BASE + ['rooms_plus_bedrooms']),
('+ rooms - bedrooms (diff)', H_BASE + ['rooms_minus_bedrooms']),
('+ rooms / bedrooms (ratio)', H_BASE + ['rooms_per_bedroom']),
('+ occupancy x population', H_BASE + ['occupancy_x_population']),
]:
print(f'{label:34}{house_r2(cols):11.6f}')
Ridge regression, R^2 to six decimals so exact ties are visible: feature set R^2 base (6 columns) 0.534771 + rooms + bedrooms (sum) 0.534770 + rooms - bedrooms (diff) 0.534770 + rooms / bedrooms (ratio) 0.552776 + occupancy x population 0.447769
Three different results in one table. The sum and difference change R-squared in the sixth decimal place, from 0.534771 to 0.534770. That's not a small improvement. It's the same model with a rounding error. The ratio really helps, reaching 0.552776, a gain of +0.018. The product hurts, dropping to 0.447769, a bigger loss than the ratio's gain.
It's worth checking that the sum's tiny change is exactly zero in theory, not just close to zero.
from sklearn.linear_model import LinearRegression
for name, model in [('Ridge(alpha=1.0)', Ridge(alpha=1.0)),
('LinearRegression (no penalty)', LinearRegression())]:
est = lambda: Pipeline([('scale', StandardScaler()), ('m', model)])
a = cross_val_score(est(), h[H_BASE], y_house, cv=KF5, scoring='r2', n_jobs=1).mean()
b = cross_val_score(est(), h[H_BASE + ['rooms_plus_bedrooms']], y_house,
cv=KF5, scoring='r2', n_jobs=1).mean()
print(f'{name:32} base {a:.8f} + sum {b:.8f} difference {b - a:+.2e}')
Ridge(alpha=1.0) base 0.53477078 + sum 0.53477033 difference -4.50e-07 LinearRegression (no penalty) base 0.53476658 + sum 0.53476658 difference +0.00e+00
With Ridge's penalty, the difference is −4.5×10⁻⁷, not quite zero. With plain LinearRegression and no penalty, the difference is exactly 0.00000000. That tiny leftover with Ridge was the penalty reacting to an extra, redundant coefficient, not new information. Giving a linear model a linear combination of columns it already has can't help. That's guaranteed by the math, not a quirk of this dataset.
The product's failure deserves the same close look as the ratio's success.
print(f'{"column":24}{"median":>12}{"max":>16}{"skew":>10}')
for col in ['AveOccup', 'Population', 'occupancy_x_population',
'rooms_per_bedroom', 'MedInc']:
s = h[col]
print(f'{col:24}{s.median():12.2f}{s.max():16.2f}{s.skew():10.2f}')
column median max skew AveOccup 2.82 1243.33 97.64 Population 1166.00 35682.00 4.94 occupancy_x_population 3277.45 9275266.67 111.63 rooms_per_bedroom 4.92 10.00 0.34 MedInc 3.53 15.00 1.65
Notebook 3 already flagged AveOccup as extremely skewed, at 97.6, with one neighborhood averaging 1,243 people per household at the top. Multiplying it by Population gives a column with a skew of 111.6, reaching 9.3 million against a median of 3,277.
After scaling, almost every row of this column sits near zero, and a few rows sit hundreds of standard deviations away. When one of those extreme rows lands in a test fold, the model multiplies its huge value by a coefficient learned mostly from ordinary rows, and the prediction can be wildly off. A few predictions like that pull the average R-squared down.
It's the same problem notebook 3 found when raw polynomial features pushed R-squared to −11.16: fix the skew before combining, or don't build the feature.
All of this used a linear model. Does a model that handles curves and scale by itself react differently?
print('HistGradientBoostingRegressor on the same columns:')
print(f'{"feature set":34}{"R^2":>11}')
for label, cols in [
('base (6 columns)', H_BASE),
('+ rooms + bedrooms (sum)', H_BASE + ['rooms_plus_bedrooms']),
('+ rooms / bedrooms (ratio)', H_BASE + ['rooms_per_bedroom']),
('+ occupancy x population', H_BASE + ['occupancy_x_population']),
]:
print(f'{label:34}{house_r2(cols, model="tree"):11.6f}')
HistGradientBoostingRegressor on the same columns: feature set R^2 base (6 columns) 0.698880 + rooms + bedrooms (sum) 0.699103 + rooms / bedrooms (ratio) 0.700395 + occupancy x population 0.697717
The tree barely notices. Every result is within 0.0015 of 0.6989, including the product that cost Ridge 0.087. Two things follow. First, trees handle badly scaled features much better, since a split like "is this above nine thousand" doesn't care if the column also reaches nine million. Second, trees gain less from good features: the ratio was worth +0.018 to Ridge and only +0.0015 here, because a tree can roughly build a ratio itself by stacking splits, though that costs extra depth.
So whether a new feature is worth building depends on the model. With a linear model, ratios are where the value is. With gradient boosting, expect small gains from rearranging columns, and spend most of your effort on features that bring in information the table doesn't already have. That's the next section.
Section 4Construction that adds information#
Section 2 failed because it rearranged three columns that only held two facts. The fix is to build from parts of the table the model hasn't seen
Telco has nine service columns the numeric baseline has never used: PhoneService, MultipleLines, InternetService and six more. Notebook 1 already counted them into num_services, and notebook 4 joined them into a 322-category bundle string. Here we look more closely at the count, and test whether the story it seems to tell is true.
SERVICES = ['PhoneService', 'MultipleLines', 'InternetService', 'OnlineSecurity',
'OnlineBackup', 'DeviceProtection', 'TechSupport', 'StreamingTV',
'StreamingMovies']
# These columns say "No", "No phone service" or "No internet service" when the
# customer does not have the feature. Everything else counts as subscribed.
NOT_SUBSCRIBED = ['No', 'No phone service', 'No internet service']
telco['num_services'] = sum(
(~telco[col].isin(NOT_SUBSCRIBED)).astype(int) for col in SERVICES
)
# A ratio built on the new information rather than on the redundant columns.
# Every customer here has at least one service, so no divide-by-zero guard is
# needed. On data where that is not guaranteed, this line would need one.
assert telco['num_services'].min() >= 1
telco['charge_per_service'] = telco['MonthlyCharges'] / telco['num_services']
profile = pd.DataFrame({
'customers': telco.groupby('num_services').size(),
'churn_rate': pd.Series(y_churn).groupby(telco['num_services']).mean().round(4),
'mean_monthly': telco.groupby('num_services')['MonthlyCharges'].mean().round(2),
'charge_per_service': telco.groupby('num_services')['charge_per_service'].mean().round(2),
})
print(profile.to_string())
customers churn_rate mean_monthly charge_per_service num_services 1 1260 0.1095 20.27 20.27 2 857 0.3104 42.04 21.02 3 846 0.4492 61.70 20.57 4 965 0.3648 69.36 17.34 5 921 0.3138 78.76 15.75 6 906 0.2561 86.37 14.39 7 674 0.2255 93.69 13.38 8 395 0.1241 99.43 12.43 9 208 0.0529 104.63 11.63
This is a more interesting shape than anything in section 2. Churn by number of services is an upside-down U: 11.0% at one service, rising to 44.9% at three, then falling to 5.3% at nine. charge_per_service falls steadily from about 20.3 to 11.6, a real volume discount.
That's different from MonthlyCharges alone: two customers paying 80 a month aren't comparable if one has two services and the other has eight. Both features come from columns the numeric model never had, so unlike section 2, they can actually add something.
print(f'{"feature set":44}{"logistic":>10}{"forest":>9}')
for label, cols in [
('base', BASE),
('base + num_services', BASE + ['num_services']),
('base + charge_per_service', BASE + ['charge_per_service']),
('base + num_services + charge_per_service',
BASE + ['num_services', 'charge_per_service']),
]:
print(f'{label:44}{churn_auc(telco, cols):10.4f}'
f'{churn_auc(telco, cols, model="forest"):9.4f}')
feature set logistic forest base 0.8087 0.8048 base + num_services 0.8196 0.8117 base + charge_per_service 0.8146 0.8126 base + num_services + charge_per_service 0.8203 0.8145
Now things move. num_services alone adds +0.0109 for the logistic regression, and both new columns together reach 0.8203. The random forest improves too, from 0.8048 to 0.8145. Compared with section 2, the difference is the whole point of this notebook:
| construction | information added | logistic AUC |
|---|---|---|
TotalCharges / tenure | nothing, it's MonthlyCharges | 0.8086 |
MonthlyCharges / avg bill | nothing, it's tenure | 0.8086 |
num_services | nine unused columns | 0.8196 |
MonthlyCharges / num_services | a ratio using those columns | 0.8146 |
Feature construction isn't about math. It's about deciding what the model hasn't been told yet. Every ratio in section 2 rearranged numbers the model could already reach. This one gives it nine columns the numeric baseline had never used.
It's tempting to tell a story about the upside-down U: one service is cheap and easy to forget, three services is a bill big enough to notice and reconsider, and nine services means the provider runs your whole house and leaving is a lot of work. Section 2 is why we shouldn't stop there. A pattern that looks like behavior has to be checked against the variables we already have before we believe a story about it.
check = pd.DataFrame({
'customers': telco.groupby('num_services').size(),
'churn_rate': pd.Series(y_churn).groupby(telco['num_services']).mean().round(3),
'pct_month_to_month': telco.groupby('num_services')['Contract']
.apply(lambda s: (s == 'Month-to-month').mean()).round(3),
'pct_phone_only': telco.groupby('num_services')['PhoneService']
.apply(lambda s: (s == 'Yes').mean()).round(3),
})
print(check.to_string())
print()
print('Churn by contract type, from notebook 1:')
print(pd.Series(y_churn).groupby(telco['Contract']).mean().round(3).to_string())
customers churn_rate pct_month_to_month pct_phone_only num_services 1 1260 0.110 0.439 0.937 2 857 0.310 0.629 0.867 3 846 0.449 0.896 0.830 4 965 0.365 0.753 0.842 5 921 0.314 0.651 0.897 6 906 0.256 0.472 0.923 7 674 0.226 0.300 0.963 8 395 0.124 0.157 1.000 9 208 0.053 0.029 1.000 Churn by contract type, from notebook 1: Contract Month-to-month 0.427 One year 0.113 Two year 0.028
The story doesn't hold up. The share of month-to-month customers follows the churn curve almost exactly, from 44% of one-service customers up to 90% at three services, then down to 3% at nine. Notebook 1 showed month-to-month customers churn at 42.7% and two-year customers at 2.8%. So the upside-down U is mostly a picture of who signs long contracts.
Customers with lots of services aren't mainly staying because switching is hard. Many of them are partway through a one- or two-year contract. The proper test keeps contract type fixed so it can't explain the pattern by itself.
One caution: this shows the two are tangled together, not which causes which. Customers may sign long contracts because they bundled many services. For prediction that doesn't matter. What matters is whether num_services still adds anything once the model knows the contract.
m2m = telco['Contract'] == 'Month-to-month'
within = pd.DataFrame({
'customers': telco.loc[m2m].groupby('num_services').size(),
'churn_rate': pd.Series(y_churn[m2m.values])
.groupby(telco.loc[m2m, 'num_services'].values).mean().round(4),
})
print('month-to-month customers only:')
print(within.to_string())
month-to-month customers only: customers churn_rate 1 553 0.2315 2 539 0.4824 3 758 0.4908 4 727 0.4608 5 600 0.4317 6 428 0.4486 7 202 0.4554 8 62 0.2258 9 6 0.5000
Among month-to-month customers only, the upside-down U flattens into a step. One service churns at 23.2%. Everything from two services up sits between 43% and 49%, with no downward slope. The eight- and nine-service rows show 22.6% and 50.0%, but they're based on only 62 and 6 customers, too few to read anything into.
That matters, because the striking 5.3% figure from the original table comes from 208 customers, only 6 of them month-to-month. What looked like a bundling effect was mostly a contract effect. What's left is real but smaller: one-service customers (93.7% phone-only, paying 20.27 against an overall average of 64.80) churn at about half the rate of everyone else, even on the same contract. Beyond that one group, the number of services says little once you know the contract type.
So is the feature still useful once contract type is accounted for?
telco['contract_ordinal'] = telco['Contract'].map(
{'Month-to-month': 0, 'One year': 1, 'Two year': 2})
print(f'{"feature set":42}{"ROC AUC":>9}')
for label, cols in [
('base', BASE),
('base + contract', BASE + ['contract_ordinal']),
('base + num_services', BASE + ['num_services']),
('base + contract + num_services', BASE + ['contract_ordinal', 'num_services']),
('base + contract + num_services + charge', BASE + ['contract_ordinal', 'num_services',
'charge_per_service']),
]:
print(f'{label:42}{churn_auc(telco, cols):9.4f}')
feature set ROC AUC base 0.8087 base + contract 0.8278 base + num_services 0.8196 base + contract + num_services 0.8324 base + contract + num_services + charge 0.8324
Yes, but less than the first table suggested. num_services alone added +0.0109 over the numeric baseline. With Contract already in the model, it still adds +0.0046, from 0.8278 to 0.8324. charge_per_service adds nothing at that point, 0.8324 either way. Its earlier +0.0007 was contract information sneaking in, not anything new.
For continuity with the series, here's the same feature added to notebook 1's engineered columns instead of the raw baseline:
telco['log_tenure'] = np.log1p(telco['tenure'])
telco['tenure_squared'] = telco['tenure'] ** 2
print('cumulative, matching the progression from notebook 1:')
for label, cols in [
('raw numeric', BASE),
('+ log_tenure + tenure_squared', BASE + ['log_tenure', 'tenure_squared']),
('+ num_services', BASE + ['log_tenure', 'tenure_squared', 'num_services']),
]:
print(f' {label:34}{churn_auc(telco, cols):9.4f}')
cumulative, matching the progression from notebook 1: raw numeric 0.8087 + log_tenure + tenure_squared 0.8159 + num_services 0.8241
fig, axes = plt.subplots(1, 3, figsize=(16, 4.3))
# Left: the marginal curve, and the contract mix that explains most of it.
ax = axes[0]
ax.bar(profile.index, profile['churn_rate'], color='#c53030', alpha=0.8,
label='churn rate')
ax.set_xlabel('number of services'); ax.set_ylabel('churn rate')
ax.set_title('The marginal curve and its confound', fontsize=10)
twin = ax.twinx()
twin.plot(check.index, check['pct_month_to_month'], 'o-', color='#2b6cb0', lw=2,
label='share month-to-month')
twin.set_ylabel('share on month-to-month'); twin.set_ylim(0, 1)
twin.grid(False)
lines = ax.get_legend_handles_labels()[0] + twin.get_legend_handles_labels()[0]
ax.legend(lines, ['churn rate', 'share month-to-month'], fontsize=8, loc='upper right')
# Middle: the same relationship with contract type held fixed.
ax = axes[1]
ax.bar(within.index, within['churn_rate'], color='#dd6b20', alpha=0.85)
ax.set_ylim(0, 0.6)
ax.set_title('Month-to-month only: a step, not a peak', fontsize=10)
ax.set_xlabel('number of services'); ax.set_ylabel('churn rate')
for x, (n, r) in zip(within.index, within[['customers', 'churn_rate']].values):
if n < 100:
ax.text(x, r + 0.02, f'n={int(n)}', ha='center', fontsize=7, color='#c53030')
# Right: the volume discount, which the confound does not touch.
ax = axes[2]
ax.plot(profile.index, profile['mean_monthly'], 'o-', color='#2b6cb0', label='monthly bill')
ax.plot(profile.index, profile['charge_per_service'], 's-', color='#dd6b20',
label='charge per service')
ax.set_title('More services costs more in total, less per service', fontsize=10)
ax.set_xlabel('number of services'); ax.set_ylabel('dollars'); ax.legend(fontsize=8)
plt.tight_layout()
plt.show()
The left chart makes the whole argument. The red bars are the churn curve that looked like a bundling story. The blue line is the share of month-to-month customers in each group, and it has the same shape.
The middle chart keeps contract fixed, and the peak disappears, leaving a step at one service and two thin bars that shouldn't be read as a trend. The right chart isn't affected: the volume discount is a simple fact about pricing, not a claim about behavior, which is why it never needed defending.
Section 5Group aggregations: many rows into one#
Now the transaction log, where construction means something quite different from the first four sections
The order table has no column that describes a customer. "How much has this person spent" isn't missing data waiting to be cleaned. It only exists as a property of a group of rows. Building it follows a three-step pattern used in almost every industry: group the log by the entity you care about, summarize each group, then join those summaries back, either onto the log rows or into a new table with one row per entity.
orders['line_total'] = (orders['quantity'] * orders['unit_price']
* (1 - orders['discount_pct']))
customers = orders.groupby('customer_id').agg(
n_orders=('order_id', 'count'),
total_spend=('line_total', 'sum'),
avg_order_value=('line_total', 'mean'),
max_order_value=('line_total', 'max'),
distinct_categories=('product_category', 'nunique'),
distinct_products=('product_id', 'nunique'),
avg_discount=('discount_pct', 'mean'),
first_order=('order_date', 'min'),
last_order=('order_date', 'max'),
).reset_index()
customers['days_active'] = (customers['last_order'] - customers['first_order']).dt.days
# Every customer here spans at least a month, so no divide-by-zero guard is
# needed. A customer whose orders all fall on one day would need one.
assert customers['days_active'].min() > 0
customers['spend_per_day'] = customers['total_spend'] / customers['days_active']
print(f'{len(orders):,} order rows became {len(customers):,} customer rows')
print()
# round() only the numeric columns; the two date columns are left alone.
numeric = customers.select_dtypes('number').columns
print(customers.head(6).round({c: 2 for c in numeric}).to_string(index=False))
2,000 order rows became 200 customer rows
customer_id n_orders total_spend avg_order_value max_order_value distinct_categories distinct_products avg_discount first_order last_order days_active spend_per_day
1 15 36786.22 2452.42 5325.96 5 11 0.20 2023-01-06 2023-12-16 344 106.94
2 8 25665.56 3208.20 6670.36 3 8 0.15 2023-01-17 2023-11-23 310 82.79
3 10 14540.49 1454.05 4343.33 4 10 0.15 2023-01-07 2023-11-24 321 45.30
4 10 13401.31 1340.13 3112.00 5 8 0.14 2023-01-12 2023-11-16 308 43.51
5 14 24437.06 1745.50 6062.36 5 12 0.15 2023-01-24 2023-11-23 303 80.65
6 9 10295.80 1143.98 2046.22 4 9 0.16 2023-01-29 2023-12-29 334 30.83Two patterns stand out. nunique counts variety rather than volume, and it's underused: distinct_categories and distinct_products answer a different question from a simple order count.
And days_active and spend_per_day are second-level features, built from the summaries rather than the log. Once first_order and last_order exist, the length of the relationship comes for free, and dividing total spend by it gives a rate you can compare between a two-month customer and a two-year one. Aggregation opens up new calculations that weren't possible on the original table.
The other kind of join keeps the log at its original length and copies the group summary onto each row. transform does this in one step.
# merge : entity-level table, one row per customer (what we built above)
# transform: same shape as the original log, group value broadcast to every row
orders['customer_total_spend'] = orders.groupby('customer_id')['line_total'].transform('sum')
orders['customer_n_orders'] = orders.groupby('customer_id')['line_total'].transform('count')
orders['share_of_customer_spend'] = orders['line_total'] / orders['customer_total_spend']
cols = ['customer_id', 'order_id', 'line_total', 'customer_total_spend',
'customer_n_orders', 'share_of_customer_spend']
print(orders[orders['customer_id'] == 1][cols].head(6).round(3).to_string(index=False))
print()
print('customer_total_spend repeats down the group; share_of_customer_spend does not.')
customer_id order_id line_total customer_total_spend customer_n_orders share_of_customer_spend
1 138 4057.560 36786.225 15 0.110
1 491 759.649 36786.225 15 0.021
1 537 2777.470 36786.225 15 0.076
1 710 1166.138 36786.225 15 0.032
1 742 745.921 36786.225 15 0.020
1 826 5325.965 36786.225 15 0.145
customer_total_spend repeats down the group; share_of_customer_spend does not.share_of_customer_spend is the most useful pattern in that output: a row measured against its own group. "This order is 11% of everything this customer has spent" means something that neither the order's value nor the customer's total says alone.
The same trick gives you a price compared with its category average, or a transaction compared with that account's normal size. One warning to hold on to: this total includes the customer's future orders too, so the January order already "knows" December's spending. That's fine for describing customers, but section 9 shows why it's a problem when the feature feeds a model.
So far this all works well. The next section is where it goes wrong.
Section 6Aggregating the target, and a model that scores 0.81 on noise#
The most common serious mistake in feature construction, and this dataset shows it very clearly
The order table has an is_returned column. It's the target. But it's also just a column that can be aggregated like any other, and "this customer's past return rate" sounds like a sensible feature. Before building it, let's check what's actually learnable here.
RAW_ORDER_COLS = ['quantity', 'unit_price', 'discount_pct',
'customer_age', 'customer_income', 'num_dependents']
y_ret = orders['is_returned'].values
def order_auc(cols):
est = Pipeline([('scale', StandardScaler()),
('clf', LogisticRegression(max_iter=3000))])
return cross_val_score(est, orders[cols], y_ret, cv=CV5,
scoring='roc_auc', n_jobs=1).mean()
forest = RandomForestClassifier(n_estimators=400, min_samples_leaf=5,
random_state=42, n_jobs=1)
forest_auc = cross_val_score(forest, orders[RAW_ORDER_COLS], y_ret, cv=CV5,
scoring='roc_auc', n_jobs=1).mean()
print(f'return rate : {y_ret.mean():.4f}')
print(f'logistic regression, raw columns : {order_auc(RAW_ORDER_COLS):.4f}')
print(f'random forest, raw columns : {forest_auc:.4f}')
print()
print('0.5 is a coin flip. There is essentially nothing to predict in this column.')
return rate : 0.0925 logistic regression, raw columns : 0.5156 random forest, raw columns : 0.5357 0.5 is a coin flip. There is essentially nothing to predict in this column.
Both models land close to 0.5. This is a synthetic file, and is_returned looks like it was generated with no real link to anything else. For this section that's ideal: we have a target we know for certain can't be predicted. Any model that beats 0.5 from here on is measuring a mistake, not a pattern.
orders['customer_return_rate'] = orders.groupby('customer_id')['is_returned'].transform('mean')
orders['product_return_rate'] = orders.groupby('product_id')['is_returned'].transform('mean')
orders['category_return_rate'] = orders.groupby('product_category')['is_returned'].transform('mean')
print(f'{"feature set":40}{"ROC AUC":>9}')
print(f'{"raw columns only":40}{order_auc(RAW_ORDER_COLS):9.4f}')
for label, extra in [
('+ customer_return_rate', ['customer_return_rate']),
('+ product_return_rate', ['product_return_rate']),
('+ all three rates', ['customer_return_rate', 'product_return_rate',
'category_return_rate']),
]:
print(f'{label:40}{order_auc(RAW_ORDER_COLS + extra):9.4f}')
feature set ROC AUC raw columns only 0.5156 + customer_return_rate 0.7831 + product_return_rate 0.6348 + all three rates 0.8089
0.5156 becomes 0.8089, on a target we just showed has no signal. Nothing was learned. The feature contains the answer. Each row's customer_return_rate was calculated from a group of rows that includes that same row, so each row's own label is built into its own feature. With a median of ten orders per customer, about a tenth of the feature is just the target, fed back to the model as if it were a normal input.
Cross-validation doesn't catch this, because the contamination happens before the split. Every fold's test rows already have a feature that saw their own labels. And the three rates in that table leak by different amounts, in a clear order.
print(f'{"grouped by":18}{"groups":>8}{"rows per group":>16}{"ROC AUC":>10}')
for key, label in [('customer_id', 'customer'), ('product_id', 'product'),
('product_category', 'category')]:
n_groups = orders[key].nunique()
print(f'{label:18}{n_groups:8d}{len(orders) / n_groups:16.0f}'
f'{order_auc(RAW_ORDER_COLS + [label + "_return_rate"]):10.4f}')
grouped by groups rows per group ROC AUC customer 200 10 0.7831 product 50 40 0.6348 category 5 400 0.5414
The leak shrinks as the groups get bigger. Customers average 10 rows per group, and the leaky feature scores 0.7831. Products average 40 rows and score 0.6348. Categories average 400 and score 0.5414, closest to the honest 0.5156. A row's own label makes up about 1 / group_size of its group's average, so the smaller the group, the more of the answer the model gets. That means the most detailed version, per customer, is the most tempting and the most dangerous.
The natural reaction is "fine, I'll use a proper hold-out set." That doesn't help.
train_idx, test_idx = train_test_split(np.arange(len(orders)), test_size=0.3,
stratify=y_ret, random_state=42)
def fit_and_score(frame, cols):
est = Pipeline([('scale', StandardScaler()),
('clf', LogisticRegression(max_iter=3000))]).fit(
frame.iloc[train_idx][cols], y_ret[train_idx])
train = roc_auc_score(y_ret[train_idx],
est.predict_proba(frame.iloc[train_idx][cols])[:, 1])
test = roc_auc_score(y_ret[test_idx],
est.predict_proba(frame.iloc[test_idx][cols])[:, 1])
return train, test
# Wrong: the rate was computed over every row, training and test together.
tr_bad, te_bad = fit_and_score(orders, RAW_ORDER_COLS + ['customer_return_rate'])
# Right: compute the rate on training rows only, then map it onto the test rows.
train_rate = orders.iloc[train_idx].groupby('customer_id')['is_returned'].mean()
fallback = y_ret[train_idx].mean()
orders['customer_return_rate_honest'] = orders['customer_id'].map(train_rate).fillna(fallback)
tr_ok, te_ok = fit_and_score(orders, RAW_ORDER_COLS + ['customer_return_rate_honest'])
print(f'{"":38}{"train AUC":>11}{"test AUC":>10}')
print(f'{"rate computed over all rows":38}{tr_bad:11.4f}{te_bad:10.4f}')
print(f'{"rate computed on train rows only":38}{tr_ok:11.4f}{te_ok:10.4f}')
train AUC test AUC rate computed over all rows 0.8061 0.7720 rate computed on train rows only 0.8453 0.5239
The leaky feature still scores 0.7720 on rows the model never trained on, because the leak is inside the feature, not the split. Splitting after the contamination can't undo it. Calculated correctly, from training rows only, the same idea scores 0.5239. That's the truth: barely above chance, and even that small gap is mostly noise from having only 200 customers.
The gap between 0.7720 and 0.5239 is the size of the lie. A team shipping the first version would report a strong model, watch it fail completely in production, and not know why, because their hold-out agreed with their cross-validation all along.
Now look at the training column, because it's the most useful part and the easiest to skip. The honest model has the highest training score, 0.8453, and the worst test score, 0.5239. That's not a contradiction. The training rows did contribute their own labels to the training-set rate, so the feature still memorizes them. It just can't carry that over to test rows, whose rate was looked up from other rows.
So this "honest" version is only half fixed. The test rows are clean, but each training row's rate still includes its own label, so the model learns to trust a feature that won't behave the same way on new rows. The full fix is cross-fitting, so no training row helps build its own value. That's what notebook 4's TargetEncoder does.
That wide train-to-test gap is the model correctly signaling that the feature doesn't generalize. The leaky model shows 0.8061 in training against 0.7720 in testing, a gap of 0.03 that looks, by every usual check, like a healthy, well-regularized model.
Three rules come out of this, and they're really the same rule. Never aggregate the target over a group that includes the row being scored. Target statistics by group need cross-fitting, which is exactly what TargetEncoder did in notebook 4.
Calculate group statistics from training rows only, then apply them to validation and test rows, with a clear plan for groups the training data never saw. And treat a big, unexplained jump as a reason to investigate, not celebrate: going from 0.5156 to 0.8089 by adding one column should make you suspicious.
# The first three bars are 5-fold cross-validated; the last two are single
# hold-out scores. They are different estimators, so the panels are kept apart
# rather than lined up on one axis.
cv_labels = ['raw columns\n(the truth)', 'with leaky\ncustomer rate',
'with all three\nleaky rates']
cv_values = [order_auc(RAW_ORDER_COLS),
order_auc(RAW_ORDER_COLS + ['customer_return_rate']),
order_auc(RAW_ORDER_COLS + ['customer_return_rate', 'product_return_rate',
'category_return_rate'])]
ho_labels = ['leaky feature', 'computed honestly']
ho_values = [te_bad, te_ok]
fig, axes = plt.subplots(1, 2, figsize=(13, 4.3),
gridspec_kw={'width_ratios': [3, 2]})
for ax, labs, vals, cols, title in [
(axes[0], cv_labels, cv_values, ['#2b6cb0', '#c53030', '#c53030'],
'5-fold cross-validation'),
(axes[1], ho_labels, ho_values, ['#c53030', '#2b6cb0'],
'single 30% hold-out'),
]:
bars = ax.bar(labs, vals, color=cols, alpha=0.88)
ax.axhline(0.5, color='#666666', ls='--', lw=1.4, label='chance (0.5)')
for b, v in zip(bars, vals):
ax.text(b.get_x() + b.get_width() / 2, v + 0.012, f'{v:.4f}',
ha='center', fontsize=9)
ax.set_ylim(0.4, 0.92)
ax.set_title(title, fontsize=10)
ax.legend(fontsize=8)
axes[0].set_ylabel('ROC AUC')
fig.suptitle('On a target with no signal, every bar above the dashed line is leakage',
fontsize=11)
plt.tight_layout()
plt.show()
This is the fifth time this rule has come up in the series, after imputers, scalers, encoders and vectorizers. Notebook 3's general version covers it without a special case: anything that learns from data must learn only from training rows. A group average learns one number per group. That counts as learning, just like the others.
Section 7Small groups produce unreliable statistics#
Even with no target and no leakage, a group average has a quality problem that depends on group size
A customer with 3 orders and one with 19 orders each end up with one average, but those two numbers aren't equally reliable. Nothing in the column shows that difference unless you add something to capture it.
per_customer = orders.groupby('customer_id').agg(
n_orders=('order_id', 'count'),
n_returns=('is_returned', 'sum'),
return_rate=('is_returned', 'mean'),
avg_order=('line_total', 'mean'),
)
size_bands = pd.cut(per_customer['n_orders'], [2, 5, 8, 12, 20])
reliability = per_customer.groupby(size_bands, observed=True).agg(
customers=('return_rate', 'size'),
mean_return_rate=('return_rate', 'mean'),
spread_of_rate=('return_rate', 'std'),
spread_of_avg_order=('avg_order', 'std'),
).round(4)
print(reliability.to_string())
print()
print(f'overall return rate: {y_ret.mean():.4f}')
customers mean_return_rate spread_of_rate spread_of_avg_order n_orders (2, 5] 16 0.1802 0.1590 738.2328 (5, 8] 54 0.0775 0.0969 630.6235 (8, 12] 84 0.0957 0.0917 549.3706 (12, 20] 46 0.0868 0.0710 524.2004 overall return rate: 0.0925
Customers with 3 to 5 orders have a mean return rate of 0.1802 and a spread of 0.1590. Customers with 13 to 20 orders are at 0.0868 with a spread of 0.0710, less than half as variable. We can test the spread part: if every customer had the same true return rate, how wide would each band's spread be from random variation alone?
rng = np.random.default_rng(0)
base_rate = y_ret.mean()
print(f'simulating under a single fixed return probability of {base_rate:.4f}')
print(f'{"band":12}{"customers":>11}{"observed sd":>13}{"simulated 95% range":>24}{"":>4}')
for lo, hi in [(3, 5), (6, 8), (9, 12), (13, 20)]:
in_band = (per_customer['n_orders'] >= lo) & (per_customer['n_orders'] <= hi)
sizes = per_customer.loc[in_band, 'n_orders'].values
observed = per_customer.loc[in_band, 'return_rate'].std()
sims = [np.array([rng.binomial(n, base_rate) / n for n in sizes]).std(ddof=1)
for _ in range(2000)]
low, high = np.percentile(sims, [2.5, 97.5])
verdict = 'consistent' if low <= observed <= high else 'OUTSIDE'
print(f'{f"{lo}-{hi} orders":12}{len(sizes):11d}{observed:13.4f}'
f'{f"[{low:.4f}, {high:.4f}]":>24} {verdict}')
simulating under a single fixed return probability of 0.0925 band customers observed sd simulated 95% range 3-5 orders 16 0.1590 [0.0774, 0.2054] consistent 6-8 orders 54 0.0969 [0.0853, 0.1358] consistent 9-12 orders 84 0.0917 [0.0747, 0.1064] consistent 13-20 orders 46 0.0710 [0.0607, 0.0947] consistent
Every band's actual spread falls within what random variation alone would produce. So the wider spread in smaller groups is explained by sample size alone. Estimate anything from four observations and the result is noisy, scattering around the truth in both directions with no bias.
The means are a separate question, and we should be careful not to assume the same answer for both.
# Pooled counts per band. Averaging per-customer rates weights a 3-order
# customer the same as a 19-order one, so for a significance test we go back
# to raw order counts.
band_totals = per_customer.groupby(size_bands, observed=True).agg(
orders=('n_orders', 'sum'),
returns=('n_returns', 'sum'),
)
band_totals['pooled_rate'] = (band_totals['returns'] / band_totals['orders']).round(4)
print(band_totals.to_string())
chi2, p, dof, _ = stats.chi2_contingency(
np.c_[band_totals['returns'], band_totals['orders'] - band_totals['returns']])
print(f'\nchi-square across all four bands : p = {p:.4f}')
rest = band_totals.iloc[1:].sum()
first = band_totals.iloc[0]
_, p_small, _, _ = stats.chi2_contingency(
[[first['returns'], first['orders'] - first['returns']],
[rest['returns'], rest['orders'] - rest['returns']]])
print(f'smallest band against all others : p = {p_small:.4f}')
_, p_rest, _, _ = stats.chi2_contingency(
np.c_[band_totals['returns'].iloc[1:],
(band_totals['orders'] - band_totals['returns']).iloc[1:]])
print(f'the other three bands among themselves : p = {p_rest:.4f}')
orders returns pooled_rate n_orders (2, 5] 70 13 0.1857 (5, 8] 384 29 0.0755 (8, 12] 867 83 0.0957 (12, 20] 679 60 0.0884 chi-square across all four bands : p = 0.0319 smallest band against all others : p = 0.0114 the other three bands among themselves : p = 0.5100
The honest answer here is more interesting than a neat one. The three larger bands can't be told apart statistically (p = 0.51). The smallest band, 13 returns out of 70 orders, is measurably higher than the rest (p = 0.011). If this were the only test we ran, that would be hard to put down to chance.
What we can't settle is why. Sixteen customers and seventy orders will never be enough to tell whether these people really return more, or whether this is the kind of result that shows up about once in a hundred tries, found here because we tested several bands at once.
That uncertainty is the real lesson, and it's stronger than a neat answer would be. The spread of a small group's average is inflated by sample size alone, predictably. The mean of a small group might also really be different, and at this sample size you can't tell which. Either way, a model given a raw customer_return_rate will find its most extreme values on the customers it knows least about.
Section 8RFM, the oldest feature set that still works#
Recency, Frequency, Monetary is decades older than machine learning, and its scoring step hides a real trap
Direct mail marketers were using RFM long before machine learning, working it out by hand, because three numbers capture most of what matters about a customer: Recency, how long since their last order; Frequency, how many orders they've placed; and Monetary, how much they've spent in total. It's a good example because it's three simple aggregations plus one scoring rule, and the scoring rule is where the trap is.
# The reference date is a modelling choice, not a fact. Section 9 returns to this.
reference_date = orders['order_date'].max() + pd.Timedelta(days=1)
rfm = orders.groupby('customer_id').agg(
recency=('order_date', lambda s: (reference_date - s.max()).days),
frequency=('order_id', 'count'),
monetary=('line_total', 'sum'),
).reset_index()
print(f'reference date: {reference_date.date()}')
print()
print(rfm[['recency', 'frequency', 'monetary']].describe().round(1).to_string())
reference date: 2024-01-01
recency frequency monetary
count 200.0 200.0 200.0
mean 39.6 10.0 22083.4
std 41.1 3.4 9740.7
min 1.0 3.0 3541.6
25% 11.0 8.0 14842.5
50% 27.0 10.0 20672.2
75% 54.2 12.0 27972.0
max 205.0 19.0 46612.9The usual next step turns each of the three into a score from 1 to 4 by quartile. pd.qcut is the obvious tool, but on frequency it quietly doesn't do what its name says.
freq = rfm['frequency']
print('frequency is an integer count with heavy ties:')
print(freq.value_counts().sort_index().to_string())
naive = pd.qcut(freq, 4, labels=[1, 2, 3, 4])
ranked = pd.qcut(freq.rank(method='first'), 4, labels=[1, 2, 3, 4])
print()
print('pd.qcut(frequency, 4) ->', naive.value_counts().sort_index().tolist())
print('pd.qcut(frequency.rank(), 4) ->', ranked.value_counts().sort_index().tolist())
print()
print('200 customers into 4 "quartiles" should be 50 each.')
frequency is an integer count with heavy ties: frequency 3 2 4 6 5 8 6 16 7 16 8 22 9 24 10 26 11 17 12 17 13 14 14 10 15 8 16 4 17 8 18 1 19 1 pd.qcut(frequency, 4) -> [70, 50, 34, 46] pd.qcut(frequency.rank(), 4) -> [50, 50, 50, 50] 200 customers into 4 "quartiles" should be 50 each.
Plain qcut returns groups of 70, 50, 34 and 46. Those aren't quartiles. The cause is ties: qcut cuts on values, and every customer with exactly 8 orders has to go in the same bucket, so the edges land wherever the ties allow, not where an even split would put them.
When ties are bad enough, qcut raises a duplicate-edge error. The common fix, duplicates='drop', hides the error but returns fewer bins than you asked for, still uneven. Ranking the values first, with method='first' to break ties by row order, gives exactly 50 in each group.
That fix has a cost of its own. Two customers with exactly 8 orders can now land in different groups, decided only by which row came first. Equal-sized groups are bought with a slightly arbitrary split inside the ties. That's usually acceptable for RFM, but you should know it's happening. Notebook 3 hit the same issue with equal-frequency binning. Watch for it any time you bin a count.
# Recency is inverted: fewer days since the last order is better.
rfm['R'] = pd.qcut(rfm['recency'].rank(method='first'), 4, labels=[4, 3, 2, 1]).astype(int)
rfm['F'] = pd.qcut(rfm['frequency'].rank(method='first'), 4, labels=[1, 2, 3, 4]).astype(int)
rfm['M'] = pd.qcut(rfm['monetary'].rank(method='first'), 4, labels=[1, 2, 3, 4]).astype(int)
rfm['rfm_score'] = rfm['R'] + rfm['F'] + rfm['M']
def segment(score):
# Thresholds on a 3 to 12 scale. These are a business choice, not a
# statistical one, and should be set with whoever acts on the output.
if score >= 10:
return 'Champion'
if score >= 7:
return 'Loyal'
if score >= 5:
return 'At risk'
return 'Lost'
rfm['segment'] = rfm['rfm_score'].apply(segment)
print('score distribution (3 = worst, 12 = best):')
print(rfm['rfm_score'].value_counts().sort_index().to_string())
print()
print('segment sizes:')
print(rfm['segment'].value_counts().to_string())
print()
print('segment profiles:')
print(rfm.groupby('segment')[['recency', 'frequency', 'monetary']].mean().round(1).to_string())
score distribution (3 = worst, 12 = best):
rfm_score
3 8
4 25
5 22
6 22
7 19
8 26
9 29
10 20
11 16
12 13
segment sizes:
segment
Loyal 74
Champion 49
At risk 44
Lost 33
segment profiles:
recency frequency monetary
segment
At risk 47.6 7.9 15480.4
Champion 16.1 14.1 33205.3
Lost 86.5 6.0 11471.5
Loyal 29.6 10.4 23377.3The profiles line up as they should, which confirms the scoring worked. Champions last ordered 16 days ago, placed 14 orders and spent 33,205. The Lost group last ordered 87 days ago, placed 6 orders and spent 11,471. Every column moves in the expected direction across all four groups. That's partly guaranteed, since the segments were built from these same three numbers. The check only confirms that no direction got flipped by mistake, like scoring recency the wrong way round.
fig, axes = plt.subplots(1, 3, figsize=(15, 4.2))
order = ['Champion', 'Loyal', 'At risk', 'Lost']
palette = {'Champion': '#2b6cb0', 'Loyal': '#38a169',
'At risk': '#dd6b20', 'Lost': '#c53030'}
for ax, col, title in zip(axes, ['recency', 'frequency', 'monetary'],
['Recency (days, lower is better)', 'Frequency (orders)',
'Monetary (total spend)']):
data = [rfm.loc[rfm['segment'] == s, col] for s in order]
bp = ax.boxplot(data, tick_labels=order, patch_artist=True, widths=0.6)
for patch, s in zip(bp['boxes'], order):
patch.set_facecolor(palette[s]); patch.set_alpha(0.65)
for med in bp['medians']:
med.set_color('black')
ax.set_title(title, fontsize=10)
ax.tick_params(axis='x', labelsize=9)
plt.tight_layout()
plt.show()
Two caveats before anyone uses this. First, the group cutoffs are a choice, not a discovery. Whether a score of 10 should count as Champion, instead of 9 or 11, depends on what the business plans to do with that group and how many VIPs it can afford.
Second, a single summed RFM score throws away information, since 4-4-2 and 2-4-4 both add up to 10 but describe different customers. Keep R, F and M as three separate columns for a model, and use the summed score only on dashboards, where a person is scanning lots of rows.
Section 9Aggregating over time without reading the future#
Section 8's reference date was the day after the last order in the file. That's fine for a one-off segmentation, but wrong as soon as the result feeds a model that predicts anything
Notebook 5 set the rule: a feature tied to a point in time can only use information available before that point. Group aggregates break this rule more easily than almost anything else, because a plain groupby().transform('sum') quietly includes every row in the group, even rows that happen after the row being scored.
orders_sorted = orders.sort_values(['customer_id', 'order_date', 'order_id']).reset_index(drop=True)
grp = orders_sorted.groupby('customer_id')['line_total']
# Whole-history: every row gets the same total, future orders included.
orders_sorted['spend_all_time'] = grp.transform('sum')
# Backward-only: strictly before the current row.
orders_sorted['spend_to_date'] = grp.cumsum() - orders_sorted['line_total']
orders_sorted['orders_to_date'] = grp.cumcount()
view = orders_sorted[orders_sorted['customer_id'] == 1][
['order_date', 'line_total', 'spend_all_time', 'spend_to_date', 'orders_to_date']]
print('one customer, in date order:')
print(view.to_string(index=False))
one customer, in date order: order_date line_total spend_all_time spend_to_date orders_to_date 2023-01-06 3801.2976 36786.225 0.0000 0 2023-01-22 4667.4360 36786.225 3801.2976 1 2023-01-23 825.3504 36786.225 8468.7336 2 2023-02-02 3278.8224 36786.225 9294.0840 3 2023-03-05 510.9588 36786.225 12572.9064 4 2023-05-14 2777.4705 36786.225 13083.8652 5 2023-07-10 937.8380 36786.225 15861.3357 6 2023-07-17 745.9212 36786.225 16799.1737 7 2023-09-02 759.6492 36786.225 17545.0949 8 2023-09-07 561.1848 36786.225 18304.7441 9 2023-10-12 1166.1385 36786.225 18865.9289 10 2023-11-02 5117.1978 36786.225 20032.0674 11 2023-11-21 5325.9648 36786.225 25149.2652 12 2023-12-07 4057.5600 36786.225 30475.2300 13 2023-12-16 2253.4350 36786.225 34532.7900 14
spend_all_time shows the same number, 36,786.22, on every row for this customer. So the order from 6 January already knows exactly what this customer will spend by December. A model built on that column will look great, because you've given it a time machine, not a feature. spend_to_date starts at zero and grows one order at a time. Each customer's first order correctly shows zero prior spend. That zero isn't a problem to fill in. It's exactly what was known at that moment.
print(f'mean spend_all_time : {orders_sorted["spend_all_time"].mean():>10,.0f}')
print(f'mean spend_to_date : {orders_sorted["spend_to_date"].mean():>10,.0f}')
print(f'rows where they agree: {np.isclose(orders_sorted["spend_all_time"], orders_sorted["spend_to_date"]).mean():.1%}')
print()
print('The two columns never agree on a single row, and the leaky one is more')
print('than twice as large on average.')
mean spend_all_time : 24,743 mean spend_to_date : 11,103 rows where they agree: 0.0% The two columns never agree on a single row, and the leaky one is more than twice as large on average.
The two columns don't match on a single row, and the leaky version is more than twice as big on average. The pattern to remember: cumsum() includes the current row, which is one row too many for a feature meant to describe "everything known before now." Subtracting the current row's value gives a strictly-prior total in one line.
The same applies to section 8's RFM table. To predict whether a customer churns in the fourth quarter, recency and frequency must be calculated as of the start of that quarter, not the end of whatever file is on disk. The general tools for this, covered in notebook 10, are lag and rolling-window features calculated over sorted dates, so the value at time t only uses rows from before t.
Section 10Automated feature construction#
Everything in this notebook was built by hand. Tools can create the same combinations automatically, and knowing when that's worth it matters as much as knowing the techniques
The best known is Featuretools, which uses a method called deep feature synthesis. You describe your tables and how they connect, and it applies aggregation functions (SUM, MEAN, COUNT, NUM_UNIQUE) and transform functions (DAY, MONTH, TIME_SINCE_PREVIOUS) across those connections, stacking them to a chosen depth. With this notebook's two tables, it would create SUM(orders.line_total), MEAN(orders.discount_pct), NUM_UNIQUE(orders.product_category) and hundreds more by itself. That's most of section 5 from a few lines of setup.
It's worth using on large database schemas, where writing every aggregation by hand is a week of repetitive work, and in early exploration, to see which table relationships carry any signal before you invest effort. It also handles section 9's time-cutoff problem properly through a cutoff_time setting, which is more careful than most hand-written pipelines on the first try.
It isn't always worth it. It creates hundreds or thousands of columns, most of them useless, so careful feature selection becomes a must, not an extra. It can't tell you that the churn curve it found is really a contract effect: it will happily create COUNT(services) and say nothing about what that count means, which was the whole point of section 4.
And SUM(orders.is_returned) is exactly the kind of feature it will build without hesitation. So section 6's disaster can be automated too, at scale, unless you exclude target columns before running the tool.
In practice, use both. Build the features your domain knowledge points to by hand, use automation to cover the combinations nobody has time to write, then select hard from everything it produces.
Wrap-upKey things to remember from this notebook#
Check whether columns depend on each other before building anything. Telco's TotalCharges is tenure × MonthlyCharges with an R-squared of 0.9991, so every ratio and product of those three columns just rebuilt one that already existed. None beat the 0.8087 baseline, and two scored 0.8086.
A bucket table can show a pattern that isn't real. bill_ratio seemed to find a group with 39.7% churn, against about 21% to 27% for the other buckets. That bucket was 44.5% brand-new customers, against 12.1% overall. The feature was tenure in disguise and added nothing to a model that already had tenure.
Sums and differences are easy to calculate and useless to a linear model. Adding AveRooms + AveBedrms moved Ridge's R-squared from 0.534771 to 0.534770. Ratios are different: AveRooms / AveBedrms added +0.018. AveOccup × Population cost 0.087, because multiplying amplifies skew, reaching 9.3 million against a median of 3,277. That's notebook 3's polynomial failure again.
The model decides whether a new feature pays off, and by how much. The same ratio was worth +0.018 to Ridge and only +0.0015 to gradient boosting, which can roughly build a ratio itself by stacking splits. With trees, expect small gains from rearranging columns and spend your effort elsewhere.
A new feature only helps if it adds information the model didn't have. Rearranging three redundant columns gained nothing. Counting nine unused service columns into num_services gained +0.0109. The question isn't "what can I calculate?" It's "what hasn't the model been told yet?"
Check for a hidden cause before telling the story. The upside-down U of churn by num_services, peaking near 45% and falling to 5%, suggests a story about bundles and lock-in. It's mostly contract mix: 90% of three-service customers are month-to-month, against 3% of nine-service customers, and among month-to-month customers the curve flattens into a step. The feature still adds +0.0046 on top of Contract, so it stays, but charge_per_service adds nothing once Contract is in.
Never aggregate the target over a group that includes the row being scored. On an is_returned column we proved can't be predicted, adding a group return rate took AUC from 0.5156 to 0.8089. The leaky version still scored 0.7720 on real hold-out rows, because the leak is inside the feature, not the split. Calculated from training rows only, the honest answer was 0.5239. Even that version still let each training row help build its own value. Group target statistics need cross-fitting, as TargetEncoder does.
A row's own share of its group's average is about 1 / group_size, so the leak shrinks as groups grow. Grouping by customer (10 rows each) gave 0.7831, by product (40 rows) 0.6348 and by category (400 rows) 0.5414. The most detailed version of a feature like this is always the most tempting and the most contaminated.
Leakage doesn't look like overfitting, which is why it gets past reviews. The leaky model scored 0.8061 in training and 0.7720 in testing, a gap that looks healthy. The honest model scored 0.8453 and 0.5239, and that wide gap correctly shows the feature doesn't generalize. Overfitting shows itself. Leakage hides by contaminating the test set too.
Small groups give unreliable averages, in two ways. Customers with 3 to 5 orders had a return-rate spread of 0.1590, against 0.0710 for those with 13 or more, and simulation shows random variation explains that gap. The smallest band's higher mean is a separate question the data can't fully answer: p = 0.011 on just 70 orders could be real or chance. Smoothing toward the overall average, as in notebook 4, handles both cases.
pd.qcut on tied whole numbers doesn't give you quartiles. Splitting 200 customers by order count gave groups of 70, 50, 34 and 46. duplicates='drop' hides the error without fixing the problem. Rank the values first with method='first'.
groupby().transform() reads the whole group, including future rows. A customer's January order carried their full-year spend of 36,786.22. Use cumsum() - current_value for a strictly-prior total, and remember the first row of each group correctly has nothing before it.
Notebook 8 is the other side of this one. After a whole notebook of adding columns, and knowing automated tools can add hundreds more, the next question is which ones to keep: variance thresholds, correlation, model-based importance, and why selection done outside cross-validation leaks just like section 6 did here.
This article goes with Notebook 7 of a feature engineering series built on the Telco Customer Churn dataset, the California Housing dataset and a synthetic e-commerce order log. Every code cell above runs exactly as shown in the notebook, with its real output.