6. Data Wrangling and Preparation

Raw data almost never has the shape, the types, or the cleanliness that an analysis assumes, and the work of closing that gap is where most project time actually goes. This chapter covers auditing and repairing dirty data, reshaping it into a form that supports computation, constructing the features a model will learn from, and coping with distributions that are skewed or badly imbalanced. The recurring theme is that every decision here is an assumption written into the data, so each one must be deliberate, documented, and fitted only on training rows.

<a id="61-data-cleaning"></a>

6.1 Data Cleaning

Cleaning is the process of turning values that a machine merely accepts into values that actually mean what they claim to mean; this section covers how to find the damage, how to decide what to do about it, and how to avoid inventing new damage in the process.

Profiling and Auditing Raw Data

Before you clean anything you have to know what is wrong, and the instinct to open the file and scroll is the wrong one. Scrolling shows you the first twenty rows, which are almost always the healthiest rows in the dataset because whoever built the export tested it on them. The problems live in the tail: the rows written during an outage, the rows from a client that misconfigured its integration, the rows from before a schema change.

Profiling is the systematic alternative. For each column you compute a small fixed set of facts, and you read them as a whole rather than inspecting values one at a time. The facts worth having are: the dtype, the count of nulls, the count of distinct values, the most frequent values with their counts, and for numeric columns the minimum, maximum, and a few quantiles. This is a small amount of information and it catches a surprising fraction of real defects.

The reason it works is that defects show up as implausible summary statistics rather than as visibly odd individual values. A single row reading age = 999 looks like a typo. A profile showing that 999 is the third most common value in age tells you it is a sentinel for "unknown", which is a completely different problem with a completely different fix. Likewise, a user_id column with 98,000 distinct values in 100,000 rows is not a key, and any join you were planning on it will silently fan out.

The analogy is a medical intake panel. A doctor does not begin by asking how you feel about each organ in turn. They take a standard battery of measurements, compare each against a reference range, and investigate only the ones that fall outside. Profiling is that battery, and the "reference range" is your own domain knowledge about what the column should look like.

Two habits make profiling pay off. First, profile before typing, on the raw strings, because the type conversion is itself lossy and you want to see what was there before something guessed. Second, save the profile. A profile from last month is the only way to distinguish "this column has always been 4% null" from "this column started being null on the 14th".

import pandas as pd
import numpy as np

raw = pd.DataFrame({
    "user_id":  ["1001", "1002", "1003", "1003", "1005", "1006", None],
    "age":      ["34", "29", "999", "999", "41", "-3", "58"],
    "country":  ["US", "us", "United States", "CA", "CA", "", "DE"],
    "signup":   ["2024-01-05", "2024-02-30", "05/03/2024", "2024-03-05",
                 "2024-04-01", "2024-04-02", "unknown"],
    "revenue":  ["19.99", "0", "1249.50", "1249.50", "N/A", "-15.00", "22.00"],
})

def profile(df: pd.DataFrame, top_n: int = 3) -> pd.DataFrame:
    """One row per column, the facts that catch most real defects."""
    rows = []
    for col in df.columns:
        s = df[col]
        # value_counts drops nulls by default, which is what we want here:
        # nulls are already counted separately and would swamp the top-n.
        vc = s.value_counts(dropna=True)
        rows.append({
            "dtype": str(s.dtype),
            "n_null": int(s.isna().sum()),
            "pct_null": round(100 * s.isna().mean(), 1),
            "n_unique": int(s.nunique(dropna=True)),
            # a ratio near 1.0 means "candidate key"; anything less is not one
            "unique_ratio": round(s.nunique(dropna=True) / max(len(s), 1), 2),
            "n_blank": int((s.astype("string").str.strip() == "").sum()),
            "top_values": dict(vc.head(top_n)),
        })
    return pd.DataFrame(rows, index=df.columns)

print(profile(raw).to_string())

# NUMERIC PROFILE: run it on a coerced copy, never on the stored column.
# errors="coerce" turns unparseable values into NaN so the summary is
# computable, and the count of new NaNs is itself a finding.
num = raw[["age", "revenue"]].apply(pd.to_numeric, errors="coerce")
new_nulls = num.isna().sum() - raw[["age", "revenue"]].isna().sum()
print("\nunparseable per column:\n", new_nulls)
print("\n", num.describe(percentiles=[.01, .5, .99]).round(2).to_string())

# SENTINEL DETECTION: repeated extreme values are codes, not measurements.
# A genuine measurement rarely repeats at the exact boundary of its range.
for col in num.columns:
    s = num[col].dropna()
    if s.empty:
        continue
    extreme = s[(s == s.max()) | (s == s.min())]
    if len(extreme) > 1:
        print(f"suspected sentinel in {col}: "
              f"{extreme.iloc[0]} appears {len(extreme)} times at a boundary")

# KEY CHECK: does the column you plan to join on actually identify rows?
dupes = raw["user_id"].dropna().duplicated().sum()
print(f"\nuser_id: {dupes} duplicate(s), {raw['user_id'].isna().sum()} null(s)"
      f" -> {'NOT a key' if dupes or raw['user_id'].isna().any() else 'usable key'}")

# PARSE-RATE CHECK on the date column, before committing to a format.
for fmt, label in [("%Y-%m-%d", "ISO"), ("%d/%m/%Y", "day-first")]:
    ok = pd.to_datetime(raw["signup"], format=fmt, errors="coerce").notna().mean()
    print(f"{label:10s} parses {ok:.0%} of rows")

  WHAT A PROFILE ROW TELLS YOU

  observation                          likely cause              action
  ----------------------------------------------------------------------------
  pct_null jumps vs last profile       upstream schema change    ask, don't fill
  one value dominates a numeric col    sentinel code             map to NA
  unique_ratio ~1.0 but n_null > 0     key with gaps             cannot join
  unique_ratio far below 1.0 on an ID  not a key / fan-out risk  fix grain first
  min < 0 on a quantity or age         sign or unit error        reject rows
  n_blank > 0 while n_null == 0        empty strings as nulls    normalise first
  n_unique huge on a categorical       free text, not a category treat as text
  dtype == object on a number column   one bad value in the col  find it, then fix

  ORDER OF OPERATIONS

    read as strings  ->  profile  ->  decide  ->  convert types  ->  re-profile
                            ^                                          |
                            +---------- compare the two ---------------+

Key Takeaways

• Scrolling through rows finds the defects that were already tested for; profiling finds the ones that were not.
• Profile the raw strings before type conversion, because conversion is itself lossy and hides what was originally present.
• Sentinel values such as 999 or -1 announce themselves as implausibly frequent extremes, not as odd single cells.
• Check the uniqueness ratio of any column you intend to join on, before you join on it.
• Store profiles over time; a null rate is only interpretable against its own history.

🧪 Practice

1. Write a profile function returning dtype, null count, distinct count, and top three values per column, then run it on any CSV you have and list every column whose profile surprises you.
2. Construct a dataset where a sentinel value shifts the mean by more than 20% while leaving the median unchanged, and show both numbers.
3. Interview: You are handed a 50-column table and one hour before a meeting. Which five profile statistics do you compute, and what does each one rule out? (Hint: think about which defects are invisible in a sample but obvious in an aggregate.)

Handling Missing Values

A missing value is not a value, and the single most common mistake in cleaning is treating it as one. Before you can decide what to do with a gap you have to answer a prior question: why is it missing? The answer determines whether dropping, filling, or modelling the gap is defensible, and getting it wrong biases every downstream number.

The standard taxonomy has three levels, and it is worth learning properly because the names are used loosely and the distinctions are not academic.

MCAR, missing completely at random. The probability of being missing is unrelated to anything, observed or unobserved. A sensor drops a reading because a packet was lost. Dropping these rows costs you precision and nothing else, because the rows that remain are a random sample of the rows you wanted.

MAR, missing at random. The probability of being missing depends on observed data but not on the missing value itself. Older respondents skip the income question more often, and you recorded age. This is the case that imputation is designed for: because the mechanism is explained by columns you have, a model conditioned on those columns can recover the distribution of what is absent.

MNAR, missing not at random. The probability of being missing depends on the missing value itself. High earners decline to state their income. No imputation from observed columns can fix this, because the very thing that predicts missingness is the thing you do not have. Here the only honest options are to model missingness explicitly, to collect better data, or to report the limitation.

The analogy: MCAR is a few pages torn at random out of a book. MAR is all the pages missing from chapters you can identify from the table of contents. MNAR is every page that mentioned a particular secret, removed by someone who knew which ones those were.

You cannot test for MNAR from the data alone, which is an uncomfortable but important fact. What you can do is test whether missingness in one column is predictable from the others. If it is, you have evidence for MAR and a route to imputation. If it is not, you have either MCAR or MNAR and you need domain knowledge to choose between them.

There is one more practical point that people skip: missingness is itself information. A missing last_login does not mean the login date is unknown, it usually means the user never logged in. Before imputing, ask whether the gap encodes a real state. If it does, encode that state explicitly and do not fill it.

import pandas as pd
import numpy as np

rng = np.random.default_rng(0)
n = 2000
df = pd.DataFrame({
    "age": rng.integers(18, 80, n),
    "tenure_months": rng.integers(0, 60, n),
})
df["income"] = 20000 + 600 * df["age"] + rng.normal(0, 8000, n)

# --- Construct the three mechanisms so the differences are visible. ---
mcar = df["income"].mask(rng.random(n) < 0.25)

# MAR: missingness depends on AGE, which we observe.
p_mar = np.where(df["age"] > 55, 0.55, 0.08)
mar = df["income"].mask(rng.random(n) < p_mar)

# MNAR: missingness depends on INCOME itself, which we do not observe
# when it is missing. This is the case no imputation can repair.
p_mnar = np.where(df["income"] > df["income"].quantile(0.75), 0.60, 0.05)
mnar = df["income"].mask(rng.random(n) < p_mnar)

truth = df["income"].mean()
print(f"true mean income          {truth:10,.0f}")
for name, s in [("MCAR", mcar), ("MAR", mar), ("MNAR", mnar)]:
    obs = s.mean()                      # complete-case (listwise) estimate
    print(f"{name:5s} observed mean       {obs:10,.0f}"
          f"   bias {100*(obs-truth)/truth:+6.1f}%   missing {s.isna().mean():.0%}")
# MCAR bias is noise. MAR bias is large but explainable by age.
# MNAR bias is large and NOT explainable by any column we hold.

# --- Is missingness predictable from other columns? ---
# A cheap, honest diagnostic: compare the observed columns between the
# rows where the target is present and the rows where it is absent.
def missingness_report(frame: pd.DataFrame, target: str) -> pd.DataFrame:
    flag = frame[target].isna()
    others = frame.drop(columns=[target]).select_dtypes("number")
    out = pd.DataFrame({
        "mean_when_present": others[~flag].mean(),
        "mean_when_missing": others[flag].mean(),
    })
    # Standardised difference: scale-free, so columns are comparable.
    pooled = others.std(ddof=1)
    out["std_diff"] = ((out["mean_when_missing"] - out["mean_when_present"])
                       / pooled).round(2)
    return out.round(1)

print("\nMAR case -- age differs sharply between present and missing:")
print(missingness_report(df.assign(income=mar), "income"))
print("\nMCAR case -- no column differs meaningfully:")
print(missingness_report(df.assign(income=mcar), "income"))
# |std_diff| above roughly 0.1 is the usual flag for a real imbalance.

# --- Missingness as a feature, which is often the whole answer. ---
events = pd.DataFrame({
    "user": [1, 2, 3, 4],
    "last_login": pd.to_datetime(["2024-05-01", None, "2024-05-03", None]),
})
# WRONG: filling with a date invents logins that never happened.
# RIGHT: the gap is a state. Name it.
events["never_logged_in"] = events["last_login"].isna()
events["days_since_login"] = (
    (pd.Timestamp("2024-05-10") - events["last_login"]).dt.days
)   # stays NaN where the state flag already carries the meaning
print("\n", events)

# --- Column-level and row-level triage. ---
frame = df.assign(income=mar, notes=np.nan)
col_null = frame.isna().mean().sort_values(ascending=False)
print("\nnull rate by column:\n", col_null.round(3))
drop_cols = col_null[col_null > 0.90].index.tolist()   # nothing to recover
print("columns beyond saving:", drop_cols)

# Drop ROWS only when the loss is small AND provably unbiased.
rowwise = frame.drop(columns=drop_cols)
keep = rowwise.dropna()
print(f"listwise deletion keeps {len(keep)}/{len(rowwise)} rows"
      f" ({len(keep)/len(rowwise):.0%})")
print("age mean before/after:",
      round(rowwise['age'].mean(), 1), round(keep['age'].mean(), 1))
# The shift in age mean IS the bias listwise deletion introduced.

  CHOOSING A RESPONSE TO MISSINGNESS

  mechanism  detectable from data?   listwise delete   impute   model the gap
  ---------------------------------------------------------------------------
  MCAR       partly (no signal)      safe, costs n     safe     unnecessary
  MAR        yes (signal in cols)    BIASED            correct  helpful
  MNAR       no                      BIASED            BIASED   only option

  DECISION PATH

    is the gap a real state?  --yes-->  encode a flag, do NOT fill
             | no
             v
    is missingness predictable from other columns?
             |                                  \
            yes (MAR)                            no
             |                                    \
      impute conditionally                 MCAR or MNAR?  -- domain call
      + add a missing-indicator            /            \
                                     drop rows      report the limitation

Key Takeaways

• The mechanism, not the percentage, determines the correct treatment.
• Listwise deletion is unbiased only under MCAR; under MAR it silently skews every estimate.
• MNAR cannot be detected from the data alone and cannot be imputed away.
• Test for MAR by checking whether observed columns differ between rows where the target is present and rows where it is absent.
• A gap that encodes a real state, such as "never logged in", should become an explicit flag rather than a filled value.

🧪 Practice

1. Simulate a column that is MCAR at 30% and one that is MAR at 30%, then show that complete-case analysis is unbiased for the first and biased for the second.
2. Write a function that, for a given target column, reports the standardised difference of every other numeric column between present and missing rows, and flag any above 0.1.
3. Interview: A model's accuracy drops sharply in production, and you find that a feature which was 2% null in training is now 40% null. What are the two distinct things that could have broken? (Hint: one is about the data and one is about what your pipeline did with the nulls.)

Imputation Strategies

Once you have decided a gap should be filled rather than flagged or dropped, the question becomes what to fill it with. Every imputation is a claim about the unobserved value, and the cost of a bad claim is not just inaccuracy: it is understated uncertainty, because the filled value looks exactly as solid as a measured one to everything downstream.

Start with the simplest methods and understand precisely what they break.

Mean or median imputation replaces gaps with a single constant. It preserves the central tendency of the column and destroys everything else. Variance shrinks, because you have added mass at the centre. Correlations with other columns shrink toward zero, because the filled rows carry no relationship at all. A histogram develops a spike. For a column that is 2% null this is usually tolerable; at 30% it corrupts the covariance structure that any linear model depends on.

Mode imputation is the categorical equivalent and has the same problem, plus one of its own: it inflates the majority class, which matters if that category is also predictive.

Conditional imputation uses other columns. Fill income from the median income within the same age band and region rather than the global median. This is the first method that respects MAR, and it is often enough. It is also cheap and easy to explain, which matters more than people admit.

Model-based imputation trains a predictor for the missing column using the others. K-nearest-neighbours imputation finds the most similar complete rows and borrows from them. Iterative imputation, sometimes called MICE, cycles through columns fitting each as a function of the rest, repeating until values stabilise. These recover relationships the simple methods destroy, at the cost of being slow and of being able to leak if fitted incorrectly.

Multiple imputation is the statistically correct answer when uncertainty matters. Instead of one filled dataset you create several, each with different plausible draws, run your analysis on each, and combine the results so that the final standard errors include the uncertainty about the imputation itself. Single imputation reports intervals that are too narrow, always.

For time series the right defaults are different. Forward fill carries the last observation forward, which is correct for a state that persists, such as a price or a configuration setting. Interpolation is correct for a quantity that moves smoothly, such as a temperature. Using the column mean for a time series is almost always wrong, because it ignores the ordering that gives the data its meaning.

Two rules override all of the above.

First: fit the imputer on training data only. The median you fill with is a statistic estimated from data. Computing it over the full dataset and then splitting leaks test information into training, and the resulting score is optimistic. This is the single most common leakage bug in practice.

Second: add a missing indicator. A boolean column recording that the value was imputed costs one column and lets the model learn that missingness itself predicted something. If missingness was uninformative the model will ignore it.

import numpy as np
import pandas as pd
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer, KNNImputer
from sklearn.experimental import enable_iterative_imputer   # noqa: F401
from sklearn.impute import IterativeImputer
from sklearn.model_selection import train_test_split
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_absolute_error

rng = np.random.default_rng(7)
n = 1500
age = rng.integers(20, 70, n)
hours = 20 + 0.3 * age + rng.normal(0, 4, n)
income = 15000 + 700 * age + 900 * hours + rng.normal(0, 6000, n)
X = pd.DataFrame({"age": age, "hours": hours.round(1)})
y = pd.Series(income, name="income")

# Punch MAR holes in `hours`: older rows are likelier to be missing.
holes = rng.random(n) < np.where(age > 50, 0.5, 0.05)
X.loc[holes, "hours"] = np.nan
print(f"missing in hours: {X['hours'].isna().mean():.0%}")

# --- What mean imputation destroys, measured. ---
complete = pd.DataFrame({"age": age, "hours": hours})
print(f"\ntrue  sd(hours)={complete['hours'].std():.2f}"
      f"  corr(age,hours)={complete['age'].corr(complete['hours']):+.3f}")

mean_filled = X.fillna(X.mean())
print(f"mean  sd(hours)={mean_filled['hours'].std():.2f}"
      f"  corr(age,hours)={mean_filled['age'].corr(mean_filled['hours']):+.3f}")

# Conditional median by age band keeps far more of the structure.
band = pd.cut(X["age"], bins=[19, 30, 40, 50, 60, 70])
cond = X.copy()
cond["hours"] = cond["hours"].fillna(cond.groupby(band, observed=True)["hours"]
                                     .transform("median"))
print(f"cond  sd(hours)={cond['hours'].std():.2f}"
      f"  corr(age,hours)={cond['age'].corr(cond['hours']):+.3f}")

# --- The leakage trap, shown concretely. ---
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0)

# WRONG: statistic computed over train AND test before splitting.
leaky_median = X["hours"].median()          # <- uses test rows
# RIGHT: statistic computed on train only, then applied to test.
clean_median = Xtr["hours"].median()
print(f"\nmedian from all data {leaky_median:.2f} vs train only {clean_median:.2f}")

# A Pipeline makes the right thing the default: fit() sees only train.
def build(imputer, add_flag: bool) -> Pipeline:
    return Pipeline([
        ("impute", ColumnTransformer([
            ("num", SimpleImputer(strategy="median", add_indicator=add_flag)
                    if isinstance(imputer, str) else imputer,
             ["age", "hours"]),
        ])),
        ("model", Ridge(alpha=1.0)),
    ])

candidates = {
    "median":            build("median", add_flag=False),
    "median + flag":     build("median", add_flag=True),
    "knn (k=5)":         build(KNNImputer(n_neighbors=5, add_indicator=True), True),
    "iterative (MICE)":  build(IterativeImputer(max_iter=10, random_state=0,
                                                add_indicator=True), True),
}
for name, pipe in candidates.items():
    pipe.fit(Xtr, ytr)                       # imputer fitted on TRAIN only
    err = mean_absolute_error(yte, pipe.predict(Xte))
    print(f"{name:18s} MAE {err:9,.0f}")

# --- Time series: ordering changes the right default. ---
ts = pd.Series(
    [10.0, np.nan, np.nan, 13.0, 14.0, np.nan, 16.0],
    index=pd.date_range("2024-01-01", periods=7, freq="D"), name="temp_c",
)
print("\n", pd.DataFrame({
    "raw": ts,
    "ffill": ts.ffill(),                        # a state that persists
    "interpolate": ts.interpolate("time"),      # a quantity that moves
    "mean": ts.fillna(ts.mean()),               # wrong: ignores order
}).to_string())

# LIMIT the reach of a forward fill, or one gap propagates indefinitely.
print("\nbounded ffill:\n", ts.ffill(limit=1).to_string())

  WHAT EACH IMPUTER PRESERVES

  method               mean   variance   correlations   cost      leak risk
  ---------------------------------------------------------------------------
  drop rows            ~      yes        yes            loses n   none
  constant (0, -1)     NO     NO         NO             trivial   none
  mean / median        yes    NO         NO             trivial   if fit on all
  mode (categorical)   ~      NO         NO             trivial   if fit on all
  conditional median   yes    partly     partly         low       if fit on all
  KNN                  yes    mostly     mostly         O(n^2)    if fit on all
  iterative / MICE     yes    mostly     yes            high      if fit on all
  multiple imputation  yes    yes        yes            highest   if fit on all

  THE VARIANCE COLLAPSE, DRAWN

    original          after mean imputation
      .  . .            .  .  |  .
    . . . . . .        . . . .|| . . .      <- spike at the mean
   .............      .........||.........
   low   mean  high   low    mean    high

  RULE: fit on train, transform everywhere; and always keep the indicator.

Key Takeaways

• Imputation trades a visible gap for an invisible assumption, and single imputation always understates uncertainty.
• Mean and median filling preserve the centre while shrinking variance and flattening correlations.
• Conditional and model-based imputation respect MAR by using the columns that explain the missingness.
• Time series call for forward fill or interpolation, chosen by whether the quantity persists or moves; bound the fill so one gap cannot spread.
• Fit every imputer inside a pipeline on training data only, and add a missing indicator so the model can use the gap itself.

🧪 Practice

1. Take a numeric column, remove 25% of it at random, impute with the mean, and report the change in standard deviation and in its correlation with another column.
2. Build two scikit-learn pipelines for the same data, one imputing with the median fitted on the full dataset and one fitted inside the pipeline, and compare their cross-validated scores.
3. For a daily time series with a five-day gap, produce forward-filled and time-interpolated versions and state which you would defend for a price series and which for a temperature series.
4. Interview: Your median imputation improves cross-validated accuracy but the model performs worse in production. What is the most likely cause? (Hint: what did the imputer learn from, and what does production look like that your validation split did not?)

Outlier Detection and Treatment

An outlier is a value far from the rest of the data. That definition is descriptive and useless on its own, because it says nothing about what you should do. The useful question is different: is this value wrong, or is it rare? Those two cases have opposite correct responses, and no statistical test can tell them apart.

A temperature reading of 450 degrees Celsius in a room is wrong. The same value in a kiln is the measurement you cared about. A transaction of ten million is an error in a coffee shop's ledger and the single most important row in a fraud dataset. Deleting outliers without establishing which case you are in is how people remove the signal and keep the noise.

So the workflow is: detect, then diagnose, then decide. Detection is mechanical. Diagnosis needs domain knowledge. Only the decision is a technique.

For detection, the methods differ in what they assume.

Z-score measures distance from the mean in standard deviations. Its flaw is that the mean and standard deviation are themselves dragged by the outliers you are trying to find, an effect called masking. With one extreme value in a small sample, no z-score can exceed a bound set by the sample size, so the test can structurally fail to flag anything.

The IQR rule uses quartiles, which extreme values cannot move, and flags anything beyond 1.5 times the interquartile range from the nearest quartile. It is robust and it is the sensible default for one column at a time. Its weakness is that on a skewed distribution it flags a large fraction of the long tail, which is not an error, it is the shape of the data.

Modified z-score uses the median and the median absolute deviation instead of the mean and standard deviation, giving a robust score you can threshold. The conventional cutoff is 3.5.

All three are univariate, and this is their shared blind spot. A person who is 1.5 metres tall is unremarkable. A person who weighs 120 kilograms is unremarkable. A person who is 1.5 metres tall and weighs 120 kilograms is a multivariate outlier that no single-column method can see. For that you need Mahalanobis distance, which measures distance in units of the data's own covariance, or a model such as isolation forest or local outlier factor, which learn what "typical" looks like jointly.

For treatment, the options in increasing aggression:

• Investigate and fix. The best outcome. A unit error, a misplaced decimal, a sentinel value: all repairable once identified.
• Leave it. Correct whenever the value is real and your method tolerates it. Tree-based models barely care about outliers in features.
• Winsorize or clip. Cap values at a percentile. Keeps the row, bounds the influence, distorts the tail.
• Transform. A log transform compresses a right tail so that extreme values stop dominating without being altered relative to each other.
• Use a robust method. Median instead of mean, quantile regression instead of least squares, Huber loss instead of squared loss. Often this is better than touching the data at all.
• Drop the row. Last resort, and only for values you have established are impossible. Always record how many you dropped.

One more thing that deserves emphasis: clip bounds are parameters fitted on training data, exactly like an imputation median. Computing the 1st and 99th percentile over your whole dataset before splitting leaks.

import numpy as np
import pandas as pd

rng = np.random.default_rng(3)
clean = rng.normal(50, 5, 200)
x = pd.Series(np.concatenate([clean, [140.0, 145.0, -20.0]]), name="value")

# --- Z-SCORE MASKING: the outliers inflate the denominator. ---
z = (x - x.mean()) / x.std(ddof=0)
print(f"mean {x.mean():.1f}  sd {x.std():.1f}")
print(f"z-score flags (|z|>3): {(z.abs() > 3).sum()}")

# --- ROBUST ALTERNATIVES: quartiles and MAD cannot be dragged. ---
q1, q3 = x.quantile([0.25, 0.75])
iqr = q3 - q1
lo, hi = q1 - 1.5 * iqr, q3 + 1.5 * iqr
print(f"\nIQR fence [{lo:.1f}, {hi:.1f}] flags {((x < lo) | (x > hi)).sum()}")

med = x.median()
mad = (x - med).abs().median()
# 0.6745 rescales MAD so the score is comparable to a normal z-score.
mod_z = 0.6745 * (x - med) / mad
print(f"modified z (>3.5) flags {(mod_z.abs() > 3.5).sum()}")
print("\ncomparison of flagged values:")
print(pd.DataFrame({"value": x, "z": z.round(2), "mod_z": mod_z.round(2)})
      .loc[mod_z.abs() > 3.5].to_string())

# --- MULTIVARIATE: neither column is extreme, the PAIR is. ---
h = rng.normal(1.70, 0.08, 300)
w = 25 * h ** 2 + rng.normal(0, 4, 300)        # weight tracks height
people = pd.DataFrame({"height_m": h.round(2), "weight_kg": w.round(1)})
people.loc[len(people)] = [1.50, 120.0]        # plausible alone, not together

for col in people.columns:                      # univariate sees nothing
    c1, c3 = people[col].quantile([0.25, 0.75])
    f = 1.5 * (c3 - c1)
    flagged = ((people[col] < c1 - f) | (people[col] > c3 + f)).sum()
    print(f"\nunivariate IQR on {col}: {flagged} flagged")

cov_inv = np.linalg.pinv(np.cov(people.values, rowvar=False))
centred = people.values - people.values.mean(axis=0)
# Mahalanobis: distance measured in units of the data's own covariance,
# so it accounts for the height/weight relationship the columns share.
md = np.sqrt(np.einsum("ij,jk,ik->i", centred, cov_inv, centred))
people["mahalanobis"] = md.round(2)
print("\ntop 3 by Mahalanobis distance:")
print(people.nlargest(3, "mahalanobis").to_string())

# --- Model-based, for many columns at once. ---
from sklearn.ensemble import IsolationForest
iso = IsolationForest(contamination=0.01, random_state=0)
people["iso_outlier"] = iso.fit_predict(people[["height_m", "weight_kg"]]) == -1
print(f"\nisolation forest flagged {people['iso_outlier'].sum()} row(s)")

# --- TREATMENT: four responses, same column, different consequences. ---
from sklearn.model_selection import train_test_split
tr, te = train_test_split(x, test_size=0.3, random_state=0)

p01, p99 = tr.quantile([0.01, 0.99])            # bounds from TRAIN only
treated = pd.DataFrame({
    "raw": te,
    "clipped": te.clip(p01, p99),               # winsorize: keep row, cap value
    "log1p": np.log1p(te.clip(lower=0)),        # compress the tail
})
print("\nmean / sd under each treatment:")
print(treated.agg(["mean", "std"]).round(2).to_string())

# Robust statistics: change the estimator, not the data.
print(f"\nmean {x.mean():.2f}  vs  median {x.median():.2f}"
      f"  vs  10% trimmed mean {x.sort_values().iloc[20:-20].mean():.2f}")

# ALWAYS LOG WHAT YOU REMOVED. An undocumented drop is unreproducible.
mask_impossible = x < 0                          # a negative measurement here
print(f"\ndropping {mask_impossible.sum()} impossible value(s): "
      f"{x[mask_impossible].tolist()}")

  DETECTION METHOD SELECTION

  method            assumes        robust?   dimensions   flags on skew
  --------------------------------------------------------------------------
  z-score > 3       normality      NO        one          moderate
  IQR fence 1.5x    nothing        yes       one          MANY (false)
  modified z > 3.5  nothing        yes       one          moderate
  Mahalanobis       elliptical     no        many         n/a
  isolation forest  nothing        yes       many         n/a
  local outlier f.  nothing        yes       many         n/a

  THE MASKING EFFECT

    without outlier:   mean=50  sd=5    a value of 70 sits at z = 4.0  FLAGGED
    with two at 140:   mean=51  sd=13   the same 70 sits at z = 1.5    MISSED

  DECIDE BEFORE YOU TREAT

    is the value physically possible?
        no  -> it is an error. fix it, or drop it and log the count.
        yes -> does your method's conclusion depend on it?
                 no  -> leave it alone.
                 yes -> transform, or clip, or switch to a robust estimator.
                        Removing it means answering a different question.

Key Takeaways

• "Wrong" and "rare" look identical statistically and demand opposite treatments, so diagnosis must come between detection and action.
• Z-scores are dragged by the very values they should flag; prefer the IQR fence or a modified z-score built on the median and MAD.
• Univariate methods cannot see a combination that is only jointly implausible, which is what Mahalanobis distance and isolation forests are for.
• Prefer a robust estimator or a transform over editing data; drop rows only for values established as impossible, and log the count.
• Clip and winsorize bounds are fitted parameters, so compute them on training data only.

🧪 Practice

1. Add two extreme values to a normal sample and show that the z-score rule flags fewer points than the modified z-score rule does.
2. Build a two-column dataset in which a row is flagged by Mahalanobis distance while passing the IQR test on both columns individually.
3. Fit the same linear regression on data that is raw, winsorized at the 1st and 99th percentile, and log-transformed, then compare the coefficients and explain which question each one answers.
4. Interview: A stakeholder asks you to remove outliers before reporting average order value. What do you ask them first? (Hint: what decision will the number be used for, and are the large orders real?)

Deduplication and Record Linkage

Duplicates are two kinds of problem wearing the same name, and conflating them wastes a lot of time. Exact duplicates are the same row recorded twice, usually from a pipeline retry or a double-submitted form. Fuzzy duplicates are two records describing the same real-world entity with different text: "Robert Smith" and "Bob Smith Jr.", or "Acme Inc." and "ACME Incorporated". The first is a data-integrity bug you fix with a key. The second is an inference problem called record linkage or entity resolution, and it never has a provably correct answer.

For exact duplicates the essential move is to define the grain of your table: the set of columns that should uniquely identify a row. An orders table might be one row per order_id. An events table might be one row per (user_id, event_type, timestamp). Once the grain is written down, duplicates are just violations of it, and you can assert the grain in the pipeline so violations fail loudly instead of quietly doubling your revenue figure.

Deduplicating is not only about which rows to remove but which to keep. If two rows share a key but differ in other columns, dropping "the duplicate" silently picks one. Usually you want the most recent, or the most complete, and you should say which. Sorting before dropping makes the choice explicit and reproducible.

Fuzzy matching is harder because comparing every record to every other record is quadratic. Ten thousand records is fifty million comparisons; a million records is five hundred billion. The standard solution is blocking: only compare records that agree on some cheap key, such as the first letter of the surname plus the postcode. Blocking trades recall for tractability, which is why you often block on several different keys and take the union.

Within a block you compute similarity on each field. The main string metrics:

Metric	Measures	Good for
Levenshtein	edit operations	typos, short strings
Jaro-Winkler	transpositions, prefix match	names
Jaccard on tokens	shared word sets	addresses, company names
Soundex / Metaphone	pronunciation	names captured by ear
Cosine on n-grams	shared character chunks	long messy strings, scalable

Then you combine field similarities into a decision. A weighted sum with a threshold is the simple version. The classical probabilistic version is the Fellegi-Sunter model, which weights each field by how informative agreement on it is: two people sharing a surname is weak evidence, two people sharing a phone number is strong evidence. Modern tooling trains a classifier on labelled pairs.

Finally, matches are transitive but your scores are not. If A matches B and B matches C, you have implicitly claimed A matches C even if their score was below threshold. Resolving this means treating matches as edges in a graph and taking connected components, which is also how you get an entity-level identifier to join on.

import pandas as pd
import numpy as np
from itertools import combinations

# --- EXACT DUPLICATES: define the grain, then assert it. ---
orders = pd.DataFrame({
    "order_id": [1, 2, 2, 3, 4, 4],
    "user_id":  [10, 11, 11, 12, 13, 13],
    "amount":   [50.0, 20.0, 20.0, 75.0, 30.0, 35.0],   # note: 4 differs
    "updated":  pd.to_datetime(["2024-01-01", "2024-01-02", "2024-01-02",
                                "2024-01-03", "2024-01-04", "2024-01-05"]),
})

def assert_grain(df: pd.DataFrame, keys: list[str]) -> None:
    """Fail loudly if the table is not one row per `keys`."""
    dup = df.duplicated(subset=keys, keep=False)
    if dup.any():
        raise ValueError(
            f"grain violated on {keys}: {dup.sum()} rows in "
            f"{df.loc[dup, keys].drop_duplicates().shape[0]} groups")

try:
    assert_grain(orders, ["order_id"])
except ValueError as e:
    print("caught:", e)

# FULL-ROW duplicates are safe to drop: no information is lost.
print("\nidentical rows:", orders.duplicated().sum())

# KEY duplicates that DIFFER need a rule. Make the rule visible.
print("\nconflicting on key:")
print(orders[orders.duplicated("order_id", keep=False)].to_string(index=False))

# Sort first so `keep="last"` means something specific and reproducible.
latest = (orders.sort_values(["order_id", "updated"])
                .drop_duplicates("order_id", keep="last"))
print("\nkeep most recent:\n", latest.to_string(index=False))

# Alternative rule: keep the most complete row.
orders["completeness"] = orders.notna().sum(axis=1)
most_complete = (orders.sort_values(["order_id", "completeness"])
                       .drop_duplicates("order_id", keep="last"))

# Or aggregate instead of choosing, when both rows are real events.
collapsed = orders.groupby("order_id", as_index=False).agg(
    user_id=("user_id", "first"),
    amount=("amount", "sum"),          # <- a DIFFERENT claim about the data
    last_seen=("updated", "max"),
)
print("\naggregated (treats both rows as real):\n", collapsed.to_string(index=False))

# --- FUZZY DUPLICATES: normalise, block, score, resolve. ---
people = pd.DataFrame({
    "rid":   [1, 2, 3, 4, 5, 6],
    "name":  ["Robert Smith", "Bob Smith", "Roberta Smythe",
              "Acme Inc.", "ACME Incorporated", "Zenith Ltd"],
    "email": ["r.smith@x.com", "rsmith@x.com", "rsmythe@y.com",
              "info@acme.com", "info@acme.com", "hi@zenith.com"],
    "zip":   ["10001", "10001", "10001", "94107", "94107", "60601"],
})

def normalise(s: pd.Series) -> pd.Series:
    """Cheap canonical form: the single highest-value step in linkage."""
    return (s.str.lower()
             .str.replace(r"[^a-z0-9\s]", "", regex=True)   # drop punctuation
             .str.replace(r"\b(inc|incorporated|ltd|llc|jr|sr)\b", "", regex=True)
             .str.replace(r"\s+", " ", regex=True)
             .str.strip())

people["key"] = normalise(people["name"])
print("\nnormalised:\n", people[["rid", "name", "key"]].to_string(index=False))

# BLOCKING: compare only within a cheap shared key, here the postcode.
# Without it, n records cost n*(n-1)/2 comparisons.
def jaccard(a: str, b: str) -> float:
    ta, tb = set(a.split()), set(b.split())
    return len(ta & tb) / len(ta | tb) if ta | tb else 0.0

def bigrams(s: str) -> set[str]:
    s = s.replace(" ", "")
    return {s[i:i + 2] for i in range(len(s) - 1)}

def bigram_sim(a: str, b: str) -> float:
    ba, bb = bigrams(a), bigrams(b)
    return len(ba & bb) / len(ba | bb) if ba | bb else 0.0

pairs = []
for _, block in people.groupby("zip"):
    for (_, a), (_, b) in combinations(block.iterrows(), 2):
        # Weighted evidence: an exact email match is far stronger than a
        # name resemblance, so it gets the larger weight.
        name_s = (0.5 * jaccard(a["key"], b["key"])
                  + 0.5 * bigram_sim(a["key"], b["key"]))
        email_s = float(a["email"] == b["email"])
        score = 0.45 * name_s + 0.55 * email_s
        pairs.append((a["rid"], b["rid"], round(name_s, 2), email_s, round(score, 2)))

scored = pd.DataFrame(pairs, columns=["a", "b", "name_sim", "email_eq", "score"])
print("\ncandidate pairs within blocks:\n", scored.to_string(index=False))

n_all = len(people) * (len(people) - 1) // 2
print(f"\nblocking compared {len(scored)} pairs instead of {n_all}")

# --- TRANSITIVE CLOSURE: matches form a graph, entities are components. ---
THRESHOLD = 0.55
parent = {r: r for r in people["rid"]}

def find(x):
    while parent[x] != x:
        parent[x] = parent[parent[x]]        # path compression
        x = parent[x]
    return x

def union(x, y):
    rx, ry = find(x), find(y)
    if rx != ry:
        parent[max(rx, ry)] = min(rx, ry)    # smallest rid becomes entity id

for _, row in scored[scored["score"] >= THRESHOLD].iterrows():
    union(int(row["a"]), int(row["b"]))

people["entity_id"] = [find(r) for r in people["rid"]]
print("\nresolved entities:\n",
      people[["rid", "name", "entity_id"]].to_string(index=False))
# Now join on entity_id, not on name. The cluster, not the pair, is the unit.

  TWO DIFFERENT PROBLEMS

  exact duplicates                    fuzzy duplicates
  ----------------------------------  ----------------------------------
  cause: retry, double submit, bad    cause: free-text entry, multiple
         join fan-out                        source systems, name changes
  detect: duplicated(subset=grain)    detect: similarity above a threshold
  fix:    drop, or aggregate          fix:    cluster into entities
  correct answer exists               correct answer is a judgement call
  assert it in the pipeline           measure precision AND recall

  THE COST OF SKIPPING BLOCKING

      records      all pairs        blocked by postcode (avg 50/block)
      10,000       50 million       ~2.5 million
      100,000      5 billion        ~25 million
      1,000,000    500 billion      ~250 million

  THRESHOLD TRADE-OFF

      low threshold  -> over-merges: two real customers become one
      high threshold -> under-merges: one customer counted twice
      Both inflate or deflate every per-customer metric you report.

Key Takeaways

• Write down the grain of every table and assert it in the pipeline, so duplicate keys fail loudly rather than doubling your totals.
• When duplicate keys carry different values, "which row to keep" is a business rule; sort explicitly rather than letting the drop choose.
• Normalising text before matching removes more duplicates than any similarity metric does.
• Blocking is what makes fuzzy matching tractable, and it trades recall for a feasible number of comparisons.
• Matches are transitive, so resolve them as connected components and join on the resulting entity identifier.

🧪 Practice

1. Write an assert_grain function and add it to a load routine, then demonstrate it catching an injected duplicate key.
2. Given a table with duplicate keys whose rows differ, produce three deduplicated versions using "keep latest", "keep most complete", and "aggregate", and explain which is right for an orders table.
3. Implement blocking on the first three characters of a normalised name, then report how many candidate pairs it produces versus the full cross-product on a thousand records.
4. Interview: Your customer count fell 12% after a deduplication release and marketing says the number is wrong. How do you determine whether the new number or the old one is correct? (Hint: what would you sample, and what does each error type look like in that sample?)

Type Conversion and Parsing

Data arrives as text far more often than people expect, and the conversion from text to typed values is where silent corruption concentrates. The reason is that parsers are built to be helpful: given something ambiguous, they guess, and a guess produces a plausible value rather than an error. A failed parse is a good day. A successful wrong parse is the expensive case.

The recurring pattern is that type inference is a heuristic applied to a sample. Pandas reads the first chunk of a CSV, decides a column is an integer, and then meets a row where it is not, which is why columns mysteriously become object or why integers become floats the moment a single null appears. The fix is always the same: declare the types you expect rather than accepting the ones you are given.

Four conversions deserve specific attention.

Numbers. Thousands separators, currency symbols, parentheses for negatives in accounting exports, European decimal commas, and percentage signs all defeat a naive conversion. Strip deliberately, then convert with errors="coerce" so that failures become nulls you can count rather than exceptions that stop the load. Then check the count, because coercion turns a loud failure into a quiet one.

Dates. The worst offender, because 03/04/2024 is a valid date under two different conventions and nothing in the file says which. Always pass an explicit format when you know it. When a column mixes formats, parse in passes: try each known format, and only fall back to inference for what remains. Also decide about time zones once, at the boundary: store UTC, convert for display, and never mix naive and aware timestamps in the same column.

Booleans. "true", "TRUE", "T", "1", "yes", "Y" all mean true somewhere. The trap is that bool("false") is True in Python, because any non-empty string is truthy. Map explicitly from a dictionary and leave anything unrecognised as null instead of guessing.

Integers with nulls. NumPy's integer types cannot hold a missing value, so pandas historically promoted such columns to float. That turns an ID into 1.0002e+15 and loses precision above 2^53. Pandas' nullable Int64 dtype solves this, and it is worth using by default for any integer that can be missing.

Underneath all of this sits encoding. Text is bytes plus an assumption, and if the assumption is wrong you get mojibake or a decode error. UTF-8 is the correct default; detect rather than guess when a source is unknown; and remember that a byte-order mark at the start of a file will otherwise become part of your first column's name.

import pandas as pd
import numpy as np
import io

messy = pd.DataFrame({
    "amount":  ["1,234.56", "$99.00", "(45.00)", "1.234,56", "12%", "N/A"],
    "count":   ["10", "20", None, "40", "50", "60"],
    "flag":    ["true", "FALSE", "Y", "0", "yes", "maybe"],
    "when":    ["2024-01-05", "05/03/2024", "2024-03-05T10:00:00Z",
                "March 7, 2024", "20240308", "not a date"],
    "zip":     ["01234", "02345", "03456", "04567", "05678", "06789"],
})

# --- NUMBERS: strip deliberately, then coerce, then COUNT the failures. ---
def parse_money(s: pd.Series, decimal_comma: bool = False) -> pd.Series:
    t = s.astype("string").str.strip()
    neg = t.str.match(r"^\(.*\)$").fillna(False)       # accounting negatives
    t = t.str.replace(r"^\((.*)\)$", r"\1", regex=True)
    t = t.str.replace(r"[^\d.,\-]", "", regex=True)    # drop $, %, spaces
    if decimal_comma:
        t = t.str.replace(".", "", regex=False).str.replace(",", ".", regex=False)
    else:
        t = t.str.replace(",", "", regex=False)        # thousands separator
    out = pd.to_numeric(t, errors="coerce")
    return out.where(~neg, -out)

amt = parse_money(messy["amount"])
print(pd.DataFrame({"raw": messy["amount"], "parsed": amt}).to_string(index=False))
failed = amt.isna() & messy["amount"].notna()
print(f"\nunparseable: {failed.sum()} -> {messy.loc[failed, 'amount'].tolist()}")
# NOTE: "12%" became 12.0, not 0.12, and "1.234,56" became 1.23456 under the
# comma-as-thousands assumption. Coercion did not fail; it was WRONG.
# This is why you inspect parsed values against raw, not just the null count.

# --- INTEGERS WITH NULLS: float promotion destroys identifiers. ---
naive = pd.to_numeric(messy["count"])                 # -> float64
nullable = messy["count"].astype("Int64")             # -> nullable integer
print(f"\nnaive dtype {naive.dtype}, values {naive.tolist()}")
print(f"Int64 dtype {nullable.dtype}, values {nullable.tolist()}")

big = pd.Series(["9007199254740993", "9007199254740994"])
print("\nfloat64 loses precision above 2^53:",
      pd.to_numeric(big).astype("int64").tolist())
print("Int64 preserves it:            ", big.astype("Int64").tolist())

# --- IDENTIFIERS ARE NOT NUMBERS. Keep them as strings. ---
print(f"\nzip as string: {messy['zip'].tolist()[:2]}")
print(f"zip as number: {pd.to_numeric(messy['zip']).tolist()[:2]}   <- zeros gone")

# --- DATES: parse in passes with explicit formats, infer only as fallback. ---
KNOWN_FORMATS = ["%Y-%m-%d", "%d/%m/%Y", "%Y-%m-%dT%H:%M:%S%z",
                 "%B %d, %Y", "%Y%m%d"]

def parse_dates(s: pd.Series, formats: list[str]) -> tuple[pd.Series, pd.Series]:
    """Return parsed timestamps plus the format that matched each row."""
    out = pd.Series(pd.NaT, index=s.index, dtype="datetime64[ns, UTC]")
    which = pd.Series(pd.NA, index=s.index, dtype="string")
    for fmt in formats:
        pending = out.isna() & s.notna()
        if not pending.any():
            break
        try:
            got = pd.to_datetime(s[pending], format=fmt, errors="coerce", utc=True)
        except ValueError:
            continue
        out[pending] = got
        which[pending & got.notna().reindex(s.index, fill_value=False)] = fmt
    return out, which

parsed, fmt_used = parse_dates(messy["when"], KNOWN_FORMATS)
print("\n", pd.DataFrame({"raw": messy["when"], "parsed": parsed,
                          "format": fmt_used}).to_string(index=False))
# The `format` column is the audit trail: if one format suddenly matches
# 90% of rows when it used to match 10%, the upstream source changed.

# --- THE AMBIGUITY THAT NO ERROR REPORTS ---
for dayfirst in (True, False):
    got = pd.to_datetime("03/04/2024", dayfirst=dayfirst)
    print(f"dayfirst={dayfirst!s:5s} -> {got.date()}")

# --- TIME ZONES: decide once, at the boundary. ---
naive_ts = pd.to_datetime(["2024-03-10 02:30", "2024-07-01 12:00"])
ny = naive_ts.tz_localize("America/New_York", nonexistent="shift_forward")
print("\nlocalised then to UTC:", ny.tz_convert("UTC").tolist())
# nonexistent="shift_forward" handles the hour that does not exist on a
# spring-forward date. Without it, localisation raises.

# --- BOOLEANS: map explicitly; never rely on truthiness. ---
print(f"\nthe trap: bool('false') is {bool('false')}")
TRUE_SET = {"true", "t", "yes", "y", "1"}
FALSE_SET = {"false", "f", "no", "n", "0"}

def parse_bool(s: pd.Series) -> pd.Series:
    t = s.astype("string").str.strip().str.lower()
    return pd.Series(
        np.select([t.isin(TRUE_SET), t.isin(FALSE_SET)], [True, False], default=None),
        index=s.index, dtype="boolean",     # nullable: unknown stays unknown
    )

flags = parse_bool(messy["flag"])
print(pd.DataFrame({"raw": messy["flag"], "parsed": flags}).to_string(index=False))

# --- PUT IT TOGETHER: a declared schema is the whole defence. ---
SCHEMA = {"zip": "string", "count": "Int64"}
csv = "zip,count,when\n01234,10,2024-01-05\n02345,,2024-01-06\n"
loaded = pd.read_csv(io.StringIO(csv), dtype=SCHEMA,
                     parse_dates=["when"], date_format="%Y-%m-%d",
                     na_values=["", "NA", "N/A", "null", "NULL", "-"])
print("\n", loaded.dtypes.to_dict())

# CATEGORY dtype: large memory win on repeated strings, plus it makes an
# unexpected new value visible instead of silently accepted.
country = pd.Series(["US"] * 5000 + ["CA"] * 5000)
print(f"\nobject {country.memory_usage(deep=True)/1e3:.0f} kB -> "
      f"category {country.astype('category').memory_usage(deep=True)/1e3:.0f} kB")

  FAILURE MODES, RANKED BY HOW LONG THEY GO UNNOTICED

  1. ambiguous date convention      never noticed; months are simply wrong
  2. decimal comma vs thousands     values off by 1000x in some rows only
  3. leading zeros stripped from ID joins miss, silently, for a subset
  4. int -> float promotion         precision loss above 2^53 on large IDs
  5. truthy string booleans         every row becomes True
  6. wrong encoding                 visible mojibake, or a loud decode error

  THE SAFE CONVERSION LOOP

    read as string
         |
         v
    strip known decorations ($, %, commas, parentheses)
         |
         v
    convert with errors="coerce"
         |
         v
    count new nulls  ------> non-zero? inspect the raw values, fix the rule
         |
         v
    spot-check PARSED against RAW  <- catches the wrong-but-successful parse
         |
         v
    assert dtype and range

Key Takeaways

• Type inference is a guess made from a sample, so declare dtypes, date formats, and null tokens rather than accepting the defaults.
• A successful but wrong parse costs far more than a failed one; compare parsed values against raw, not only the null count.
• Identifiers are strings, not numbers, or leading zeros and precision above 2^53 disappear.
• Use pandas' nullable Int64 and boolean dtypes so a missing value stays missing instead of promoting the column or becoming True.
• Normalise time zones to UTC at the ingestion boundary and keep a record of which date format matched each row.

🧪 Practice

1. Write a parse_money function handling thousands separators, currency symbols, and accounting parentheses, and prove it on a column containing all three.
2. Build a date column that mixes three formats, parse it in passes, and output a table of which format matched how many rows.
3. Show that reading a zero-padded identifier without a declared dtype breaks a join, then fix it with a schema.
4. Interview: A nightly load has produced correct row counts for months, but a monthly report is suddenly wrong for 3% of customers. Where do you look first? (Hint: which conversions fail for only a subset of values rather than for all of them?)

Standardizing Text and Categories

Two strings that a human reads as identical are different values to a computer, and every such difference splits a group, breaks a join, or invents a category. "United States", "united states", "USA", "U.S.", and " US " are five categories where there should be one. Standardizing means collapsing surface variation so that equal things compare equal.

The work divides into three layers, and doing them in order matters.

Unicode normalization comes first, because it is invisible. The character é can be stored as one code point or as e followed by a combining accent. The two look the same and are not equal. Normalizing to NFC makes them equal. Similarly, non-breaking spaces, zero-width characters, and typographic quotes all masquerade as ordinary characters. Normalize before you compare anything, or your comparisons are unreliable in ways that are extremely hard to debug from a printed table.

Mechanical cleaning comes second: trim whitespace, collapse internal runs of whitespace, and case-fold. Use casefold rather than lower if you handle non-English text, because it handles cases like the German sharp s that lower does not. Decide whether to strip punctuation and accents based on the domain: removing accents merges "Muller" and "Müller", which is right for fuzzy matching and wrong for a mailing label.

Semantic mapping comes third and is the only layer that requires knowledge. This is where "USA" becomes "US". The right structure is an explicit dictionary, version-controlled, applied deterministically, with anything unmatched routed to a review list rather than silently dropped or bucketed as "Other". Wherever a standard exists, use it: ISO 3166 for countries, ISO 4217 for currencies, ISO 639 for languages. A canonical external standard is a mapping you do not have to maintain or defend.

For the messy remainder, fuzzy grouping helps. Computing string similarity within a column and clustering near-identical values surfaces the "Acme Inc" / "Acme Inc." pairs quickly. Treat its output as a proposal for a human to approve, not an automatic merge, because the same technique that merges typos will happily merge "Springfield, IL" and "Springfield, MA".

Finally, once a column is standardized, pin the allowed values. Convert to a categorical dtype with an explicit list of categories, so that a new value arriving next month appears as a null or an error rather than enlarging your category set without notice. That turns a silent data drift into a failing check.

import pandas as pd
import unicodedata
import re
from difflib import SequenceMatcher

raw = pd.Series([
    "United States", "united states", "USA", "U.S.", " US ",
    "  Canada", "CANADA", "canada\u00a0", "Deutschland", "Germany",
    "Café", "Cafe\u0301", "", None, "acme inc", "Acme Inc.",
], name="value")

# --- LAYER 1: Unicode normalization, the invisible problem. ---
a, b = "Café", "Cafe\u0301"                   # composed vs decomposed
print(f"look identical: {a == b}   lengths: {len(a)} vs {len(b)}")
nfc = unicodedata.normalize
print(f"after NFC:      {nfc('NFC', a) == nfc('NFC', b)}")

INVISIBLES = {
    "\u00a0": " ",    # non-breaking space
    "\u200b": "",     # zero-width space
    "\ufeff": "",     # byte-order mark
    "\u2019": "'",    # right single quote
    "\u2014": "-",    # em dash
}

def normalise_text(s: pd.Series, strip_accents: bool = False) -> pd.Series:
    t = s.astype("string")
    t = t.map(lambda v: unicodedata.normalize("NFC", v) if pd.notna(v) else v)
    for bad, good in INVISIBLES.items():
        t = t.str.replace(bad, good, regex=False)
    if strip_accents:
        # NFKD splits a letter from its accent, then we drop the accents.
        t = t.map(lambda v: "".join(
            c for c in unicodedata.normalize("NFKD", v)
            if not unicodedata.combining(c)) if pd.notna(v) else v)
    # LAYER 2: mechanical cleaning. casefold beats lower on non-English text.
    t = t.str.strip().str.replace(r"\s+", " ", regex=True).str.casefold()
    return t.replace("", pd.NA)               # empty string is a null

key = normalise_text(raw, strip_accents=True)
print("\n", pd.DataFrame({"raw": raw, "key": key}).to_string())
print(f"\ndistinct before {raw.nunique()} -> after {key.nunique()}")

# --- LAYER 3: semantic mapping against a standard. ---
# Version-control this dictionary. It is business logic, not a constant.
COUNTRY_MAP = {
    "united states": "US", "usa": "US", "us": "US", "united states of america": "US",
    "canada": "CA", "ca": "CA",
    "germany": "DE", "deutschland": "DE", "de": "DE",
}

mapped = key.map(COUNTRY_MAP)
unmatched = key[mapped.isna() & key.notna()]
print("\nmapped:\n", pd.DataFrame({"key": key, "iso": mapped}).dropna().to_string())
print("\nTO REVIEW (never silently bucket these):", unmatched.unique().tolist())
# Routing the remainder to a review queue is the difference between a
# mapping you can trust and an "Other" bucket that grows forever.

# --- FUZZY GROUPING: propose merges, do not apply them. ---
def propose_merges(values: pd.Series, threshold: float = 0.88) -> pd.DataFrame:
    uniq = sorted(values.dropna().unique())
    out = []
    for i, a in enumerate(uniq):
        for b in uniq[i + 1:]:
            r = SequenceMatcher(None, a, b).ratio()
            if r >= threshold:
                out.append({"a": a, "b": b, "ratio": round(r, 3)})
    return pd.DataFrame(out).sort_values("ratio", ascending=False)

companies = normalise_text(pd.Series(
    ["Acme Inc", "Acme Inc.", "ACME INC", "Acme Corp",
     "Springfield IL", "Springfield MA", "Zenith Ltd"]))
print("\nproposed merges for human review:")
print(propose_merges(companies).to_string(index=False))
# Note the Springfield pair: high similarity, and merging them would be a
# factual error. This is exactly why the output is a proposal.

# --- PIN THE ALLOWED VALUES so drift fails a check instead of passing. ---
ALLOWED = ["US", "CA", "DE", "GB", "FR"]
col = pd.Series(["US", "CA", "DE", "MX"]).astype(
    pd.CategoricalDtype(categories=ALLOWED))
print("\npinned categories:\n", col.to_string())
print(f"unexpected values became null: {col.isna().sum()}")

def assert_allowed(s: pd.Series, allowed: list[str], name: str) -> None:
    unexpected = set(s.dropna().unique()) - set(allowed)
    if unexpected:
        raise ValueError(f"{name}: unexpected categories {sorted(unexpected)}")

try:
    assert_allowed(pd.Series(["US", "CA", "MX"]), ALLOWED, "country")
except ValueError as e:
    print("caught:", e)

# --- RARE-CATEGORY COLLAPSING: a threshold, applied from TRAIN counts. ---
cats = pd.Series(["a"] * 500 + ["b"] * 300 + ["c"] * 8 + ["d"] * 3 + ["e"] * 1)
freq = cats.value_counts(normalize=True)
keep = freq[freq >= 0.01].index                  # 1% floor, fitted on train
collapsed = cats.where(cats.isin(keep), "__rare__")
print("\nbefore:", cats.nunique(), "after:", collapsed.nunique())
print(collapsed.value_counts().to_string())
# Store `keep` as a fitted parameter. Recomputing it on test data leaks.

  THE THREE LAYERS, IN ORDER

  1. UNICODE           NFC normalize, strip invisibles      no knowledge needed
     "Cafe\u0301" -> "Café",  "US\u00a0" -> "US "

  2. MECHANICAL        trim, collapse spaces, casefold      no knowledge needed
     " United States " -> "united states"

  3. SEMANTIC          map to a canonical value             DOMAIN KNOWLEDGE
     "united states" -> "US"  (ISO 3166)

  Skipping layer 1 makes layers 2 and 3 unreliable in ways a printed
  table cannot show you.

  WHERE TO GET CANONICAL VALUES

    countries   ISO 3166-1 alpha-2      currencies  ISO 4217
    languages   ISO 639-1               time zones  IANA tz database
    industries  NAICS / SIC             units       UCUM

  Using a standard means the mapping is someone else's job to maintain.

Key Takeaways

• Normalize Unicode first, because composed and decomposed characters look identical while comparing unequal.
• Separate mechanical cleaning, which needs no knowledge, from semantic mapping, which needs domain expertise and belongs in a version- controlled dictionary.
• Map to an external standard such as ISO 3166 whenever one exists.
• Fuzzy grouping proposes merges for review; applied automatically it will merge genuinely distinct values.
• Pin allowed categories so a new value fails a check instead of quietly expanding the category set, and fit rare-category thresholds on training data.

🧪 Practice

1. Build a series containing a composed and a decomposed version of the same accented word, show they are unequal, and fix it with NFC normalization.
2. Write a standardization function with all three layers and report the distinct-value count after each layer on a real categorical column.
3. Implement a rare-category collapse at a 1% threshold as a fitted transform, then apply it to a held-out set without recomputing counts.
4. Interview: A dashboard's country breakdown shows both "US" and "USA" with separate totals. Where in the pipeline would you fix it, and what do you do about the historical rows? (Hint: fixing forward and fixing backward are different projects with different risks.)

<a id="62-data-transformation"></a>

6.2 Data Transformation

Transformation is the vocabulary of data manipulation: a small set of operations that compose into any reshaping you need, and whose failure modes are mostly about accidentally changing the grain of a table without realising it.

Filtering, Sorting, Slicing

These three operations are the ones everybody learns first and the ones that produce the most subtle bugs, precisely because they look too simple to get wrong. The bugs come from three places: nulls do not behave like values in comparisons, sorting is not stable unless you ask, and selecting by position and selecting by label are different operations that sometimes give the same answer.

Filtering is selecting rows by a condition. The critical fact is that a comparison against a null is neither true nor false, and so a filter silently drops null rows from both sides. df[df.age > 30] excludes nulls, and so does df[df.age <= 30]. Their row counts do not add up to the total, and whoever is reading your output will not know why. If you want the nulls, you must say so.

The second filtering trap is chained comparisons. In Python, a < x < b on a Series does not work as written because the ampersand and pipe operators bind more tightly than comparison, so every compound condition needs parentheses. between avoids the issue entirely and is clearer, but be aware that it is inclusive on both ends by default.

Sorting matters more than it seems because downstream operations depend on order: forward fills, rolling windows, drop_duplicates(keep="first"), and head. The trap is ties. Two rows with the same sort key can come out in either order, and the order can change between library versions or when the data changes upstream. That makes results non-reproducible in a way that is maddening to debug. The fix is to always sort by enough columns to break every tie, and to use a stable sort algorithm when you cannot.

Null placement in sorts is also a choice, not a default worth trusting. Say where you want them.

Slicing distinguishes label-based and position-based access. loc takes labels and includes the endpoint; iloc takes positions and excludes it. When your index happens to be 0..n-1 they coincide, which is exactly how the bug hides until the index is a date or has been filtered.

Finally, a practical point about copies and views. Filtering usually returns a new object, but assigning into the result of a chained selection may write into a temporary and vanish, or may raise a warning about setting on a copy. Assign through a single loc call, or call copy() explicitly when you intend an independent frame.

import pandas as pd
import numpy as np

df = pd.DataFrame({
    "user":   ["a", "b", "c", "d", "e", "f"],
    "age":    [25, 34, np.nan, 41, 34, 19],
    "plan":   ["free", "pro", "pro", "free", "pro", None],
    "spend":  [0.0, 120.0, 80.0, 15.0, 120.0, 5.0],
    "joined": pd.to_datetime(["2024-01-05", "2024-02-11", "2024-02-11",
                              "2024-03-02", "2024-03-02", "2024-04-19"]),
})

# --- NULLS VANISH FROM BOTH SIDES OF A FILTER ---
over = df[df["age"] > 30]
under = df[df["age"] <= 30]
print(f"total {len(df)}, over30 {len(over)}, under30 {len(under)},"
      f" sum {len(over) + len(under)}   <- does not add up")
print("the missing row:\n", df[df["age"].isna()].to_string(index=False))

# Be explicit about which side the nulls belong on.
over_or_unknown = df[(df["age"] > 30) | df["age"].isna()]
print(f"\nover30 including unknown age: {len(over_or_unknown)}")

# --- COMPOUND CONDITIONS NEED PARENTHESES ---
# WRONG: `df["age"] > 20 & df["age"] < 40` binds as `20 & df["age"]`.
mask = (df["age"] > 20) & (df["age"] < 40)
print(f"\nparenthesised mask selects {mask.sum()} rows")
print(f"between(20, 40) selects {df['age'].between(20, 40).sum()} rows"
      f"   <- inclusive on BOTH ends")
print(f"between exclusive: {df['age'].between(20, 40, inclusive='neither').sum()}")

# `isin` for set membership; `~` to negate. Note nulls are never `isin`.
print(f"\npro or free: {df['plan'].isin(['pro', 'free']).sum()} of {len(df)}")
print("query() reads well for long conditions:\n",
      df.query("age > 20 and spend < 100").to_string(index=False))

# --- SORTING: TIES ARE NON-DETERMINISTIC UNLESS YOU BREAK THEM ---
ambiguous = df.sort_values("age")                  # rows b and e tie at 34
deterministic = df.sort_values(["age", "user"], kind="mergesort")  # stable
print("\nambiguous tie order:", ambiguous["user"].tolist())
print("deterministic order: ", deterministic["user"].tolist())
# kind="mergesort" is stable: equal keys keep their original relative order.
# The default quicksort is not, so the same code can reorder on new data.

# Null placement is a decision.
print("\nnulls last: ", df.sort_values("age", na_position="last")["user"].tolist())
print("nulls first:", df.sort_values("age", na_position="first")["user"].tolist())

# Multi-key with mixed direction: newest first, then cheapest.
print("\n", df.sort_values(["joined", "spend"], ascending=[False, True])
        [["user", "joined", "spend"]].to_string(index=False))

# WHY ORDER MATTERS DOWNSTREAM: the same dedupe, two answers.
for direction in (True, False):
    kept = (df.sort_values("joined", ascending=direction)
              .drop_duplicates("plan", keep="first")["user"].tolist())
    print(f"dedupe after ascending={direction}: kept {kept}")

# --- LOC VS ILOC: labels include the endpoint, positions do not ---
indexed = df.set_index("user")
print("\nloc 'b':'d' (inclusive):", indexed.loc["b":"d"].index.tolist())
print("iloc 1:3   (exclusive):  ", indexed.iloc[1:3].index.tolist())

filtered = df[df["spend"] > 10]                   # index is now 1,2,3,4
print("\nindex after filtering:", filtered.index.tolist())
print("iloc[0] ->", filtered.iloc[0]["user"], " (position)")
print("loc[1]  ->", filtered.loc[1]["user"], " (label)")
# reset_index(drop=True) after filtering, or these two keep diverging.

# --- BOOLEAN MASKS ARE COMPOSABLE AND SELF-DOCUMENTING ---
is_paying = df["spend"] > 0
is_recent = df["joined"] >= "2024-03-01"
print(f"\npaying {is_paying.sum()}, recent {is_recent.sum()},"
      f" both {(is_paying & is_recent).sum()}")
# Naming masks turns an unreadable one-liner into a readable filter.

# --- ASSIGNMENT: one loc call, not a chained selection ---
safe = df.copy()
safe.loc[safe["plan"].isna(), "plan"] = "unknown"   # single indexing operation
print("\n", safe[["user", "plan"]].to_string(index=False))

  FILTER TRUTH TABLE WITH NULLS

    value    x > 30    x <= 30    x.isna()
    41       True      False      False
    19       False     True       False
    NaN      False     False      True        <- excluded by BOTH comparisons

    Any partition built from comparisons alone loses the null rows.
    Row counts that do not sum to the total are the symptom.

  LOC VS ILOC

    index:   'a'  'b'  'c'  'd'  'e'
    position  0    1    2    3    4

    loc['b':'d']  -> b, c, d      (labels, endpoint INCLUDED)
    iloc[1:3]     -> b, c         (positions, endpoint EXCLUDED)

Key Takeaways

• A comparison against a null is false, so filters silently drop null rows from both sides of a partition; decide explicitly where they belong.
• Parenthesise every compound mask, and remember between is inclusive by default.
• Break every sort tie with an additional key, and use a stable sort, or results change when the data does.
• loc is label-based and endpoint-inclusive; iloc is position-based and endpoint-exclusive, and they coincide only on a default index.
• Name boolean masks as variables; it documents intent and makes composition readable.

🧪 Practice

1. Build a column with nulls and show that the counts of x > k and x <= k do not sum to the row count, then write a partition that does.
2. Create a table with tied sort keys and demonstrate two different drop_duplicates outcomes depending on sort direction.
3. Filter a frame, then show that iloc[0] and loc[0] return different rows, and explain when this bites in production code.
4. Interview: A colleague reports that a filtered report's totals do not match the unfiltered table's total. What are the two most likely causes? (Hint: one involves nulls, the other involves the grain.)

Grouping and Aggregation

Grouping is the single most useful operation in analysis, and it works by a pattern worth naming explicitly: split, apply, combine. Split the rows into groups by some key, apply a function to each group, and combine the results into one output. Once you see the pattern, SQL's GROUP BY, pandas' groupby, and Spark's aggregations are all the same thing with different syntax.

The thing that turns grouping from a formula into a tool is understanding the three shapes the "apply" step can produce, because choosing the wrong one is the most common source of confusion.

Aggregate reduces each group to a single row. Ten thousand orders become fifty rows, one per country. This changes the grain of the table, which is the point.

Transform returns a value for every original row, broadcast from its group. The group mean attached to each row, so you can compute a deviation from the group's own average. The grain is unchanged, which is what makes this composable with the original frame.

Filter keeps or discards whole groups based on a group-level condition. Keep only countries with more than a hundred orders. Rows survive or not together.

The analogy: aggregate is writing one summary line per class on a report. Transform is writing each student's class average next to their own mark. Filter is deciding which classes appear in the report at all.

Beyond the shapes, several practical points separate reliable aggregation from misleading aggregation.

Nulls are skipped by default, which means mean divides by the count of non-null values, not by the group size. That is usually what you want and occasionally a serious misrepresentation, so report the count alongside any mean. Similarly count counts non-nulls while size counts rows, and the difference between them is exactly the null count.

Empty groups disappear. If a category has no rows in this period, it simply does not appear in the output, and a chart built from that output shows no gap. Reindexing against the full set of expected keys turns a missing row into an explicit zero.

Averages of averages are wrong. The mean of per-country conversion rates is not the overall conversion rate unless every country has the same number of users. Aggregate the numerator and denominator separately, then divide. This mistake, sometimes met as Simpson's paradox when it reverses a conclusion, appears constantly in real reporting.

Finally, named aggregation is the syntax worth defaulting to. It states the output column name, the input column, and the function on one line, which means the result has readable flat column names instead of a multi-level index you then have to flatten.

import pandas as pd
import numpy as np

rng = np.random.default_rng(1)
orders = pd.DataFrame({
    "country": rng.choice(["US", "CA", "DE"], 400, p=[0.6, 0.3, 0.1]),
    "plan":    rng.choice(["free", "pro"], 400),
    "amount":  rng.gamma(2, 40, 400).round(2),
    "items":   rng.integers(1, 6, 400),
})
orders.loc[rng.random(400) < 0.08, "amount"] = np.nan     # some missing values

# --- AGGREGATE: one row per group. Named aggregation gives flat names. ---
summary = orders.groupby("country", as_index=False).agg(
    n_rows=("amount", "size"),          # ALL rows in the group
    n_amounts=("amount", "count"),      # NON-NULL rows only
    revenue=("amount", "sum"),
    avg_amount=("amount", "mean"),
    p90_amount=("amount", lambda s: s.quantile(0.90)),
    distinct_plans=("plan", "nunique"),
)
summary["n_missing"] = summary["n_rows"] - summary["n_amounts"]
print(summary.round(2).to_string(index=False))
# size - count IS the null count. Reporting a mean without the count that
# produced it hides how much data the mean rests on.

# --- TRANSFORM: same shape as the input, broadcast from the group. ---
orders["country_avg"] = orders.groupby("country")["amount"].transform("mean")
orders["vs_country"] = (orders["amount"] - orders["country_avg"]).round(2)
# Group-relative z-score: a standard feature-engineering move.
orders["z_in_country"] = orders.groupby("country")["amount"].transform(
    lambda s: (s - s.mean()) / s.std(ddof=0)).round(2)
print("\n", orders.head(5)[
    ["country", "amount", "country_avg", "vs_country", "z_in_country"]
].round(2).to_string(index=False))

# Rank within group, another transform that is hard to write any other way.
orders["rank_in_country"] = orders.groupby("country")["amount"].rank(
    method="dense", ascending=False)

# --- FILTER: keep or drop entire groups. ---
big = orders.groupby("country").filter(lambda g: len(g) >= 50)
print(f"\nkept {big['country'].nunique()} of {orders['country'].nunique()}"
      f" countries, {len(big)} of {len(orders)} rows")

# --- THE AVERAGE-OF-AVERAGES ERROR ---
per_country = orders.groupby("country").agg(
    rev=("amount", "sum"), n=("amount", "count"))
per_country["avg"] = per_country["rev"] / per_country["n"]

wrong = per_country["avg"].mean()                       # unweighted
right = per_country["rev"].sum() / per_country["n"].sum()  # pooled
print(f"\nmean of group means {wrong:.2f}  vs  overall mean {right:.2f}")
print("weighted by group size:",
      round(np.average(per_country["avg"], weights=per_country["n"]), 2))
# RULE: aggregate the numerator and the denominator, THEN divide.

# --- EMPTY GROUPS DISAPPEAR unless you reindex. ---
EXPECTED = ["US", "CA", "DE", "FR", "JP"]
counts = orders.groupby("country").size()
print("\nas produced:\n", counts.to_string())
print("\nreindexed against the expected set:\n",
      counts.reindex(EXPECTED, fill_value=0).to_string())
# The second version shows FR and JP as zero. The first shows nothing, and
# a chart built on it silently omits them.

# Same issue with two keys: observed=False materialises absent combinations
# for categorical keys.
cat = orders.assign(
    country=orders["country"].astype("category"),
    plan=orders["plan"].astype("category"))
full_grid = cat.groupby(["country", "plan"], observed=False).size()
present = len(cat.groupby(["country", "plan"], observed=True).size())
print(f"\ncombinations present {present}"
      f", full grid {len(full_grid)}")

# --- MULTI-KEY GROUPING and flattening a MultiIndex result. ---
multi = orders.groupby(["country", "plan"]).agg({"amount": ["sum", "mean"],
                                                 "items": "sum"})
multi.columns = ["_".join(c) for c in multi.columns]     # flatten
print("\n", multi.round(2).to_string())

# --- SIMPSON'S PARADOX: the aggregate reverses the subgroups. ---
trial = pd.DataFrame({
    "group":    ["A"] * 2 + ["B"] * 2,
    "severity": ["mild", "severe"] * 2,
    "treated":  [100, 900, 900, 100],
    "cured":    [90, 270, 800, 25],
})
trial["rate"] = trial["cured"] / trial["treated"]
print("\nby severity:\n", trial.to_string(index=False))
overall = trial.groupby("group").apply(
    lambda g: g["cured"].sum() / g["treated"].sum(), include_groups=False)
print("\noverall rate:\n", overall.round(3).to_string())
# A beats B in BOTH severity strata and loses overall, because the groups
# have different severity mixes. Never aggregate across a confounded split.

  THE THREE SHAPES OF SPLIT-APPLY-COMBINE

  input: 6 rows, 2 groups                   AGGREGATE (grain changes)
  +-----+--------+                          +-----+-------+
  | grp | value  |                          | grp | mean  |
  +-----+--------+   groupby("grp")         +-----+-------+
  |  A  |   10   | ------------------>      |  A  |  20   |
  |  A  |   30   |     .agg("mean")         |  B  |  50   |
  |  A  |   20   |                          +-----+-------+
  |  B  |   40   |
  |  B  |   60   |   .transform("mean")     TRANSFORM (grain preserved)
  |  B  |   50   | ------------------>      +-----+-------+------+
  +-----+--------+                          | grp | value | mean |
                                            |  A  |  10   |  20  |
                     .filter(len >= 4)      |  A  |  30   |  20  |
                   ------------------>      |  ...       |      |
                     FILTER: both groups    +-----+-------+------+
                     have 3 rows, so the
                     output is EMPTY

  COUNT VS SIZE

    size()  = rows in the group          (nulls included)
    count() = non-null values            (nulls excluded)
    size - count = the null count, and the gap your mean ignored

Key Takeaways

• Aggregate changes the grain, transform preserves it, and filter keeps or drops whole groups; picking the wrong one is the usual confusion.
• Aggregations skip nulls, so report size alongside any mean to show how much data it rests on.
• Averaging group averages is wrong unless group sizes are equal; sum the numerator and denominator, then divide.
• Groups with no rows vanish from the output, so reindex against the full expected key set to make zeros explicit.
• Prefer named aggregation for flat, readable output columns.

🧪 Practice

1. For a grouped table, compute the mean of group means and the pooled mean, and construct group sizes that make them differ by more than 20%.
2. Add a column holding each row's deviation from its group's median using transform, then use it to find the most unusual row per group.
3. Produce a monthly count by category that includes zero rows for categories absent in a given month.
4. Interview: A weekly report shows conversion rate rising in every channel while the overall rate falls. Explain how that is possible and what you would present instead. (Hint: what changed about the mix?)

Joins and Merges

A join combines two tables by matching rows on a key. It is the operation that makes relational data useful and the operation most likely to silently change your row count, which is why the single most valuable habit in this whole chapter is: check the row count before and after every join, and know what you expected it to be.

The join types differ in which non-matching rows survive.

Join	Keeps	Use when
inner	rows matching in both	you need both sides populated
left	all left rows	enriching a table you must not shrink
right	all right rows	rarely; write it as a left join instead
outer (full)	all rows from both	reconciling two sources
cross	every pair	building a complete grid deliberately
anti	left rows with no match	finding what is missing

Join type is only half the story. The other half is cardinality, the relationship between key multiplicities, and this is where the damage happens.

A one-to-one join leaves row counts alone. A many-to-one join, the normal case of looking up an attribute, also leaves the left row count alone, because each left row finds exactly one match. A many-to-many join multiplies: three rows on the left sharing a key and four on the right produce twelve rows. If you expected an enrichment you now have duplicated revenue, and because every row looks individually plausible, nothing alerts you except the total.

Pandas will tell you if you ask. The validate parameter asserts the cardinality and raises if the data violates it. Using it converts a silent fan-out into a stack trace, which is precisely the trade you want.

Several more traps are worth internalising.

Key type mismatches produce zero matches with no error. An integer 101 does not equal the string "101", and a key with trailing whitespace does not equal one without. A join returning suspiciously few rows is almost always this.

Nulls do not match nulls. Two rows with a null key are not joined together, which is correct behaviour and surprising the first time.

Indicator columns are the diagnostic tool. Merging with indicator=True adds a column saying whether each row came from the left, the right, or both, and the counts of those three values tell you immediately what your join did.

For time series there is a distinct operation: as-of joins, which match each left row to the most recent right row at or before its timestamp. This is how you attach a price to a trade or a configuration to an event, and doing it with a regular join plus a filter is both slower and easier to get wrong. Crucially, an as-of join with the direction reversed is a leakage bug: matching to the next value attaches future information.

import pandas as pd
import numpy as np

orders = pd.DataFrame({
    "order_id": [1, 2, 3, 4, 5],
    "user_id":  [10, 11, 12, 11, 99],
    "amount":   [50.0, 20.0, 75.0, 30.0, 10.0],
})
users = pd.DataFrame({
    "user_id": [10, 11, 12, 13],
    "country": ["US", "CA", "DE", "FR"],
})

# --- THE HABIT: assert what you expect, every time. ---
def safe_merge(left, right, expect_rows=None, **kw):
    out = left.merge(right, indicator=True, **kw)
    counts = out["_merge"].value_counts().to_dict()
    print(f"  {len(left)} x {len(right)} -> {len(out)} rows   {counts}")
    if expect_rows is not None and len(out) != expect_rows:
        raise AssertionError(f"expected {expect_rows} rows, got {len(out)}")
    return out.drop(columns="_merge")

print("inner:")
inner = safe_merge(orders, users, on="user_id", how="inner")
print("left:")
left = safe_merge(orders, users, on="user_id", how="left", expect_rows=len(orders))
print("outer:")
outer = safe_merge(orders, users, on="user_id", how="outer")

# ANTI-JOIN: what is on the left and missing on the right. This is the
# question "which orders reference a user we do not have?" and it is the
# check that finds broken referential integrity.
merged = orders.merge(users, on="user_id", how="left", indicator=True)
orphans = merged[merged["_merge"] == "left_only"]
print("\norphan orders:\n", orphans[["order_id", "user_id"]].to_string(index=False))

# --- THE FAN-OUT: many-to-many silently multiplies rows and money. ---
addresses = pd.DataFrame({          # a user can have several addresses
    "user_id": [10, 10, 11, 11, 11],
    "kind":    ["home", "work", "home", "work", "billing"],
})
fanned = orders.merge(addresses, on="user_id", how="left")
print(f"\nrows {len(orders)} -> {len(fanned)}")
print(f"revenue {orders['amount'].sum():.2f} -> {fanned['amount'].sum():.2f}"
      f"   <- INFLATED, and every row looks fine")

# The fix is to decide the grain BEFORE joining: collapse the right side.
one_per_user = (addresses[addresses["kind"] == "home"]
                .drop_duplicates("user_id"))
fixed = orders.merge(one_per_user, on="user_id", how="left")
print(f"after collapsing right side: {len(fixed)} rows,"
      f" revenue {fixed['amount'].sum():.2f}")

# `validate` turns the silent fan-out into an exception.
try:
    orders.merge(addresses, on="user_id", how="left", validate="many_to_one")
except Exception as e:
    print("\nvalidate caught it:", type(e).__name__, "-", str(e)[:60])

# --- KEY TYPE MISMATCH: zero matches, no error. ---
a = pd.DataFrame({"id": [1, 2, 3], "x": ["a", "b", "c"]})
b = pd.DataFrame({"id": ["1", "2", "3"], "y": ["p", "q", "r"]})
print(f"\nint vs string key -> {len(a.merge(b, on='id'))} matches")
print(f"after aligning dtypes -> "
      f"{len(a.assign(id=a['id'].astype(str)).merge(b, on='id'))} matches")

# Whitespace is the same class of bug and even less visible.
c = pd.DataFrame({"k": ["us ", "ca"], "v": [1, 2]})
d = pd.DataFrame({"k": ["us", "ca"], "w": [3, 4]})
print(f"with whitespace -> {len(c.merge(d, on='k'))} matches")
print(f"after strip     -> "
      f"{len(c.assign(k=c['k'].str.strip()).merge(d, on='k'))} matches")

# --- NULL KEYS DO NOT MATCH EACH OTHER ---
e = pd.DataFrame({"k": [1, None], "v": ["x", "y"]})
f = pd.DataFrame({"k": [1, None], "w": ["p", "q"]})
print(f"\nnull key matches: {len(e.merge(f, on='k'))} (not 2)")

# --- COLUMN NAME COLLISIONS: always name the suffixes. ---
g = pd.DataFrame({"id": [1], "value": [10], "updated": ["2024-01-01"]})
h = pd.DataFrame({"id": [1], "value": [99], "updated": ["2024-02-01"]})
collided = g.merge(h, on="id", suffixes=("_order", "_invoice"))
print("\n", collided.to_string(index=False))

# --- AS-OF JOIN: the correct tool for "the latest value at that time". ---
trades = pd.DataFrame({
    "time":  pd.to_datetime(["2024-01-01 09:31", "2024-01-01 09:35",
                             "2024-01-01 09:47"]),
    "shares": [100, 250, 75],
}).sort_values("time")
prices = pd.DataFrame({
    "time":  pd.to_datetime(["2024-01-01 09:30", "2024-01-01 09:34",
                             "2024-01-01 09:40", "2024-01-01 09:50"]),
    "price": [10.0, 10.5, 11.0, 12.0],
}).sort_values("time")

# direction="backward" is the only safe default: it can only use the past.
asof = pd.merge_asof(trades, prices, on="time", direction="backward")
print("\nas-of backward (safe):\n", asof.to_string(index=False))

leaky = pd.merge_asof(trades, prices, on="time", direction="forward")
print("\nas-of forward (LEAKS the future):\n", leaky.to_string(index=False))

# tolerance bounds how stale a match may be; beyond it, the value is null,
# which is more honest than silently using an hour-old price.
bounded = pd.merge_asof(trades, prices, on="time", direction="backward",
                        tolerance=pd.Timedelta("2min"))
print("\nwith a 2-minute tolerance:\n", bounded.to_string(index=False))

  CARDINALITY DECIDES THE ROW COUNT, NOT THE JOIN TYPE

  left keys    right keys    result            validate=
  ---------------------------------------------------------------
  unique       unique        len(left) rows    "one_to_one"
  repeated     unique        len(left) rows    "many_to_one"   <- lookups
  unique       repeated      MORE than left    "one_to_many"
  repeated     repeated      product per key   (no safe option)

  THE FAN-OUT, DRAWN

    left                right                 result
    k=A amount=50       k=A kind=home         k=A 50 home
                        k=A kind=work         k=A 50 work     <- 50 counted twice
    sum(amount)=50                            sum(amount)=100

  DIAGNOSING A JOIN WITH indicator=True

    both        matched as intended
    left_only   right side is missing keys -> referential integrity gap
    right_only  left side is missing keys  -> filtered upstream? wrong grain?

  ZERO MATCHES CHECKLIST
    dtype mismatch (int vs str) | whitespace | case | leading zeros stripped
    | null keys | key is actually two columns

Key Takeaways

• Check row counts before and after every join against an expectation you stated in advance.
• The join type decides which non-matches survive; cardinality decides the row count, and many-to-many fan-out inflates every sum.
• Use validate to turn an unexpected fan-out into an exception.
• Zero matches almost always means a key dtype, whitespace, or case mismatch rather than genuinely absent data.
• Use an as-of join for "latest value at that time", always backward, and set a tolerance so stale matches become null.

🧪 Practice

1. Build two frames whose join fans out, show that the sum of a numeric column inflates, then fix it by collapsing the right side first.
2. Merge with indicator=True on data with orphan keys on both sides and interpret all three indicator counts.
3. Use merge_asof to attach the most recent price to each trade, then repeat with direction="forward" and explain why the second version would inflate a backtest.
4. Interview: A join that should enrich 1 million orders returns 1.4 million rows. Walk through your diagnosis. (Hint: which side has duplicate keys, and what is the intended grain of the output?)

Reshaping: Wide and Long Formats

The same data can be laid out in two fundamentally different ways, and moving between them is one of the most common operations in practice. In wide format, each variable has its own column: one row per entity, with jan, feb, mar as separate columns. In long format, one column holds the variable name and another holds the value: one row per entity-month pair.

Neither is correct in general. They serve different purposes.

Wide format is for reading and reporting. A human comparing months scans across a row. A spreadsheet pivot table is wide. Correlation matrices and many statistical routines expect wide input.

Long format is for computing and storing. Adding a new month means adding rows, not altering a schema. Grouping, filtering, and plotting libraries all prefer it, because "the variable" is a value you can filter on rather than a column name you have to enumerate. Databases store long for the same reason.

The mental model: wide is a printed table, long is a list of facts. A fact is (entity, variable, value), and any wide table is just a pivoted view of the underlying facts.

The operations, by name:

• melt or unpivot goes wide to long. You name the columns that stay fixed, the identifiers, and the rest collapse into a name column and a value column.
• pivot goes long to wide. It requires that the combination of index and column keys be unique, and it raises if it is not, which is a useful safety check.
• pivot_table goes long to wide with aggregation, so duplicates are summarised rather than rejected. Useful, but be aware that it silently aggregates where pivot would have warned you about unexpected duplicates.
• stack and unstack move levels between the index and the columns, which is the same reshape expressed on a MultiIndex.

Two practical hazards. First, pivoting produces columns whose names come from data, which means a new value upstream changes your schema. Code that indexes those columns by name breaks in production, so either pin the expected set or work in long format until the last step. Second, pivoting creates nulls for combinations that do not exist, and those nulls mean "no observation", not "zero". Filling them with zero is a claim, sometimes right and sometimes badly wrong.

import pandas as pd
import numpy as np

wide = pd.DataFrame({
    "store": ["A", "B", "C"],
    "region": ["north", "north", "south"],
    "jan": [100, 150, 90],
    "feb": [110, 140, 95],
    "mar": [120, 160, np.nan],
})
print("WIDE (good for reading):\n", wide.to_string(index=False))

# --- WIDE -> LONG with melt ---
long = wide.melt(
    id_vars=["store", "region"],      # columns that stay as identifiers
    value_vars=["jan", "feb", "mar"],  # columns to collapse
    var_name="month",
    value_name="sales",
)
print("\nLONG (good for computing):\n", long.to_string(index=False))

# Why long is easier to compute on: one expression, no column enumeration.
print("\ngrouped in long form:\n",
      long.groupby("month", as_index=False)["sales"].agg(["sum", "count"])
      .to_string(index=False))
# In wide form the same result needs the month columns listed by name, and
# adding April means editing the code.

# --- LONG -> WIDE with pivot; it REJECTS duplicate index/column pairs. ---
back = long.pivot(index=["store", "region"], columns="month",
                  values="sales").reset_index()
back.columns.name = None
print("\npivoted back:\n", back.to_string(index=False))

dupes = pd.concat([long, long.iloc[[0]]])       # inject a duplicate fact
try:
    dupes.pivot(index="store", columns="month", values="sales")
except ValueError as e:
    print("\npivot refuses duplicates:", str(e)[:58])

# pivot_table AGGREGATES instead of refusing. Convenient, and it hides the
# duplicate you might have wanted to know about.
pt = dupes.pivot_table(index="store", columns="month", values="sales",
                       aggfunc="sum", margins=True, margins_name="total")
print("\npivot_table with margins:\n", pt.to_string())

# --- NULLS FROM PIVOTING MEAN "NO ROW", NOT "ZERO" ---
sparse = pd.DataFrame({
    "store": ["A", "A", "B"],
    "month": ["jan", "feb", "jan"],
    "sales": [100, 110, 150],
})
grid = sparse.pivot(index="store", columns="month", values="sales")
print("\nsparse pivot:\n", grid.to_string())
print("\nfilled with 0 (a CLAIM that B sold nothing in feb):\n",
      grid.fillna(0).to_string())
# If the store did not exist in February, 0 is wrong and NaN is right.
# If it existed and sold nothing, 0 is right. The data cannot tell you.

# --- STACK / UNSTACK: the same reshape on a MultiIndex ---
mi = long.dropna().set_index(["store", "month"])["sales"]
print("\nunstack moves an index level to columns:\n",
      mi.unstack("month").to_string())
print("\nstack moves it back:\n", mi.unstack("month").stack().head().to_string())

# --- WIDE COLUMN NAMES THAT ENCODE TWO VARIABLES ---
messy = pd.DataFrame({
    "id": [1, 2],
    "sales_2023": [10, 20], "sales_2024": [11, 21],
    "cost_2023": [5, 9],    "cost_2024": [6, 10],
})
# wide_to_long splits a stub name from a suffix in one call.
tidy = pd.wide_to_long(messy, stubnames=["sales", "cost"], i="id",
                       j="year", sep="_").reset_index()
print("\nwide_to_long:\n", tidy.to_string(index=False))

# The general version: melt, then split the name column.
melted = messy.melt(id_vars="id", var_name="metric_year", value_name="value")
melted[["metric", "year"]] = melted["metric_year"].str.rsplit("_", n=1, expand=True)
print("\nmelt + split:\n",
      melted.drop(columns="metric_year").to_string(index=False))

# --- SCHEMA STABILITY: pin expected columns after any pivot. ---
EXPECTED_MONTHS = ["jan", "feb", "mar", "apr"]
stable = grid.reindex(columns=EXPECTED_MONTHS)
print("\npinned columns (apr present and null, not absent):\n",
      stable.to_string())
# Now downstream code that reads `apr` works whether or not April had data.

  THE SAME FOUR FACTS, TWO LAYOUTS

  WIDE                             LONG
  store  jan  feb                  store  month  sales
  A      100  110                  A      jan    100
  B      150  140                  A      feb    110
                                   B      jan    150
                                   B      feb    140

  adding March:  ALTER the schema   adding March:  INSERT rows
  human reading: easy               human reading: harder
  grouping:      enumerate columns  grouping:      one expression
  storage:       sparse -> nulls    storage:       dense, no nulls

  THE ROUND TRIP

       melt / unpivot / stack
   WIDE -------------------------> LONG
       <-------------------------
       pivot / unstack

    pivot       raises on duplicate (index, column) pairs   <- a safety net
    pivot_table aggregates them silently                    <- convenient

  RULE OF THUMB
    compute and store in long; pivot to wide as the LAST step, for humans.

Key Takeaways

• Wide is for human reading and some statistical routines; long is for computing, storing, and plotting.
• melt goes wide to long, pivot goes long to wide and refuses duplicates, and pivot_table aggregates them instead.
• Pivoting turns data values into column names, so pin the expected column set or the schema changes when the data does.
• Nulls created by pivoting mean "no observation", which is not the same as zero; filling them is a substantive claim.
• Keep pipelines in long format and pivot only at the final presentation step.

🧪 Practice

1. Convert a wide monthly-sales table to long, compute a per-month total, and pivot it back, verifying the round trip is lossless.
2. Take a wide table whose column names encode two variables, such as sales_2024, and produce a tidy long version with separate metric and year columns.
3. Pivot a sparse long table, then justify in one sentence each whether the resulting nulls should stay null or become zero.
4. Interview: A downstream job broke because a pivoted table gained a column. How would you have prevented it? (Hint: where did the column name come from, and what could you have declared instead?)

Window and Rolling Calculations

Some questions cannot be answered one row at a time. "What is this customer's spend relative to their previous order?" and "What is the seven-day average?" both need a row to see its neighbours. A window function computes a value for each row using a set of related rows, without collapsing the table the way an aggregation would.

The difference from a group-by aggregation is the key idea: an aggregation returns one row per group, a window function returns one value per original row. You keep the detail and gain the context.

Windows come in three useful varieties.

Offset windows look at a fixed number of rows away. shift(1) gives the previous value, which is how you compute differences, growth rates, and "time since last event". This is the simplest and most used.

Rolling windows cover a fixed span ending at the current row: a moving average over the last seven rows, or over the last seven days. The distinction between "rows" and "days" matters enormously for irregular data. Seven rows of a sensor that sometimes drops readings is not seven days. When the index is a timestamp, use a time-based window and the gaps take care of themselves.

Expanding and cumulative windows cover everything from the start up to the current row: a running total, a cumulative maximum, an all-time average to date. These are the honest way to build features from history, because by construction they can only see the past.

The three parameters that decide whether a window is correct:

1. Ordering. A window without an explicit sort is meaningless, and pandas will happily compute one on whatever order the rows happen to be in.
2. Partitioning. A rolling average over a table containing many customers must restart at each customer, or the first rows of one customer average in the last rows of another. In pandas this is a groupby before the rolling; in SQL it is PARTITION BY.
3. Alignment. By default a rolling window ends at the current row, which is what you want for features. A centred window uses rows on both sides, which is fine for smoothing a chart and is leakage if the output becomes a model feature, because it contains the future.

That last point is the one to burn in. rolling(7, center=True) and shift(-1) both reach forward. In a visualisation that is a legitimate choice. In a training feature it produces a model that scores brilliantly in validation and fails in production, because at prediction time the future rows do not exist.

import pandas as pd
import numpy as np

rng = np.random.default_rng(2)
days = pd.date_range("2024-01-01", periods=14, freq="D")
df = pd.DataFrame({
    "date": list(days) * 2,
    "store": ["A"] * 14 + ["B"] * 14,
    "sales": np.r_[100 + rng.normal(0, 8, 14).round(1),
                   200 + rng.normal(0, 15, 14).round(1)],
}).sort_values(["store", "date"]).reset_index(drop=True)

# --- OFFSET WINDOWS: the previous value, per partition. ---
g = df.groupby("store")["sales"]
df["prev"] = g.shift(1)                       # previous row WITHIN the store
df["delta"] = (df["sales"] - df["prev"]).round(1)
df["pct_change"] = (g.pct_change() * 100).round(1)
print(df.head(4).to_string(index=False))
print("\nfirst row of each store is null by construction:")
print(df.groupby("store").head(1)[["store", "date", "prev"]].to_string(index=False))

# WITHOUT the groupby, store B's first row borrows store A's last row.
df["prev_wrong"] = df["sales"].shift(1)
boundary = df.index[df["store"].ne(df["store"].shift(1))][1]
print(f"\nat the store boundary (row {boundary}):")
print(df.loc[[boundary], ["store", "sales", "prev", "prev_wrong"]]
        .to_string(index=False))

# --- ROLLING WINDOWS: rows vs time, and why it matters. ---
df["ma3"] = g.rolling(3).mean().round(1).reset_index(level=0, drop=True)
df["ma3_min1"] = (g.rolling(3, min_periods=1).mean().round(1)
                  .reset_index(level=0, drop=True))
print("\n", df[df["store"] == "A"].head(4)[
    ["date", "sales", "ma3", "ma3_min1"]].to_string(index=False))
# min_periods=1 trades a full window for no leading nulls. Which you want
# depends on whether a partial average is meaningful for your metric.

# TIME-BASED windows handle irregular spacing correctly.
irregular = pd.DataFrame({
    "sales": [10, 20, 30, 40],
}, index=pd.to_datetime(["2024-01-01", "2024-01-02", "2024-01-09", "2024-01-10"]))
print("\nrows-based 3:\n", irregular["sales"].rolling(3).sum().to_string())
print("\ntime-based 3D:\n", irregular["sales"].rolling("3D").sum().to_string())
# The rows window spans 9 days at one point. The time window does not,
# because it asks "the last 3 days", which is what the metric meant.

# --- EXPANDING / CUMULATIVE: only ever sees the past. ---
df["cum_sales"] = g.cumsum().round(1)
df["running_avg"] = g.expanding().mean().round(1).reset_index(level=0, drop=True)
df["running_max"] = g.cummax()
print("\n", df[df["store"] == "A"].head(4)[
    ["date", "sales", "cum_sales", "running_avg", "running_max"]]
    .to_string(index=False))

# --- LEAKAGE: centred and forward-looking windows reach into the future. ---
a = df[df["store"] == "A"].set_index("date")["sales"]
comparison = pd.DataFrame({
    "sales": a,
    "trailing_3": a.rolling(3).mean().round(1),            # SAFE as a feature
    "centred_3": a.rolling(3, center=True).mean().round(1),  # LEAKS
    "next_day": a.shift(-1),                                # LEAKS outright
})
print("\n", comparison.head(5).to_string())
print("\nrow 2's centred mean uses row 3, which has not happened yet.")

# EWM: an exponentially weighted mean, trailing and smooth. A good default
# for features when a hard window boundary feels arbitrary.
print("\nEWM (halflife=3):\n", a.ewm(halflife=3).mean().round(1).head().to_string())

# --- RANK AND PERCENTILE WITHIN A WINDOW ---
df["rank_in_store"] = g.rank(method="dense", ascending=False)
df["pct_in_store"] = g.rank(pct=True).round(2)

# --- LAG FEATURES FOR FORECASTING: build them explicitly and check them. ---
def add_lags(frame, col, group, lags=(1, 7), windows=(7, 14)):
    out = frame.sort_values([group, "date"]).copy()
    gb = out.groupby(group)[col]
    for k in lags:
        out[f"{col}_lag{k}"] = gb.shift(k)
    for w in windows:
        # shift(1) BEFORE rolling so the window excludes the current row.
        out[f"{col}_ma{w}"] = (gb.shift(1).rolling(w, min_periods=2).mean()
                               .reset_index(level=0, drop=True).round(1))
    return out

feat = add_lags(df[["date", "store", "sales"]], "sales", "store")
print("\nlag features (note the leading nulls are correct, not a bug):")
print(feat[feat["store"] == "A"].head(9).to_string(index=False))

# The audit: every feature column must be null wherever it would have
# required a row the model could not have seen.
assert feat.groupby("store")["sales_lag7"].apply(
    lambda s: s.head(7).isna().all()).all()
print("\nlag audit passed: no feature uses a row from the future")

  WINDOW ALIGNMENT AND LEAKAGE

  rows:        t-3   t-2   t-1    t    t+1   t+2
  trailing(3)   [     x     x     x ]              SAFE: ends at t
  centred(3)               [ x     x     x ]       LEAKS: includes t+1
  shift(-1)                        ->    x         LEAKS: IS t+1
  expanding     [ x  x  x  x  x ]                  SAFE: starts at the beginning

  Use centred windows for charts. Never for features.

  THE THREE THINGS A WINDOW NEEDS

    1. ORDER        sort_values("date")     or the window is arbitrary
    2. PARTITION    groupby("store")        or windows bleed across entities
    3. ALIGNMENT    trailing, not centred   or the feature sees the future

  ROWS VS TIME

    index:  Jan 1   Jan 2   Jan 9   Jan 10
    rolling(3)     spans 9 days at Jan 9   <- probably not the intent
    rolling("3D")  spans 3 days always     <- what "3-day average" means

Key Takeaways

• Window functions add context to every row instead of collapsing rows into groups, which is what makes them the backbone of feature building.
• A window is only meaningful with an explicit sort, and it must be partitioned by entity or values bleed across entities.
• Use time-based rolling windows for irregularly spaced data, because a count of rows is not a span of time.
• Centred windows and negative shifts include the future; they are fine for smoothing a chart and are leakage in a model feature.
• Shift by one before rolling when a feature must exclude the current row, and assert that early rows are null.

🧪 Practice

1. Compute a seven-day trailing average per entity and verify that the first six rows of each entity are null.
2. Build an irregular daily series with gaps and show that a three-row window and a three-day window give different answers.
3. Produce trailing, centred, and lead versions of a rolling mean, and mark which of them could be used as a model feature.
4. Interview: A sales forecasting model achieves a mean absolute error of 2% in backtesting and 25% live. Which feature would you inspect first? (Hint: which of your windows could have touched a row after the prediction timestamp?)

Binning and Discretization

Binning converts a continuous variable into a small number of ordered categories: age into age brackets, income into deciles, a score into low, medium, and high. It seems like throwing information away, and it is, so the question is when the trade is worth it.

There are four genuine reasons to bin.

Non-linearity without a non-linear model. If risk rises sharply after age sixty and is flat before, a linear term fits neither region well. Binning lets a linear model capture the step for free.

Robustness to outliers and noise. A bin label is unmoved by a value being 100 or 10,000, so heavy tails stop dominating.

Communication. Stakeholders reason about "the top decile" more readily than about a coefficient per unit. Regulated scorecards are built from bins for exactly this reason.

Privacy and sparsity. Reporting an age bracket instead of a birthdate reduces re-identification risk, and coarse buckets have enough rows per cell to be estimated stably.

Against that, binning discards within-bin variation, creates arbitrary cliffs where a value one unit either side of a boundary lands in different categories, and usually costs predictive accuracy compared with letting a flexible model use the raw variable. Modern gradient boosting finds its own splits, so for pure prediction it rarely needs help.

The methods differ in where they put the boundaries.

Equal-width divides the range into intervals of the same size. Simple, interpretable, and terrible on skewed data, where almost every observation lands in the first bin.

Equal-frequency, or quantile binning, puts the same number of observations in each bin. Robust to skew and the usual default. Its boundary values are data-dependent, which means they must be treated as fitted parameters. It also breaks down when a value is very common, because a single repeated value cannot be split across bins.

Domain boundaries come from knowledge: legal age thresholds, clinical reference ranges, pricing tiers. When they exist they beat both of the above, because they mean something to the reader.

Supervised binning chooses boundaries to maximise separation of the target, for example by fitting a shallow decision tree on the one variable and reading off its splits. Powerful and the most prone to overfitting, so it must be fitted on training data only and validated.

Two implementation details that cause real bugs. First, bin edges are fitted parameters: compute them on training data and reuse them, or your test set gets different bins and your deployed model gets different bins again. Second, handle values outside the training range. A test value above the largest training edge falls outside every bin and becomes null unless you extend the outermost edges to infinity.

import numpy as np
import pandas as pd

rng = np.random.default_rng(5)
income = pd.Series(rng.lognormal(10.5, 0.8, 1000).round(0), name="income")
print(f"income: min {income.min():,.0f}  median {income.median():,.0f}"
      f"  max {income.max():,.0f}  skew {income.skew():.2f}")

# --- EQUAL WIDTH: intuitive edges, useless on skewed data. ---
eq_width = pd.cut(income, bins=5)
print("\nequal width:\n", eq_width.value_counts().sort_index().to_string())
# Almost everything lands in bin 1. The bins are "correct" and uninformative.

# --- EQUAL FREQUENCY: equal counts, data-dependent edges. ---
eq_freq, edges = pd.qcut(income, q=5, retbins=True,
                         labels=["Q1", "Q2", "Q3", "Q4", "Q5"])
print("\nequal frequency:\n", eq_freq.value_counts().sort_index().to_string())
print("edges:", np.round(edges, 0).tolist())

# THE DUPLICATE-EDGE PROBLEM: a very common value cannot be split.
spiky = pd.Series([0] * 600 + list(rng.integers(1, 100, 400)))
try:
    pd.qcut(spiky, q=5)
except ValueError as e:
    print("\nqcut on a spiky column:", str(e)[:56])
print("with duplicates='drop':",
      pd.qcut(spiky, q=5, duplicates="drop").cat.categories.size, "bins produced")
# Better still: separate the spike explicitly before binning the remainder.

# --- DOMAIN BOUNDARIES: meaningful, stable, explainable. ---
age = pd.Series([17, 18, 24, 25, 39, 40, 64, 65, 80])
AGE_EDGES = [0, 18, 25, 40, 65, np.inf]
AGE_LABELS = ["minor", "18-24", "25-39", "40-64", "65+"]
banded = pd.cut(age, bins=AGE_EDGES, labels=AGE_LABELS, right=False)
print("\n", pd.DataFrame({"age": age, "band": banded}).to_string(index=False))
# right=False makes intervals [lo, hi), so 18 is "18-24" not "minor".
# Off-by-one at a legal threshold is a real and embarrassing class of bug.

# --- FITTED EDGES: compute on train, apply everywhere. ---
from sklearn.model_selection import train_test_split
tr, te = train_test_split(income, test_size=0.3, random_state=0)

_, fitted = pd.qcut(tr, q=4, retbins=True)
# EXTEND the outer edges, or test values beyond the training range go null.
safe_edges = fitted.copy()
safe_edges[0], safe_edges[-1] = -np.inf, np.inf

naive_te = pd.cut(te, bins=fitted)
safe_te = pd.cut(te, bins=safe_edges)
print(f"\ntest rows binned to null: naive {naive_te.isna().sum()},"
      f" extended edges {safe_te.isna().sum()}")

# The scikit-learn equivalent, which stores the edges for you.
from sklearn.preprocessing import KBinsDiscretizer
disc = KBinsDiscretizer(n_bins=4, encode="ordinal", strategy="quantile",
                        subsample=None)
disc.fit(tr.to_frame())
print("learned edges:", np.round(disc.bin_edges_[0], 0).tolist())
print("transformed test head:",
      disc.transform(te.to_frame())[:5].ravel().tolist())

# --- SUPERVISED BINNING: let a shallow tree choose the cut points. ---
from sklearn.tree import DecisionTreeClassifier
x = rng.normal(45, 15, 2000).clip(18, 90)
# A genuine step in risk at 60, which no unsupervised binning would find.
y = rng.random(2000) < np.where(x > 60, 0.45, 0.10)

tree = DecisionTreeClassifier(max_depth=2, min_samples_leaf=100, random_state=0)
tree.fit(x.reshape(-1, 1), y)
cuts = sorted(tree.tree_.threshold[tree.tree_.feature == 0].round(1))
print(f"\ntree-selected cut points: {cuts}")

sup_edges = [-np.inf, *cuts, np.inf]
sup_bins = pd.cut(x, bins=sup_edges)
rates = pd.DataFrame({"bin": sup_bins, "y": y}).groupby(
    "bin", observed=True)["y"].agg(["mean", "size"])
print("\nevent rate per supervised bin:\n", rates.round(3).to_string())
# Compare with quantile bins, which spread the step across two bins:
q_rates = pd.DataFrame({"bin": pd.qcut(x, 4), "y": y}).groupby(
    "bin", observed=True)["y"].mean()
print("\nevent rate per quantile bin:\n", q_rates.round(3).to_string())

# --- WHAT BINNING COSTS: the cliff at a boundary. ---
pair = pd.Series([59.9, 60.1])
print("\ntwo nearly identical values, different bins:")
print(pd.DataFrame({"x": pair, "bin": pd.cut(pair, bins=sup_edges)})
      .to_string(index=False))

  WHERE THE EDGES GO

  skewed data, 5 bins

  EQUAL WIDTH            EQUAL FREQUENCY         DOMAIN
  |#########|  920       |##|  200               |###|  legal thresholds
  |#|         60         |##|  200               |###|  clinical ranges
  ||          15         |##|  200               |###|  pricing tiers
  ||           4         |##|  200
  ||           1         |##|  200
  interpretable edges,   equal power per bin,    meaningful to the reader,
  useless counts         edges are parameters    needs domain knowledge

  BIN OR NOT?

    goal                          bin?     why
    ---------------------------------------------------------------
    linear model, known step      yes      captures non-linearity
    gradient boosting             no       it finds its own splits
    stakeholder-facing report     yes      deciles communicate
    regulated scorecard           yes      required for explainability
    maximum predictive accuracy   no       discarding variation costs
    privacy / small cells         yes      coarser cells, less risk

  ALWAYS: edges fitted on train, outer edges at +/- infinity.

Key Takeaways

• Binning trades information for non-linearity, robustness, interpretability, and privacy, so bin only when one of those is the goal.
• Equal-width bins are meaningless on skewed data; quantile bins are the robust default, and domain boundaries beat both when they exist.
• Supervised binning finds the real step in the target and overfits the most, so fit it on training data and validate.
• Bin edges are fitted parameters, and the outermost edges must run to plus and minus infinity or unseen values become null.
• Every boundary creates a cliff where nearly identical values are treated as different; check that no such cliff sits where a decision is made.

🧪 Practice

1. Bin a log-normal column with equal-width and equal-frequency methods and compare the resulting counts per bin.
2. Fit quantile edges on a training split, apply them to a test split without extending the outer edges, and count the nulls produced.
3. Fit a depth-two decision tree on one feature to derive supervised bins, then compare event rates per bin against quantile bins.
4. Interview: A credit model bins income into deciles and a customer complains that a fifty-dollar raise reduced their approved limit. What happened, and how would you change the feature? (Hint: what happens at a boundary, and what continuous alternative keeps monotonicity?)

Tidy Data Principles

"Tidy data" is a specific definition, not a synonym for clean, and having the definition makes a large class of messy-data problems mechanical instead of improvised. A dataset is tidy when:

1. Each variable forms a column.
2. Each observation forms a row.
3. Each type of observational unit forms a table.

That is it. The payoff is that once data is in this shape, the same small set of tools works on it. You stop writing bespoke code per dataset, because grouping, joining, plotting, and modelling all assume this layout. The idea comes from Hadley Wickham's 2014 paper and underlies both the tidyverse in R and most of the pandas idioms you already use.

The value is in the diagnosis. Messy data is messy in a handful of recognisable ways, and each has a standard fix.

Column headers are values, not variable names. Columns named 2022, 2023, 2024 mean year is a variable stored as headers. Fix: melt.

Multiple variables in one column. A column named m_18-25 encodes both sex and age band. Fix: melt, then split the name.

Variables in both rows and columns. A element column containing tmax and tmin for each date means the measurement type is a variable in rows while dates sit in columns. Fix: melt, then pivot the measurement type back out.

Multiple observational units in one table. A table with one row per song play, repeating the song's title and artist on every row, mixes two units. Fix: normalise into a plays table and a songs table with a key. This is the tidy-data version of database normalisation, and it matters because repeated attributes get updated inconsistently.

One observational unit split across tables. Twelve monthly files with identical schemas. Fix: concatenate, adding a column for whatever the file name encoded, because that is a variable too.

Two caveats worth stating, because tidy is a principle and not a law. Tidy data is often less efficient to store than wide data, since identifier columns repeat on every row. And some legitimate operations want wide input: correlation matrices, many matrix-oriented algorithms, and anything a human reads. The principle applies to your working format, and pivoting for a specific purpose is not a violation.

The one place tidiness is worth defending strictly is the boundary of a pipeline, because a table whose grain is undefined is a table nobody can join to correctly.

import pandas as pd
import numpy as np

# --- MESS 1: column headers are values. ---
messy1 = pd.DataFrame({
    "religion": ["Agnostic", "Atheist", "Buddhist"],
    "<$10k": [27, 12, 27], "$10-20k": [34, 27, 21], "$20-30k": [60, 37, 30],
})
print("headers are income bands, i.e. a variable:\n", messy1.to_string(index=False))
tidy1 = messy1.melt(id_vars="religion", var_name="income", value_name="count")
print("\ntidy:\n", tidy1.to_string(index=False))

# --- MESS 2: two variables inside one column name. ---
messy2 = pd.DataFrame({
    "country": ["AD", "AE"],
    "m_014": [0, 2], "m_1524": [10, 16], "f_014": [1, 3], "f_1524": [8, 12],
})
tidy2 = messy2.melt(id_vars="country", var_name="sex_age", value_name="cases")
tidy2[["sex", "age"]] = tidy2["sex_age"].str.split("_", expand=True)
tidy2 = tidy2.drop(columns="sex_age")
print("\ntwo variables split out:\n", tidy2.to_string(index=False))

# --- MESS 3: variables in both rows and columns. ---
messy3 = pd.DataFrame({
    "station": ["MX17", "MX17", "MX17", "MX17"],
    "month":   ["2010-01", "2010-01", "2010-02", "2010-02"],
    "element": ["tmax", "tmin", "tmax", "tmin"],   # a variable stored in rows
    "d1": [27.8, 14.5, 29.6, 15.8],
    "d2": [27.3, 14.4, 29.2, 15.9],
})
print("\nelement is a variable in rows, days are variables in columns:")
print(messy3.to_string(index=False))

long3 = messy3.melt(id_vars=["station", "month", "element"],
                    var_name="day", value_name="value")
long3["date"] = pd.to_datetime(
    long3["month"] + "-" + long3["day"].str.lstrip("d").str.zfill(2))
tidy3 = (long3.pivot(index=["station", "date"], columns="element",
                     values="value")
              .reset_index().rename_axis(columns=None))
print("\ntidy: one row per station-date, one column per measurement:")
print(tidy3.to_string(index=False))
# Note this needed BOTH operations: melt the days, then pivot the element.

# --- MESS 4: two observational units in one table. ---
plays = pd.DataFrame({
    "played_at": pd.to_datetime(["2024-01-01 10:00", "2024-01-01 11:00",
                                 "2024-01-02 09:00"]),
    "user": ["u1", "u2", "u1"],
    "track": ["Hey Jude", "Hey Jude", "Yesterday"],
    "artist": ["The Beatles", "The Beatles", "The Beatles"],
    "duration_s": [431, 431, 125],
})
print("\ntrack attributes repeat on every play:\n", plays.to_string(index=False))

tracks = (plays[["track", "artist", "duration_s"]].drop_duplicates()
          .reset_index(drop=True).rename_axis("track_id").reset_index())
fact = plays.merge(tracks[["track_id", "track"]], on="track").drop(
    columns=["track", "artist", "duration_s"])
print("\ntracks (one row per track):\n", tracks.to_string(index=False))
print("\nplays (one row per play):\n", fact.to_string(index=False))
# Now a duration correction is one edit, not one edit per play row.

# --- MESS 5: one unit split across files. ---
jan = pd.DataFrame({"store": ["A", "B"], "sales": [100, 150]})
feb = pd.DataFrame({"store": ["A", "B"], "sales": [110, 140]})
# The file name encoded the month, so it becomes a column.
combined = pd.concat([jan.assign(month="2024-01"), feb.assign(month="2024-02")],
                     ignore_index=True)
print("\nconcatenated with the filename variable recovered:")
print(combined.to_string(index=False))

# --- THE PAYOFF: one vocabulary works on any tidy table. ---
print("\nthe same three operations on a tidy frame:")
print(tidy2.groupby(["sex"], as_index=False)["cases"].sum().to_string(index=False))
print(tidy2[tidy2["age"] == "1524"].to_string(index=False))
print(tidy2.pivot_table(index="country", columns="sex", values="cases",
                        aggfunc="sum").to_string())

# --- DECLARE THE GRAIN. A table without a stated grain cannot be joined. ---
def declare_grain(df: pd.DataFrame, keys: list[str], name: str) -> pd.DataFrame:
    if df.duplicated(keys).any():
        raise ValueError(f"{name} is not one row per {keys}")
    print(f"{name}: grain verified as one row per {keys}")
    return df

declare_grain(tidy3, ["station", "date"], "weather")
declare_grain(tracks, ["track_id"], "tracks")

  THE THREE RULES

    variable -> column     observation -> row     unit -> table

  MESSY PATTERN                            FIX
  ---------------------------------------------------------------------
  headers are values (2022, 2023)          melt
  one column holds two variables           melt, then split the name
  variables in rows AND columns            melt, then pivot
  two units in one table                   split into two tables + key
  one unit across many files               concat, keep the file variable
  values stored in column names            melt

  UNIT vs TABLE

    NOT TIDY                         TIDY
    plays                            plays              tracks
    +-------+--------+----------+    +-------+------+    +------+--------+
    | time  | track  | duration |    | time  | t_id |    | t_id | durat. |
    | 10:00 | HeyJude| 431      |    | 10:00 | 1    |    | 1    | 431    |
    | 11:00 | HeyJude| 431      |    | 11:00 | 1    |    +------+--------+
    +-------+--------+----------+    +-------+------+
    duration repeated -> can          duration stored once -> cannot
    disagree with itself              disagree with itself

Key Takeaways

• Tidy means one column per variable, one row per observation, one table per observational unit, and it is a definition rather than a vibe.
• Recognising which of the standard messy patterns you have turns cleanup into a mechanical melt, split, pivot, or normalise.
• Repeated attributes across rows are a sign that two observational units share a table, and they will eventually disagree with themselves.
• Information encoded in file names is a variable and belongs in a column.
• Tidy applies to the working format; pivoting to wide for reading or for a matrix algorithm is not a violation.

🧪 Practice

1. Take a table with years as column headers and tidy it, then reproduce the original with a single pivot.
2. Tidy a table that has a measurement-type column in rows and day numbers in columns, which requires a melt followed by a pivot.
3. Split a denormalised event log into a fact table and a dimension table, then show that correcting one attribute now requires a single edit.
4. Interview: A teammate argues that tidy data wastes storage because identifier columns repeat. How do you respond? (Hint: what does columnar compression do with a repeated value, and what is the cost of the alternative when the schema must change?)

<a id="63-feature-engineering"></a>

6.3 Feature Engineering

Feature engineering is the work of converting what you know about a problem into numbers a model can use; this section covers the standard transformations, the encodings that categorical data requires, and the leakage traps that make a good validation score meaningless.

Scaling and Normalization

Many algorithms compare features to each other, and features measured in different units are not comparable. Age runs from 18 to 90. Income runs from 20,000 to 400,000. A distance calculation between two people is dominated almost entirely by income, not because income matters more but because its numbers are larger. Scaling removes that accident of units.

The set of algorithms that care is specific, and knowing which ones saves a lot of pointless preprocessing.

Distance-based methods care: k-nearest neighbours, k-means, support vector machines, and anything using a radial basis kernel. Distance is literally a sum over features, so the largest-magnitude feature wins.

Gradient-descent methods care: linear and logistic regression trained iteratively, and all neural networks. With features on wildly different scales the loss surface becomes a long narrow valley, and gradient descent zigzags down it slowly. Scaling makes it round and convergence fast.

Regularised models care: ridge and lasso penalise the size of coefficients, and a coefficient's size depends on its feature's units. Without scaling the penalty falls unequally on features for no principled reason.

PCA cares, because it finds directions of maximum variance and variance is unit-dependent. Unscaled PCA on data containing income will return income as the first component every time.

Tree-based models do not care. A decision tree asks whether a feature exceeds a threshold, and that question is unchanged by any monotonic rescaling. Scaling before a random forest or gradient boosting is harmless and useless.

The main methods:

Method	Formula	Output range	Outlier-robust
Standardisation	(x - mean) / sd	unbounded, mean 0	no
Min-max	(x - min) / (max - min)	[0, 1]	no, very sensitive
Robust	(x - median) / IQR	unbounded	yes
Max-abs	x / max abs	[-1, 1], keeps zeros	no
Quantile	rank to uniform or normal	bounded	yes, strongly
L2 normalisation	x / norm(x) per row	unit length row	n/a, row-wise

Standardisation is the default. It centres at zero and gives unit variance. It does not bound the output, which is fine, and it does not make the distribution normal, which is a common misunderstanding: a skewed variable remains skewed after standardisation.

Min-max is what you use when you need a bounded range, for example for image pixel values or a neural network with a saturating activation. Its weakness is severe: a single extreme value sets the maximum and squashes everything else into a narrow band near zero.

Robust scaling uses the median and interquartile range, so outliers do not move the parameters. This is the right default when you know the data has heavy tails and you do not want to remove them.

Quantile transformation maps values to their ranks and then to a uniform or normal distribution. It is the most aggressive option: it will normalise almost anything, at the cost of destroying the original spacing between values entirely.

Two rules matter more than the choice of method. First, fit on training data only, for the same reason as every other fitted transform: the mean and standard deviation are statistics, and computing them over the test set leaks. Second, note that unseen values can land outside the training range. Min-max scaling on a test value above the training maximum gives a number above one. That is usually fine, but if your pipeline asserts a range, it will fail.

import numpy as np
import pandas as pd
from sklearn.preprocessing import (StandardScaler, MinMaxScaler, RobustScaler,
                                   MaxAbsScaler, QuantileTransformer)
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier

rng = np.random.default_rng(4)
n = 800
X = pd.DataFrame({
    "age": rng.integers(18, 90, n).astype(float),
    "income": rng.lognormal(10.8, 0.6, n).round(0),
})
X.loc[0, "income"] = 5_000_000            # one extreme value, on purpose
y = ((X["age"] > 55) ^ (X["income"] > X["income"].median())).astype(int)

# --- WHY IT MATTERS: distance is dominated by the larger-scale column. ---
p, q = X.iloc[1], X.iloc[2]
d_age = (p["age"] - q["age"]) ** 2
d_inc = (p["income"] - q["income"]) ** 2
print(f"squared distance contribution: age {d_age:,.0f}  income {d_inc:,.0f}")
print(f"income accounts for {100 * d_inc / (d_age + d_inc):.4f}% of the distance")

# --- THE METHODS, SIDE BY SIDE ---
scalers = {
    "standard": StandardScaler(),
    "minmax": MinMaxScaler(),
    "robust": RobustScaler(),
    "maxabs": MaxAbsScaler(),
    "quantile-normal": QuantileTransformer(output_distribution="normal",
                                           n_quantiles=200, random_state=0),
}
rows = {}
for name, sc in scalers.items():
    z = pd.Series(sc.fit_transform(X[["income"]]).ravel())
    rows[name] = {"min": z.min(), "p50": z.median(), "p99": z.quantile(0.99),
                  "max": z.max(), "skew": z.skew()}
print("\nincome after each scaler:")
print(pd.DataFrame(rows).T.round(2).to_string())
# Read the minmax row: the outlier pushed the 99th percentile to ~0.03, so
# 99% of the data now occupies 3% of the range. Standardisation keeps the
# spread but not the skew. Only the quantile transform removes the skew.

std_income = pd.Series(StandardScaler().fit_transform(X[["income"]]).ravel())
print(f"\nstandardisation does NOT normalise: "
      f"raw skew {X['income'].skew():.2f}, "
      f"standardised skew {std_income.skew():.2f}")

# --- WHO CARES, MEASURED ---
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0)
for model_name, model in [("kNN", KNeighborsClassifier(15)),
                          ("random forest", RandomForestClassifier(
                              n_estimators=80, random_state=0))]:
    bare = cross_val_score(model, Xtr, ytr, cv=5).mean()
    scaled = cross_val_score(
        Pipeline([("s", StandardScaler()), ("m", model)]), Xtr, ytr, cv=5).mean()
    print(f"{model_name:14s} raw {bare:.3f}   scaled {scaled:.3f}"
          f"   delta {scaled - bare:+.3f}")
# kNN improves substantially. The forest does not move, because a threshold
# test is invariant to any monotonic rescaling.

# --- LEAKAGE: the scaler's parameters are statistics. ---
leaky_mean = X["income"].mean()                     # uses test rows
clean_mean = Xtr["income"].mean()                   # train only
print(f"\nmean from all data {leaky_mean:,.0f} vs train only {clean_mean:,.0f}")

# A Pipeline makes the correct behaviour automatic inside cross-validation:
# the scaler is re-fitted on each training fold, never on the held-out fold.
pipe = Pipeline([("scale", StandardScaler()), ("model", KNeighborsClassifier(15))])
print("pipeline CV (no leakage):",
      round(cross_val_score(pipe, Xtr, ytr, cv=5).mean(), 3))

# --- UNSEEN VALUES ESCAPE THE TRAINING RANGE ---
mm = MinMaxScaler().fit(Xtr[["income"]])
te_scaled = mm.transform(Xte[["income"]])
print(f"\ntest range after min-max fitted on train: "
      f"[{te_scaled.min():.3f}, {te_scaled.max():.3f}]")
print(f"values outside [0,1]: {((te_scaled < 0) | (te_scaled > 1)).sum()}")
# clip=True bounds the output, at the cost of hiding that it happened.

# --- ROW-WISE NORMALISATION is a different operation entirely. ---
from sklearn.preprocessing import Normalizer
counts = pd.DataFrame({"cat": [3, 30], "dog": [1, 10], "bird": [0, 0]})
print("\nraw counts:\n", counts.to_string(index=False))
print("\nL2-normalised rows (direction, not magnitude):")
print(pd.DataFrame(Normalizer().fit_transform(counts),
                   columns=counts.columns).round(3).to_string(index=False))
# Both rows have the same composition, so after row normalisation they are
# identical. Use this for term counts, where proportions carry the meaning.

# --- SPARSE DATA: centring destroys sparsity. ---
print("\nStandardScaler(with_mean=False) or MaxAbsScaler preserve zeros,")
print("which matters when a matrix is 99% zeros and would not fit dense.")

  WHO NEEDS SCALING

  algorithm                    needs it?   why
  ----------------------------------------------------------------------
  kNN, k-means, SVM/RBF        YES         distance sums over features
  linear/logistic + gradient   YES         conditioning of the loss surface
  ridge, lasso, elastic net    YES         penalty depends on units
  PCA, factor analysis         YES         variance depends on units
  neural networks              YES         activation saturation, convergence
  decision tree                no          threshold test is scale-invariant
  random forest, boosting      no          same reason
  Naive Bayes (Gaussian)       no          fits per-feature parameters

  WHAT AN OUTLIER DOES TO MIN-MAX

    raw:      20k ....... 400k                          5,000k
    minmax:   |####################|                         |
              0                 0.08                         1
              99% of the data is squeezed into 8% of the range.

    robust:   median-centred, IQR-scaled -> the outlier is far out,
              and everything else keeps its spread.

  ALWAYS: fit on train inside a Pipeline. Never fit on the full dataset.

Key Takeaways

• Scaling matters for distance-based, gradient-based, regularised, and variance-based methods, and not at all for tree ensembles.
• Standardisation is the default; it centres and rescales but does not make a skewed distribution normal.
• Min-max is for bounded output and is destroyed by a single extreme value; robust scaling uses the median and IQR instead.
• Quantile transformation will normalise almost anything at the cost of the original spacing between values.
• Fit scalers inside a pipeline so each cross-validation fold re-fits on its own training rows, and expect test values outside the training range.

🧪 Practice

1. Add one extreme value to a column, apply min-max and robust scaling, and report where the 99th percentile lands under each.
2. Compare cross-validated accuracy of k-nearest neighbours with and without standardisation, then repeat for a random forest and explain the difference.
3. Fit a scaler on the full dataset and then inside a pipeline, and show the gap between the two cross-validated scores.
4. Interview: Your model uses both age and annual revenue and is a k-nearest-neighbours classifier. A colleague says scaling is a detail. What number would you show them? (Hint: decompose a single pairwise distance into per-feature contributions.)

Encoding Categorical Variables

Models consume numbers, and a category is not a number. Encoding is how you bridge that gap, and the choice of encoding is one of the highest-leverage decisions in feature engineering, because a bad one either invents structure that does not exist or explodes your feature count.

The first distinction is whether the categories have a meaningful order.

Nominal categories do not: red, green, blue. There is no sense in which green is between red and blue. Ordinal categories do: small, medium, large. The order is real and a model can use it.

The classic mistake is label encoding a nominal variable: mapping red to 0, green to 1, blue to 2. A linear model now believes green is halfway between red and blue, and that blue is twice green. That relationship is fabricated. For ordinal variables the same encoding is correct and useful, provided the integers follow the real order rather than alphabetical order, which is what the default gives you.

One-hot encoding is the standard answer for nominal variables: one binary column per category. It asserts no order and no distance, which is exactly right. Its cost is width, one column per category, and that cost becomes prohibitive for high-cardinality variables.

Three details on one-hot that matter in practice.

The dummy trap. If you include a column for every category plus an intercept, the columns are perfectly collinear, because they sum to one. For linear models with an intercept, drop one category as a reference. For regularised models and trees it does not matter, and for interpretation keeping all columns is sometimes clearer.

Unseen categories. A category appearing in production but not in training has no column. The encoder must be told what to do: handle_unknown="ignore" produces all zeros, which is the sensible default. Silently failing is not.

Sparsity. One-hot output is mostly zeros, and for wide categoricals a sparse matrix representation is the difference between fitting in memory and not.

Ordinal encoding for genuinely ordered categories is a single column, and you should pass the order explicitly rather than letting the library sort alphabetically. ["low", "medium", "high"] sorts to ["high", "low", "medium"], which is nonsense.

Binary and hashing encodings compress: binary represents the category index in base two, using about log2(k) columns; hashing maps categories to a fixed number of columns via a hash function, accepting collisions. Both trade interpretability for width and are mostly superseded by target encoding for tabular work, though hashing remains useful for streaming data where the category set is not known up front.

A final, practical note: tree-based models tolerate ordinal encoding of nominal variables better than linear models do, because a tree can carve out individual values with successive splits. It is still not ideal, and modern implementations such as LightGBM and CatBoost accept categorical columns natively, which is better than either.

import numpy as np
import pandas as pd
from sklearn.preprocessing import OneHotEncoder, OrdinalEncoder
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LinearRegression

df = pd.DataFrame({
    "colour": ["red", "green", "blue", "red", "green"],
    "size":   ["small", "large", "medium", "medium", "small"],
    "price":  [10.0, 30.0, 20.0, 18.0, 11.0],
})

# --- THE FABRICATED-ORDER BUG ---
naive = df["colour"].map({"red": 0, "green": 1, "blue": 2})
print("label-encoded nominal:", naive.tolist())
print("this asserts green - red == blue - green, which is meaningless")

# A linear model will use that fake ordering as if it were real.
lm = LinearRegression().fit(naive.to_frame(), df["price"])
print(f"fitted slope on the fake scale: {lm.coef_[0]:.2f} per 'colour unit'")

# --- ONE-HOT: no order, no distance. ---
ohe = OneHotEncoder(sparse_output=False, handle_unknown="ignore")
encoded = pd.DataFrame(ohe.fit_transform(df[["colour"]]),
                       columns=ohe.get_feature_names_out(["colour"]))
print("\none-hot:\n", encoded.astype(int).to_string(index=False))

# THE DUMMY TRAP: the columns sum to 1, so with an intercept they are
# perfectly collinear and the coefficients are not identified.
print("\nrow sums:", encoded.sum(axis=1).tolist(), "-> collinear with intercept")
dropped = OneHotEncoder(sparse_output=False, drop="first")
cols = dropped.fit_transform(df[["colour"]])
print(f"drop='first' gives {cols.shape[1]} columns instead of {encoded.shape[1]};"
      f" the omitted category is the reference level")

# UNSEEN CATEGORIES: decide the behaviour, do not discover it.
print("\nunseen 'purple' with handle_unknown='ignore':",
      ohe.transform(pd.DataFrame({"colour": ["purple"]})).astype(int).tolist())
try:
    OneHotEncoder(handle_unknown="error").fit(df[["colour"]]).transform(
        pd.DataFrame({"colour": ["purple"]}))
except ValueError as e:
    print("with 'error':", str(e)[:52])

# --- ORDINAL: pass the real order. Alphabetical is not the order. ---
print("\nalphabetical sort of sizes:", sorted(df["size"].unique()))
SIZE_ORDER = ["small", "medium", "large"]
oe = OrdinalEncoder(categories=[SIZE_ORDER],
                    handle_unknown="use_encoded_value", unknown_value=-1)
df["size_ord"] = oe.fit_transform(df[["size"]]).astype(int)
print(df[["size", "size_ord"]].to_string(index=False))
# Now the number carries the real ordering, and a monotone model can use it.

# Pandas' ordered Categorical does the same and supports comparisons.
df["size_cat"] = pd.Categorical(df["size"], categories=SIZE_ORDER, ordered=True)
print("\ncomparison works on an ordered Categorical:",
      (df["size_cat"] > "small").tolist())

# --- WIDTH: one-hot cost grows with cardinality. ---
card = pd.Series(np.arange(1, 6)).map(lambda k: 10 ** k)
print("\ncolumns produced by each encoding:")
print(pd.DataFrame({
    "categories": card,
    "one_hot": card,
    "binary": np.ceil(np.log2(card)).astype(int),
    "hashing_256": 256,
    "target": 1,
}).to_string(index=False))

# --- BINARY ENCODING: the category index written in base two. ---
def binary_encode(s: pd.Series, prefix: str) -> pd.DataFrame:
    codes = pd.Categorical(s).codes + 1              # 0 is reserved for unseen
    width = int(np.ceil(np.log2(codes.max() + 1)))
    bits = ((codes[:, None] >> np.arange(width)) & 1)
    return pd.DataFrame(bits, columns=[f"{prefix}_b{i}" for i in range(width)],
                        index=s.index)

many = pd.Series([f"cat_{i % 12}" for i in range(12)], name="c")
print(f"\nbinary encoding of 12 categories -> "
      f"{binary_encode(many, 'c').shape[1]} columns")
print(binary_encode(many, "c").head(4).to_string(index=False))

# --- HASHING: fixed width, no fitted vocabulary, accepts collisions. ---
from sklearn.feature_extraction import FeatureHasher
hasher = FeatureHasher(n_features=8, input_type="string")
hashed = hasher.transform([[v] for v in many]).toarray()
print(f"\nhashed to 8 columns; distinct rows: {len(np.unique(hashed, axis=0))}"
      f" of {len(many)}  <- the shortfall is collisions")
# Useful for streaming, where the category set is not known in advance.

# --- PUT IT IN A PIPELINE so encodings are fitted on train only. ---
prep = ColumnTransformer([
    ("nominal", OneHotEncoder(handle_unknown="ignore", sparse_output=False),
     ["colour"]),
    ("ordinal", OrdinalEncoder(categories=[SIZE_ORDER]), ["size"]),
], remainder="drop")
model = Pipeline([("prep", prep), ("lm", LinearRegression())])
model.fit(df[["colour", "size"]], df["price"])
print("\nfeature names out:", list(prep.get_feature_names_out()))

  CHOOSING AN ENCODING

  cardinality   ordered?   linear model        tree model
  --------------------------------------------------------------------
  2             n/a        one binary column   one binary column
  3-15          no         one-hot (drop one)  one-hot or native categorical
  3-15          yes        ordinal, real order ordinal, real order
  15-100        no         target/frequency    native categorical, or target
  100+          no         target + smoothing  native categorical, or hashing
  unbounded     no         hashing             hashing

  THE FABRICATED ORDER, DRAWN

    label encoding of {red, green, blue} as {0, 1, 2}
    implies:   red ---- green ---- blue
               |<-- 1 -->|<-- 1 -->|
    and therefore "the average of red and blue is green", plus a single
    coefficient forced to be monotone across all three. All invented.

    one-hot instead:  each category is its own axis, no distances implied.

  THE DUMMY TRAP
    red + green + blue == 1 for every row, so with an intercept the design
    matrix is singular. Drop one level for linear models with an intercept.

Key Takeaways

• Label encoding a nominal variable invents an order and a distance that a linear model will use as if they were real.
• One-hot encoding asserts nothing about order, at a cost of one column per category.
• Drop one level for linear models with an intercept, and always set explicit behaviour for unseen categories.
• Pass the real order to an ordinal encoder, because the default is alphabetical and alphabetical is not an order.
• Hashing gives fixed width without a fitted vocabulary and pays for it in collisions, which makes it the tool for streaming category sets.

🧪 Practice

1. Fit a linear regression on a label-encoded nominal variable and on its one-hot version, then compare the fitted values and explain the gap.
2. Ordinal-encode a size column with the default alphabetical order and with the true order, and show how the fitted coefficient changes.
3. One-hot encode a column, then transform a frame containing an unseen category under both ignore and error, and state which you would ship.
4. Interview: A feature has 50,000 distinct values. Walk through your options. (Hint: what breaks first, memory or the model's ability to estimate a parameter per level?)

Interaction and Polynomial Features

A linear model adds up the effects of its features. That is a strong assumption: it says the effect of one feature does not depend on any other, and it says the effect is constant across the feature's range. Both are often false, and interaction and polynomial features are how you relax them without abandoning the model.

An interaction is the product of two features, and it lets one feature's effect depend on another's. The canonical example is medicine: a drug helps younger patients and harms older ones. A model with separate age and dose terms cannot express that, because it must pick one sign for dose. Add age × dose and it can. The interaction coefficient is precisely the amount by which the dose effect changes per year of age.

A polynomial feature is a power of a single feature, and it lets the effect curve. Adding x² lets a model represent a relationship that rises then falls, such as the yield of a process against temperature.

The critical insight about both is that the model is still linear in its parameters. You are not changing the algorithm, you are changing the input space. This is why you can get non-linear behaviour out of ordinary least squares, and it is the same idea a kernel method implements implicitly.

Three practical cautions.

Combinatorial explosion. Degree-two polynomial features on 100 inputs give roughly 5,000 terms; degree three gives about 170,000. You almost never want a blanket expansion. Generate interactions you have a reason to believe in, or restrict to interaction_only degree two over a small chosen subset.

Collinearity. x and x² are highly correlated over a positive range, which makes coefficients unstable and hard to interpret. Centring the feature before squaring reduces this a lot, and is worth doing routinely.

Interpretation changes. Once an interaction is present, the coefficient on dose is no longer "the effect of dose". It is the effect of dose when age is zero, which may be a value that never occurs. Centring makes it "the effect at the mean age", which is interpretable. Forgetting this produces confident statements about coefficients that mean something else entirely.

Finally, note that tree ensembles learn interactions implicitly by splitting on one feature inside a branch defined by another. You rarely need to hand-build interactions for gradient boosting. You very often need to for linear and logistic regression, and that is precisely where a well-chosen interaction adds the most.

import numpy as np
import pandas as pd
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score

rng = np.random.default_rng(11)
n = 600
age = rng.integers(20, 80, n).astype(float)
dose = rng.uniform(0, 10, n)
# TRUE model: the dose effect flips sign with age. No additive model can
# represent this, because it must choose one sign for the dose term.
effect = 8.0 - 0.18 * age                      # positive young, negative old
y = 50 + 0.2 * age + effect * dose + rng.normal(0, 6, n)
X = pd.DataFrame({"age": age, "dose": dose})

additive = LinearRegression().fit(X, y)
print("additive model:")
print(f"  R2 {additive.score(X, y):.3f}  dose coef {additive.coef_[1]:+.2f}")
print("  a single dose coefficient averages a helpful and a harmful effect")

Xi = X.assign(age_x_dose=X["age"] * X["dose"])
inter = LinearRegression().fit(Xi, y)
print("\nwith the interaction:")
print(f"  R2 {inter.score(Xi, y):.3f}")
for name, c in zip(Xi.columns, inter.coef_):
    print(f"  {name:12s} {c:+.3f}")
print(f"  recovered dose effect at age 40:"
      f" {inter.coef_[1] + 40 * inter.coef_[2]:+.2f}  (true {8 - 0.18*40:+.2f})")

# --- CENTRING MAKES THE MAIN EFFECT INTERPRETABLE ---
Xc = pd.DataFrame({"age_c": age - age.mean(), "dose_c": dose - dose.mean()})
Xc["inter"] = Xc["age_c"] * Xc["dose_c"]
centred = LinearRegression().fit(Xc, y)
print(f"\nuncentred dose coef {inter.coef_[1]:+.2f} = effect at age 0 (impossible)")
print(f"centred   dose coef {centred.coef_[1]:+.2f} = effect at the MEAN age")

# It also cuts the collinearity that makes coefficients unstable.
print(f"\ncorr(age, age*dose)     {np.corrcoef(age, age * dose)[0,1]:+.3f}")
print(f"corr(age_c, centred int) "
      f"{np.corrcoef(Xc['age_c'], Xc['inter'])[0,1]:+.3f}")

# --- POLYNOMIAL: curvature in one feature. ---
temp = rng.uniform(80, 220, n)
yield_ = 100 - 0.01 * (temp - 150) ** 2 + rng.normal(0, 2, n)   # peaks at 150

lin = LinearRegression().fit(temp.reshape(-1, 1), yield_)
quad = Pipeline([("poly", PolynomialFeatures(2, include_bias=False)),
                 ("lm", LinearRegression())]).fit(temp.reshape(-1, 1), yield_)
print(f"\nlinear R2 {lin.score(temp.reshape(-1,1), yield_):.3f}"
      f"   quadratic R2 {quad.score(temp.reshape(-1,1), yield_):.3f}")

coefs = quad.named_steps["lm"].coef_
peak = -coefs[0] / (2 * coefs[1])
print(f"quadratic implies an optimum at {peak:.0f} degrees (true 150)")

# --- THE EXPLOSION: count terms before you generate them. ---
print("\nfeatures produced by PolynomialFeatures:")
rows = []
for p in (5, 20, 100):
    dummy = np.zeros((2, p))
    full2 = PolynomialFeatures(2, include_bias=False).fit(dummy)
    only2 = PolynomialFeatures(2, include_bias=False,
                               interaction_only=True).fit(dummy)
    full3 = PolynomialFeatures(3, include_bias=False).fit(dummy)
    rows.append({"inputs": p,
                 "degree2": full2.n_output_features_,
                 "degree2_interaction_only": only2.n_output_features_,
                 "degree3": full3.n_output_features_})
print(pd.DataFrame(rows).to_string(index=False))

# --- HIGH DEGREE OVERFITS; regularisation and honest CV are the defences. ---
x1 = rng.uniform(-3, 3, 120)
y1 = np.sin(x1) + rng.normal(0, 0.25, 120)
print("\ncross-validated R2 by polynomial degree:")
for deg in (1, 3, 5, 9, 15):
    plain = Pipeline([("p", PolynomialFeatures(deg, include_bias=False)),
                      ("s", StandardScaler()), ("m", LinearRegression())])
    ridged = Pipeline([("p", PolynomialFeatures(deg, include_bias=False)),
                       ("s", StandardScaler()), ("m", Ridge(alpha=10.0))])
    a = cross_val_score(plain, x1.reshape(-1, 1), y1, cv=5, scoring="r2").mean()
    b = cross_val_score(ridged, x1.reshape(-1, 1), y1, cv=5, scoring="r2").mean()
    print(f"  degree {deg:2d}   plain {a:+.3f}   ridge {b:+.3f}")
# Note the plain model's CV score collapsing at high degree while the
# in-sample fit keeps improving. That gap IS overfitting.

# --- DOMAIN-DRIVEN FEATURES BEAT BLIND EXPANSION ---
listings = pd.DataFrame({
    "sqft": [1200, 2400, 900], "bedrooms": [2, 4, 1],
    "lot_sqft": [4000, 6000, 2500], "price": [300_000, 620_000, 210_000],
})
listings["sqft_per_bedroom"] = listings["sqft"] / listings["bedrooms"]
listings["building_ratio"] = listings["sqft"] / listings["lot_sqft"]
listings["price_per_sqft"] = listings["price"] / listings["sqft"]
print("\nratios a blanket polynomial expansion would never produce:")
print(listings.round(1).to_string(index=False))
# A ratio is a quotient, not a product. Domain knowledge supplies the form.

# --- TREES FIND INTERACTIONS THEMSELVES ---
from sklearn.ensemble import GradientBoostingRegressor
print("\ngradient boosting on raw features only:",
      round(cross_val_score(GradientBoostingRegressor(random_state=0),
                            X, y, cv=5, scoring="r2").mean(), 3))
print("linear with explicit interaction:      ",
      round(cross_val_score(LinearRegression(), Xi, y, cv=5,
                            scoring="r2").mean(), 3))

  WHAT AN INTERACTION BUYS YOU

  additive model:  y = a + b*age + c*dose
    the dose effect is ONE number, the same at every age.

      effect of dose
        |______________________   c
        +----------------------> age

  with interaction: y = a + b*age + c*dose + d*(age*dose)
    the dose effect is a LINE in age, slope d.

      effect of dose
        |\
        | \________              c + d*age
        |          \
        +-----------\----------> age
                     crosses zero: helpful below, harmful above

  TERM COUNT

    inputs   degree 2   degree 2 (interactions only)   degree 3
    5        20         15                             55
    20       230        210                            1,770
    100      5,150      5,050                          176,850

  INTERPRETATION AFTER AN INTERACTION

    uncentred: coef(dose) = effect of dose when age == 0        (meaningless)
    centred:   coef(dose) = effect of dose at the MEAN age      (usable)

Key Takeaways

• An interaction term lets one feature's effect depend on another's, which a purely additive model cannot express at any coefficient value.
• Polynomial terms let the effect curve, and the model stays linear in its parameters throughout.
• Centre features before forming products and powers, both to reduce collinearity and to keep the main effects interpretable.
• Blanket expansion explodes combinatorially and overfits; pick terms you have a reason to believe in, and regularise.
• Tree ensembles discover interactions by splitting, so hand-built interactions pay off mainly for linear and logistic models.

🧪 Practice

1. Simulate data where a feature's effect reverses sign across another feature's range, then show that an additive model cannot fit it and one interaction term can.
2. Fit a quadratic to a curved relationship and recover the location of the optimum from the fitted coefficients.
3. Compare cross-validated scores for polynomial degrees 1 through 15, with and without ridge regularisation, and identify where overfitting begins.
4. Interview: You add a product term and the main effect's coefficient flips sign. Has the relationship changed? (Hint: what value of the other feature is the main effect now measured at?)

Date and Time Features

A timestamp is one of the richest columns you will ever get and one of the least useful in raw form. As a number of seconds since 1970 it is a monotonically increasing value that a model can only use as a trend, and the first thing it will learn is that larger timestamps mean later, which is not a pattern that generalises to the future at all.

The work is decomposing it into the components that actually carry signal.

Calendar components. Hour of day, day of week, day of month, week of year, month, quarter, year. These capture the cycles that drive most human-generated data: the commute peak, the weekend, the month-end close, the December holidays.

Boolean flags. Weekend, holiday, month start, month end, quarter end. Holidays matter disproportionately and cannot be derived from the calendar alone, because they depend on jurisdiction and some of them move. Use a maintained holiday calendar.

Elapsed durations. Days since signup, hours since last login, time until contract renewal. These are frequently the strongest features in the whole model, because a behavioural state usually depends on how long since something happened rather than on the absolute date.

Cyclical encoding. This is the one people miss. Hour 23 and hour 0 are adjacent in reality and maximally distant as integers. Encoding an hour as a single number tells the model that late night is far from early morning. The fix is to map the cycle onto a circle with a sine and a cosine pair, so that adjacency in time becomes adjacency in the feature space. For trees this matters less, because a tree can split out hour 23 on its own, but for linear models and neural networks it matters a lot.

Time zones. If your users are spread across zones, "9am" in UTC is breakfast for some and midnight for others. Local time is what drives behaviour, so store UTC and derive local-time features using each entity's zone. Getting this wrong flattens every daily pattern in your data.

And the leakage point, which is the most important one: a feature must only use information available at prediction time. "Days until churn" is the label. "Total sessions in the user's lifetime" computed over the full dataset includes sessions that happen after the prediction moment. Every temporal feature needs an explicit as-of timestamp, and the discipline of writing that timestamp down catches most leakage before it happens.

import numpy as np
import pandas as pd

ts = pd.Series(pd.to_datetime([
    "2024-01-01 00:15", "2024-03-15 09:30", "2024-07-04 18:45",
    "2024-11-29 23:50", "2024-12-31 12:00",
]), name="event_at")

# --- CALENDAR DECOMPOSITION ---
feat = pd.DataFrame({
    "hour": ts.dt.hour,
    "dayofweek": ts.dt.dayofweek,            # Monday = 0
    "day": ts.dt.day,
    "month": ts.dt.month,
    "quarter": ts.dt.quarter,
    "weekofyear": ts.dt.isocalendar().week.astype(int),
    "is_weekend": ts.dt.dayofweek.isin([5, 6]),
    "is_month_end": ts.dt.is_month_end,
    "is_quarter_end": ts.dt.is_quarter_end,
    "days_in_month": ts.dt.days_in_month,
})
print(pd.concat([ts, feat], axis=1).to_string(index=False))

# --- THE RAW TIMESTAMP AS A FEATURE: a trend that cannot extrapolate. ---
print("\nepoch seconds:", (ts.astype("int64") // 10**9).tolist())
print("a model fitted on these learns 'later = bigger', which says nothing")
print("about next month, because every future value is out of range.")

# --- CYCLICAL ENCODING: hour 23 and hour 0 are neighbours. ---
def cyclical(values: pd.Series, period: int, name: str) -> pd.DataFrame:
    theta = 2 * np.pi * values / period
    return pd.DataFrame({f"{name}_sin": np.sin(theta),
                         f"{name}_cos": np.cos(theta)})

hours = pd.Series([0, 6, 12, 18, 23])
cyc = cyclical(hours, 24, "hour").round(3)
print("\n", pd.concat([hours.rename("hour"), cyc], axis=1).to_string(index=False))

def dist(a, b, frame):
    pa, pb = frame.loc[a], frame.loc[b]
    return float(np.hypot(pa.iloc[0] - pb.iloc[0], pa.iloc[1] - pb.iloc[1]))

print(f"\nraw distance  hour 23 to 0: {abs(23 - 0)}   hour 11 to 12: {abs(11-12)}")
print(f"cyclic distance 23 to 0: {dist(4, 0, cyc):.3f}"
      f"   12 to 18: {dist(2, 3, cyc):.3f}")
# On the circle, 23 and 0 are as close as any other one-hour step. On the
# integer scale they are 23 apart, which is simply false.

# Apply the same trick to any cycle: weekday period 7, month period 12.
month_cyc = cyclical(ts.dt.month, 12, "month")

# --- HOLIDAYS: not derivable from the calendar. ---
from pandas.tseries.holiday import USFederalHolidayCalendar
cal = USFederalHolidayCalendar()
holidays = cal.holidays(start="2024-01-01", end="2024-12-31")
hol_feat = pd.DataFrame({
    "date": ts.dt.normalize(),
    "is_holiday": ts.dt.normalize().isin(holidays),
})
# Proximity matters as much as the day itself: retail moves before a holiday.
nearest = ts.dt.normalize().map(
    lambda d: min((abs((d - h).days) for h in holidays), default=np.nan))
hol_feat["days_to_nearest_holiday"] = nearest
print("\n", hol_feat.to_string(index=False))

# --- ELAPSED TIME: usually the strongest temporal feature. ---
users = pd.DataFrame({
    "user": ["u1", "u2", "u3"],
    "signup": pd.to_datetime(["2023-01-15", "2024-02-01", "2024-05-20"]),
    "last_login": pd.to_datetime(["2024-05-01", "2024-05-25", None]),
})
AS_OF = pd.Timestamp("2024-06-01")            # the prediction moment, stated
users["tenure_days"] = (AS_OF - users["signup"]).dt.days
users["days_since_login"] = (AS_OF - users["last_login"]).dt.days
users["never_logged_in"] = users["last_login"].isna()
users["logins_per_tenure_month"] = (30 / users["days_since_login"]).round(2)
print("\n", users.to_string(index=False))
# Naming the as-of timestamp as a constant is the single best defence
# against temporal leakage: every duration is measured from one stated point.

# --- TIME ZONES: local time drives behaviour. ---
events = pd.DataFrame({
    "utc": pd.to_datetime(["2024-06-01 13:00", "2024-06-01 13:00"], utc=True),
    "tz": ["America/New_York", "Asia/Tokyo"],
})
events["local_hour"] = [
    t.tz_convert(z).hour for t, z in zip(events["utc"], events["tz"])]
print("\n", events.to_string(index=False))
print("the same UTC hour is morning in one zone and late night in the other;")
print("using the UTC hour as a feature averages those away.")

# --- LEAKAGE: a lifetime aggregate computed over the whole table. ---
log = pd.DataFrame({
    "user": ["u1"] * 4,
    "at": pd.to_datetime(["2024-05-10", "2024-05-20", "2024-06-10", "2024-06-20"]),
})
print(f"\nlifetime session count (LEAKS): {len(log)}")
print(f"count as of {AS_OF.date()} (safe): {(log['at'] < AS_OF).sum()}")
# The first number includes June rows. A model predicting on 1 June could
# not have known them, and validation using them will look great and lie.

# --- RESAMPLING: aligning irregular events to a regular grid. ---
irregular = pd.Series([1, 1, 1, 1], index=pd.to_datetime(
    ["2024-06-01 09:00", "2024-06-01 09:40", "2024-06-03 11:00",
     "2024-06-03 11:30"]), name="events")
daily = irregular.resample("D").sum()          # absent days become 0, correctly
print("\ndaily counts:\n", daily.to_string())

  WHY CYCLICAL ENCODING

  hour as an integer                 hour on a circle
  0  1  2 ............ 22  23                 0
  |--|--|--------------|--|              23        1
  0 and 23 are 23 apart                22            2
                                        |            |
                                       18     ...    6
                                          12
  0 and 23 are one step apart, which is what the clock means.

  FEATURE FAMILIES FROM ONE TIMESTAMP

  family        examples                            strongest for
  ------------------------------------------------------------------------
  calendar      hour, dow, month, quarter           seasonality
  flags         weekend, holiday, month_end          calendar effects
  cyclical      sin/cos of hour, dow, month          linear models, NNs
  elapsed       days since signup / last event       behavioural state
  counts as-of  events in the last 7 / 30 days       recent intensity
  local time    hour in the entity's own time zone   anything human

  THE AS-OF RULE

    AS_OF = the moment the prediction is made.
    Every feature answers: "what was true strictly BEFORE AS_OF?"
    A feature that cannot be written that way is a label in disguise.

Key Takeaways

• A raw timestamp is a trend a model cannot extrapolate; the value is in the decomposition.
• Encode cycles with a sine and cosine pair so that hour 23 and hour 0 are neighbours, which matters most for linear models and neural networks.
• Elapsed durations such as days since last login are often the strongest temporal features, and they need an explicit as-of timestamp.
• Holidays require a maintained calendar and proximity to a holiday often matters as much as the day itself.
• Derive behavioural features in the entity's local time while storing UTC, or daily patterns average away.

🧪 Practice

1. Decompose a timestamp column into at least eight features and report which of them vary in your data and which are constant.
2. Encode hour of day cyclically and verify numerically that hours 23 and 0 are closer than hours 0 and 12.
3. Compute a "sessions in the last 30 days" feature for a fixed as-of date and show that the naive lifetime count differs.
4. Interview: A churn model uses "days since last login" and performs nearly perfectly in validation. Why are you suspicious? (Hint: what is the last login of a user who has already churned, and when was it measured relative to the label?)

Text Feature Extraction

Text is a sequence of characters and a model needs a vector, so something has to bridge that gap. The methods form a ladder, and the useful thing is knowing where on the ladder your problem sits, because the cheapest rung that works is almost always the right answer.

The bottom rung is manual features: length in characters, word count, count of digits, proportion of capitals, presence of a URL, count of exclamation marks. These sound trivial and they are unreasonably effective for tasks like spam detection, because the signal in "ALL CAPS WITH!!! LINKS" is mostly not in the words.

The next rung is bag of words, which counts token occurrences and discards order. "The dog bit the man" and "The man bit the dog" produce identical vectors, which is either an acceptable simplification or a fatal flaw depending on the task. For topic classification it is acceptable. For sentiment it is often not, which is why you add n-grams: contiguous token pairs or triples, so that "not good" survives as a unit.

TF-IDF adjusts the counts. The intuition is that a word appearing in every document tells you nothing about which document you are looking at, while a word appearing in three of ten thousand documents is highly informative. Term frequency measures how often the word appears here; inverse document frequency downweights words that appear everywhere. The product is the standard representation for classical text classification and retrieval, and it is a strong baseline that people skip past too quickly.

Before any of this comes preprocessing, and each step is a deliberate loss of information. Lowercasing merges "Apple" and "apple", which is right for topic detection and wrong for named-entity work. Removing stop words shrinks the vocabulary and destroys negation. Stemming chops suffixes crudely ("running" to "run", but also "university" to "univers"); lemmatisation does it properly with a dictionary and costs more. Choose each based on the task, and be aware that defaults are choices too.

Character n-grams deserve a mention because they solve problems word tokens cannot: misspellings, languages without spaces, product codes, and fuzzy name matching. They produce more features and are remarkably robust.

The top rung is embeddings: dense vectors where semantic similarity becomes geometric proximity, so that "cat" and "kitten" are close without ever appearing in the same document. Word-level embeddings such as word2vec and GloVe give one vector per word type. Contextual and sentence-level embeddings from transformer models give a vector per occurrence or per document, capture word sense, and are the current default when accuracy matters more than cost and interpretability.

The practical guidance is simple: start with TF-IDF and a linear model. It is fast, it is interpretable feature by feature, and it is frequently within a few points of anything more elaborate. Move up the ladder when you have measured that you need to.

import numpy as np
import pandas as pd
from sklearn.feature_extraction.text import CountVectorizer, TfidfVectorizer
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

docs = pd.Series([
    "Free money!!! Click here NOW http://spam.example",
    "Meeting moved to 3pm, see the attached agenda",
    "URGENT: your account will be CLOSED, verify at http://x.example",
    "Can you review the quarterly report before Friday?",
    "Win a FREE prize, claim now!!!",
    "Lunch tomorrow? I booked the usual place",
], name="text")
labels = pd.Series([1, 0, 1, 0, 1, 0], name="spam")

# --- RUNG 1: manual features. Cheap, interpretable, often enough. ---
manual = pd.DataFrame({
    "n_chars": docs.str.len(),
    "n_words": docs.str.split().str.len(),
    "n_exclaim": docs.str.count("!"),
    "pct_upper": (docs.str.count(r"[A-Z]") / docs.str.len()).round(3),
    "has_url": docs.str.contains(r"https?://").astype(int),
    "n_digits": docs.str.count(r"\d"),
    "avg_word_len": (docs.str.replace(r"\s", "", regex=True).str.len()
                     / docs.str.split().str.len()).round(2),
})
print(pd.concat([labels, manual], axis=1).to_string(index=False))
print("\ncorrelation with the label:")
print(manual.corrwith(labels).round(2).sort_values(ascending=False).to_string())
# pct_upper and n_exclaim alone separate these classes. No vocabulary needed.

# --- RUNG 2: bag of words. Order is discarded. ---
cv = CountVectorizer(lowercase=True, min_df=1)
counts = cv.fit_transform(docs)
print(f"\nbag of words: {counts.shape[1]} features,"
      f" density {counts.nnz / np.prod(counts.shape):.1%}")

a = cv.transform(["the dog bit the man"]).toarray()
b = cv.transform(["the man bit the dog"]).toarray()
print(f"'dog bit man' and 'man bit dog' identical? {np.array_equal(a, b)}")

# N-GRAMS restore the local order that matters most.
bigrams = CountVectorizer(ngram_range=(1, 2))
bigrams.fit(["the service was not good", "the service was good"])
print("\nvocabulary with bigrams:", sorted(bigrams.vocabulary_)[:8], "...")
print("'not good' exists as one token, so negation survives")

# --- RUNG 3: TF-IDF. Rare words weigh more. ---
tfidf = TfidfVectorizer(lowercase=True, sublinear_tf=True, min_df=1)
M = tfidf.fit_transform(docs)
vocab = np.array(tfidf.get_feature_names_out())
idf = pd.Series(tfidf.idf_, index=vocab).sort_values()
print("\nlowest IDF (common, uninformative):", idf.head(4).round(2).to_dict())
print("highest IDF (rare, informative):   ", idf.tail(4).round(2).to_dict())

# The top-weighted terms of one document are a readable explanation.
row = pd.Series(M[0].toarray().ravel(), index=vocab)
print("\ntop terms in document 0:", row.nlargest(4).round(3).to_dict())

# --- PREPROCESSING CHOICES ARE INFORMATION LOSS. Make them deliberately. ---
sample = ["The service was NOT good at all", "The service was good"]
configs = {
    "defaults": TfidfVectorizer(),
    "no stop words removed": TfidfVectorizer(stop_words=None, ngram_range=(1, 2)),
    "stop words removed": TfidfVectorizer(stop_words="english"),
    "char 3-5 grams": TfidfVectorizer(analyzer="char_wb", ngram_range=(3, 5)),
}
print("\nfeature count and whether the two sentences stay distinguishable:")
for name, vec in configs.items():
    V = vec.fit_transform(sample)
    sim = float((V[0] @ V[1].T).toarray().ravel()[0])
    print(f"  {name:24s} {V.shape[1]:5d} features   cosine {sim:.3f}")
# "stop words removed" drops "not", so the two opposite sentences become
# nearly identical. That is a correctness bug, not a tuning choice.

# --- CHARACTER N-GRAMS survive misspelling. ---
pair = ["accommodation available", "acommodation avalable"]   # typos
word_v = TfidfVectorizer(analyzer="word").fit(pair)
char_v = TfidfVectorizer(analyzer="char_wb", ngram_range=(3, 4)).fit(pair)
Wm, Cm = word_v.transform(pair), char_v.transform(pair)
print(f"\nmisspelled pair cosine: word tokens "
      f"{float((Wm[0] @ Wm[1].T).toarray().ravel()[0]):.3f}"
      f"   char n-grams {float((Cm[0] @ Cm[1].T).toarray().ravel()[0]):.3f}")

# --- THE BASELINE TO BEAT: TF-IDF plus a linear model, in a pipeline. ---
pipe = Pipeline([
    ("tfidf", TfidfVectorizer(ngram_range=(1, 2), min_df=1, sublinear_tf=True)),
    ("clf", LogisticRegression(max_iter=1000)),
])
print("\nCV accuracy (tiny sample, illustrative):",
      round(cross_val_score(pipe, docs, labels, cv=3).mean(), 3))
# CRITICAL: the vectorizer is inside the pipeline, so the vocabulary and the
# IDF weights are fitted per fold. Fitting the vectorizer on all documents
# before splitting leaks the test set's vocabulary and document frequencies.

pipe.fit(docs, labels)
weights = pd.Series(pipe.named_steps["clf"].coef_.ravel(),
                    index=pipe.named_steps["tfidf"].get_feature_names_out())
print("\nmost spam-indicative terms:", weights.nlargest(4).round(2).to_dict())
print("most ham-indicative terms: ", weights.nsmallest(4).round(2).to_dict())
# This readable weight list is why TF-IDF plus linear remains valuable even
# when a transformer scores higher.

# --- DIMENSION CONTROL on real corpora ---
print("\nvocabulary controls that matter at scale:")
print("  min_df=5        drop terms in fewer than 5 documents (typos, noise)")
print("  max_df=0.5      drop terms in over half the documents (boilerplate)")
print("  max_features=N  keep the N most frequent terms")
print("  HashingVectorizer: fixed width, no vocabulary, streaming-friendly")

# --- RUNG 4: embeddings, for when meaning must transfer across words. ---
print("\nTF-IDF cosine between 'cat sat on the mat' and 'kitten rested on rug':")
z = TfidfVectorizer().fit_transform(["cat sat on the mat", "kitten rested on rug"])
print(f"  {float((z[0] @ z[1].T).toarray().ravel()[0]):.3f}"
      "  <- no shared content words, so near zero")
print("  A sentence embedding scores these as similar, because it encodes")
print("  meaning rather than surface tokens. That is the reason to climb.")

  THE LADDER

  rung  representation       captures            cost     interpretable
  ---------------------------------------------------------------------------
  1     manual features      surface form        trivial  fully
  2     bag of words         vocabulary          low      per token
  3     + n-grams            short phrases       medium   per phrase
  4     TF-IDF               informativeness     low      per token, weighted
  5     char n-grams         morphology, typos   medium   partly
  6     word embeddings      word similarity     medium   no
  7     contextual/sentence  word sense, context high     no

  Start at 4. Climb only when a measurement says to.

  WHAT PREPROCESSING COSTS

  step             merges                     destroys
  ---------------------------------------------------------------
  lowercase        Apple / apple              proper nouns, acronyms
  remove stop wd   shrinks vocabulary         negation ("not good")
  stemming         run/running/ran            readability of features
  lemmatisation    run/running/ran (properly) nothing, costs compute
  strip punctuation                           "!!!", emoticons, URLs

  TF-IDF INTUITION

    tf(t,d)  = how often t occurs in d          -> local importance
    idf(t)   = log(N / documents containing t)  -> global rarity
    "the" is in every document -> idf ~ 0 -> weight ~ 0, automatically.

Key Takeaways

• Manual surface features are unreasonably strong for tasks where form rather than meaning carries the signal.
• Bag of words discards order; n-grams restore the local order that matters, such as negation.
• TF-IDF downweights terms that appear everywhere and is the baseline any more expensive method has to beat.
• Every preprocessing step is deliberate information loss, and removing stop words destroys negation.
• Fit the vectorizer inside the pipeline, or the test set's vocabulary and document frequencies leak into training.

🧪 Practice

1. Build five manual text features and report their correlation with a binary label on any labelled text you have.
2. Show that a unigram vectorizer maps two sentences differing only in "not" to nearly identical vectors, and that adding bigrams separates them.
3. Compare word-token and character n-gram similarity for a misspelled pair of strings.
4. Interview: Your TF-IDF model scores 0.94 in cross-validation and 0.71 in production. Name two causes specific to text. (Hint: one is about where the vocabulary came from, the other about how language changes over time.)

Target and Frequency Encoding

One-hot encoding a variable with ten thousand categories gives you ten thousand columns, most of which are almost always zero. That is unusable, so you need an encoding whose width does not grow with cardinality. The two main options both replace a category with a statistic about it.

Frequency encoding replaces each category with how often it occurs. It is unsupervised, so it cannot leak the target, and it is surprisingly useful: rarity is often predictive in its own right, since rare merchant codes, rare user agents, and rare error messages all tend to behave differently from common ones. Its weakness is that two categories with the same count become indistinguishable.

Target encoding, also called mean encoding or likelihood encoding, replaces each category with the mean of the target within that category. A postcode becomes its historical default rate. This is enormously powerful because it compresses all the information the category carries about the target into one number, and it is enormously dangerous for exactly the same reason: you have put the answer into the feature.

Two things go wrong, and both have standard fixes.

Overfitting on rare categories. A category with one row gets that row's target as its encoded value, which is a perfect predictor of that row and pure noise for any other. The fix is smoothing: blend the category mean toward the global mean, weighted by how many rows the category has. With many rows, trust the category; with few, fall back to the prior. The weight is usually n / (n + m) for a smoothing constant m, which is a simple form of empirical Bayes shrinkage.

Leakage through the row's own target. If a row's encoded value was computed using that row's own label, the feature contains the answer. Cross-validated or out-of-fold encoding fixes this: compute each row's encoding from the other folds only, so no row ever sees its own target. The version of this used by CatBoost goes further, computing each row's encoding from rows that precede it in a random permutation, which removes the leakage almost entirely.

The failure mode to recognise is distinctive: a target-encoded model shows an excellent validation score and a much worse test score, and the gap grows with the cardinality of the encoded variable. If you see that, suspect the encoding before anything else.

There are variants worth knowing. Count-and-sum encoding keeps the numerator and denominator as separate features, letting the model combine them itself. Weight-of-evidence encoding is the credit-scoring standard, using a log-odds ratio and inheriting the same smoothing needs. For multiclass problems you need one encoded column per class, which starts to cost width again.

import numpy as np
import pandas as pd
from sklearn.model_selection import KFold, cross_val_score, train_test_split
from sklearn.ensemble import GradientBoostingClassifier

rng = np.random.default_rng(9)
n = 4000
# 300 categories with very uneven support: the long tail is where it breaks.
cats = rng.choice([f"z{i:03d}" for i in range(300)], n,
                  p=(w := rng.dirichlet(np.full(300, 0.3))))
true_rate = dict(zip([f"z{i:03d}" for i in range(300)],
                     rng.beta(2, 5, 300)))
y = pd.Series((rng.random(n) < np.array([true_rate[c] for c in cats])).astype(int))
X = pd.DataFrame({"cat": cats})
print(f"{X['cat'].nunique()} categories, "
      f"{(X['cat'].value_counts() <= 3).sum()} of them with 3 rows or fewer")

# --- FREQUENCY ENCODING: unsupervised, so it cannot leak. ---
freq = X["cat"].value_counts(normalize=True)
X["cat_freq"] = X["cat"].map(freq)
X["cat_count"] = X["cat"].map(X["cat"].value_counts())
print(f"\nfrequency encoding: 1 column, correlation with y "
      f"{X['cat_freq'].corr(y):+.3f}")

# --- NAIVE TARGET ENCODING: the leak, demonstrated. ---
naive_map = y.groupby(X["cat"]).mean()
X["te_naive"] = X["cat"].map(naive_map)
print(f"\nnaive target encoding correlation with y: {X['te_naive'].corr(y):+.3f}")

singles = X["cat"].map(X["cat"].value_counts()) == 1
print(f"rows in singleton categories: {singles.sum()}")
print("for those rows the encoded value EQUALS their own label:",
      bool((X.loc[singles, 'te_naive'] == y[singles]).all()))
# The feature contains the answer for exactly the rows where it helps least.

# --- SMOOTHING: shrink small categories toward the global mean. ---
def smooth_target_encode(keys: pd.Series, target: pd.Series, m: float = 20.0):
    """Empirical-Bayes shrinkage: weight n/(n+m) on the category mean."""
    prior = target.mean()
    stats = target.groupby(keys).agg(["sum", "count"])
    smoothed = (stats["sum"] + m * prior) / (stats["count"] + m)
    return smoothed, prior

mapping, prior = smooth_target_encode(X["cat"], y, m=20.0)
X["te_smooth"] = X["cat"].map(mapping).fillna(prior)

show = pd.DataFrame({
    "count": y.groupby(X["cat"]).count(),
    "raw_mean": y.groupby(X["cat"]).mean().round(3),
    "smoothed": mapping.round(3),
}).sort_values("count")
print(f"\nglobal prior {prior:.3f}")
print("\nsmallest categories are pulled hard toward the prior:")
print(show.head(4).to_string())
print("\nlargest categories keep their own mean:")
print(show.tail(3).to_string())

# --- OUT-OF-FOLD ENCODING: no row sees its own target. ---
def oof_target_encode(keys: pd.Series, target: pd.Series,
                      n_splits: int = 5, m: float = 20.0, seed: int = 0):
    out = pd.Series(np.nan, index=keys.index, dtype=float)
    kf = KFold(n_splits=n_splits, shuffle=True, random_state=seed)
    for fit_idx, apply_idx in kf.split(keys):
        mp, pr = smooth_target_encode(keys.iloc[fit_idx], target.iloc[fit_idx], m)
        out.iloc[apply_idx] = keys.iloc[apply_idx].map(mp).fillna(pr).values
    return out

X["te_oof"] = oof_target_encode(X["cat"], y)
print("\ncorrelation with the target:")
print(pd.Series({
    "frequency (safe)": X["cat_freq"].corr(y),
    "naive TE (leaks)": X["te_naive"].corr(y),
    "smoothed TE (leaks less)": X["te_smooth"].corr(y),
    "out-of-fold TE (safe)": X["te_oof"].corr(y),
}).round(3).to_string())
# The naive correlation is inflated by construction. The out-of-fold value
# is the honest estimate of how much the category actually tells you.

# --- THE VALIDATION GAP, MEASURED ---
Xtr, Xte, ytr, yte = train_test_split(X[["cat"]], y, test_size=0.3, random_state=1)

# WRONG: encode with a map fitted on everything, then split.
bad_map = y.groupby(X["cat"]).mean()
bad_tr = Xtr["cat"].map(bad_map).to_frame("te")
bad_te = Xte["cat"].map(bad_map).fillna(y.mean()).to_frame("te")

# RIGHT: out-of-fold on train, and a train-only map applied to test.
good_tr = oof_target_encode(Xtr["cat"], ytr).to_frame("te")
mp, pr = smooth_target_encode(Xtr["cat"], ytr, m=20.0)
good_te = Xte["cat"].map(mp).fillna(pr).to_frame("te")

clf = GradientBoostingClassifier(random_state=0)
for label, (a, b) in {"leaky encoding": (bad_tr, bad_te),
                      "out-of-fold encoding": (good_tr, good_te)}.items():
    clf.fit(a, ytr)
    print(f"\n{label}")
    print(f"  train-fit score {clf.score(a, ytr):.3f}"
          f"   held-out score {clf.score(b, yte):.3f}")
# The leaky version shows a wide train/test gap. The honest version's
# numbers are lower and are the ones that will hold in production.

# --- UNSEEN CATEGORIES fall back to the prior, always. ---
unseen = pd.Series(["brand_new_code"])
print(f"\nunseen category encodes to the prior: "
      f"{unseen.map(mp).fillna(pr).iloc[0]:.3f}")

# --- COUNT-AND-SUM: give the model the parts, let it combine them. ---
parts = y.groupby(X["cat"]).agg(["sum", "count"])
X["cat_pos"] = X["cat"].map(parts["sum"])
X["cat_n"] = X["cat"].map(parts["count"])
print("\ncount-and-sum keeps the sample size visible to the model, so it can")
print("distrust a 1.0 rate that rests on two rows.")

  SHRINKAGE: HOW MUCH TO TRUST A CATEGORY

  encoded = (sum_of_target + m * prior) / (count + m)

    count    weight on the category mean (m=20)
    1        0.05      almost entirely the prior
    5        0.20
    20       0.50      half and half
    100      0.83
    1000     0.98      almost entirely its own mean

  larger m = more conservative. Tune it; it is a real hyperparameter.

  THE THREE WAYS TARGET ENCODING LEAKS

  1. own row      the encoding used this row's label      -> out-of-fold
  2. whole set    the map was fitted before the split     -> fit on train only
  3. time         the map used rows after the prediction  -> fit on past only

  ALL THREE must be closed. Closing one is not enough.

  SYMPTOM TO RECOGNISE

    validation AUC 0.94, test AUC 0.71, and the gap widens with cardinality
    -> suspect the target encoding before the model.

Key Takeaways

• Frequency encoding is unsupervised and cannot leak, and rarity is often predictive on its own.
• Target encoding compresses a high-cardinality category into one strong feature and puts the answer into the input if done carelessly.
• Smooth toward the global mean with weight n / (n + m) so rare categories do not memorise their own rows.
• Compute encodings out of fold so no row ever sees its own target, and fit the final map on training rows only.
• A wide validation-to-test gap that grows with cardinality is the signature of leaky target encoding.

🧪 Practice

1. Target-encode a high-cardinality column naively and show that singleton categories receive their own row's label as the feature value.
2. Implement smoothed encoding and plot or tabulate the shrinkage weight against category size for two values of the smoothing constant.
3. Compare held-out scores for a model trained on a leaky encoding and one trained on out-of-fold encoding.
4. Interview: A teammate target-encodes before splitting and reports 0.97 AUC. What do you predict happens in production, and how do you show them without a deployment? (Hint: which rows contributed to their own feature's value?)

Feature Selection Methods

More features are not better. Past a point they cost accuracy, not just compute, and the reasons are worth understanding rather than taking on faith.

Irrelevant features add variance: the model fits noise in them, and the noise does not repeat. Redundant features destabilise coefficients, since two near-identical columns can trade arbitrarily large opposite weights and the model has no basis to choose. Sheer dimensionality thins the data: in high dimensions every point is far from every other, distances become uninformative, and the volume the model must cover grows exponentially. This last effect is the curse of dimensionality, and it is why a distance-based method degrades sharply as you add columns.

There are also operational costs. Every feature is a pipeline dependency, a potential source of drift, and something that can be null in production for a reason nobody anticipated. Fewer features is less to monitor.

Selection methods fall into three families.

Filter methods score each feature against the target independently of any model. Correlation, mutual information, chi-squared, ANOVA F-statistic. They are fast, they scale to very wide data, and they share a blind spot: they evaluate features one at a time, so they keep two redundant features that both score well and discard a feature that is useless alone and essential in combination. The XOR pattern is the standard counterexample, where each input has zero marginal association with the target and the pair determines it completely.

Wrapper methods train the model repeatedly on different feature subsets and keep what improves the score. Forward selection adds features one at a time; backward elimination removes them; recursive feature elimination drops the weakest feature and refits, repeatedly. These account for interactions because the model does, and they are expensive, and they overfit the selection criterion if the score you optimise is the same one you report.

Embedded methods select as part of fitting. Lasso's L1 penalty drives coefficients to exactly zero, which is selection. Tree ensembles produce importances. These are the best value in practice: one fit, interactions accounted for, and the selection is consistent with the model you will ship.

Two warnings about importances specifically. Impurity-based importance is biased toward high-cardinality and continuous features, because they offer more places to split, and it can rank a pure-noise column with many distinct values above a genuinely useful binary one. Permutation importance is measured on held-out data by shuffling one column and observing the drop in score, which is model-agnostic and far more trustworthy, though it splits credit oddly between correlated features, since shuffling one leaves the other to compensate.

The overriding procedural rule: selection is part of the model, so it must happen inside cross-validation. Selecting features using the whole dataset and then cross-validating the reduced set is leakage, and it produces impressively wrong scores. The demonstration below shows it selecting pure noise and reporting high accuracy.

import numpy as np
import pandas as pd
from sklearn.datasets import make_classification
from sklearn.feature_selection import (SelectKBest, f_classif,
                                       mutual_info_classif, RFECV)
from sklearn.linear_model import LassoCV, LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.preprocessing import StandardScaler

Xa, ya = make_classification(n_samples=1200, n_features=30, n_informative=5,
                            n_redundant=5, n_repeated=0, shuffle=False,
                            random_state=0)
X = pd.DataFrame(Xa, columns=[f"f{i:02d}" for i in range(30)])
y = pd.Series(ya)
# By construction f00-f04 are informative, f05-f09 are redundant copies,
# and f10-f29 are noise.

# --- FILTER: fast, univariate, blind to interactions. ---
f_scores = pd.Series(f_classif(X, y)[0], index=X.columns)
mi = pd.Series(mutual_info_classif(X, y, random_state=0), index=X.columns)
print("top 8 by ANOVA F:", f_scores.nlargest(8).index.tolist())
print("top 8 by mutual information:", mi.nlargest(8).index.tolist())

# THE BLIND SPOT: XOR. Each feature alone is uninformative.
rng = np.random.default_rng(0)
a, b = rng.integers(0, 2, 2000), rng.integers(0, 2, 2000)
xor = pd.DataFrame({"a": a, "b": b})
yx = pd.Series(a ^ b)
print("\nXOR: univariate mutual information per feature:",
      np.round(mutual_info_classif(xor, yx, random_state=0), 4).tolist())
print("accuracy of a model using both:",
      round(cross_val_score(RandomForestClassifier(random_state=0),
                            xor, yx, cv=5).mean(), 3))
# Both scores are ~0 and the pair predicts perfectly. Any filter method
# would discard both features.

# --- EMBEDDED: L1 drives coefficients to exactly zero. ---
lasso = Pipeline([("s", StandardScaler()),
                  ("m", LassoCV(cv=5, random_state=0, max_iter=5000))]).fit(X, y)
coefs = pd.Series(lasso.named_steps["m"].coef_, index=X.columns)
print(f"\nlasso kept {(coefs != 0).sum()} of {len(coefs)} features:"
      f" {coefs[coefs != 0].index.tolist()[:8]}")

# --- IMPURITY IMPORTANCE IS BIASED toward many-valued features. ---
Xb = X.copy()
Xb["noise_continuous"] = rng.normal(size=len(Xb))    # pure noise, many values
Xb["noise_binary"] = rng.integers(0, 2, len(Xb))     # pure noise, two values
rf = RandomForestClassifier(n_estimators=250, random_state=0).fit(Xb, y)
imp = pd.Series(rf.feature_importances_,
                index=Xb.columns).sort_values(ascending=False)
rank = {c: i + 1 for i, c in enumerate(imp.index)}
print(f"\nimpurity rank of noise_continuous: {rank['noise_continuous']}/{len(imp)}")
print(f"impurity rank of noise_binary:     {rank['noise_binary']}/{len(imp)}")
print("the continuous noise column outranks the binary one purely because")
print("it offers more candidate split points.")

# --- PERMUTATION IMPORTANCE: measured on held-out data. ---
Xtr, Xte, ytr, yte = train_test_split(Xb, y, test_size=0.3, random_state=0)
rf2 = RandomForestClassifier(n_estimators=200, random_state=0).fit(Xtr, ytr)
perm = permutation_importance(rf2, Xte, yte, n_repeats=10, random_state=0)
pi = pd.Series(perm.importances_mean, index=Xb.columns).sort_values(ascending=False)
print("\ntop 6 by permutation importance:", pi.head(6).round(4).to_dict())
print(f"noise_continuous permutation importance: "
      f"{pi['noise_continuous']:+.4f}  <- at or below zero, correctly")

# --- WRAPPER: RFE with cross-validation picks the subset size. ---
rfecv = RFECV(LogisticRegression(max_iter=2000), step=2, cv=5,
              scoring="accuracy", min_features_to_select=2)
rfecv.fit(StandardScaler().fit_transform(X), y)
print(f"\nRFECV chose {rfecv.n_features_} features")

# --- THE LEAKAGE, DEMONSTRATED ON PURE NOISE ---
noise = pd.DataFrame(rng.normal(size=(200, 2000)))    # NOTHING is predictive
rand_y = pd.Series(rng.integers(0, 2, 200))

# WRONG: select on all the data, then cross-validate the survivors.
best = SelectKBest(f_classif, k=10).fit(noise, rand_y).get_support(indices=True)
leaky = cross_val_score(LogisticRegression(max_iter=1000),
                        noise.iloc[:, best], rand_y, cv=5).mean()

# RIGHT: selection inside the pipeline, re-run on every training fold.
honest = cross_val_score(
    Pipeline([("sel", SelectKBest(f_classif, k=10)),
              ("m", LogisticRegression(max_iter=1000))]),
    noise, rand_y, cv=5).mean()

print(f"\nPURE NOISE, true accuracy must be ~0.50")
print(f"  select-then-validate: {leaky:.3f}   <- leakage")
print(f"  selection in pipeline: {honest:.3f}  <- honest")
# With 2000 noise columns, 10 of them correlate with the label by chance in
# THIS sample. Selecting them using the labels you then score against is
# how a pipeline reports skill it does not have.

# --- REDUNDANCY: drop one of a near-duplicate pair. ---
corr = X.corr().abs()
upper = corr.where(np.triu(np.ones(corr.shape), k=1).astype(bool))
drop = [c for c in upper.columns if (upper[c] > 0.95).any()]
print(f"\nfeatures above 0.95 correlation with an earlier one: {drop}")

  THE THREE FAMILIES

  family     how                    cost    interactions   ships with model
  ---------------------------------------------------------------------------
  filter     score each feature     O(p)    NO             no
  wrapper    refit on subsets       O(p^2)  yes            no
  embedded   penalty or importance  O(1)    yes            YES

  WHY FILTERS MISS XOR

    a  b  y=a^b       marginal association of a with y: zero
    0  0  0           marginal association of b with y: zero
    0  1  1           joint: perfectly determined
    1  0  1
    1  1  0

  THE SELECTION LEAKAGE LOOP

    WRONG                                RIGHT
    select on all data                   for each fold:
        |                                    select on THIS fold's train
    split into folds                         fit on THIS fold's train
        |                                    score on THIS fold's test
    cross-validate                       average

    The wrong version lets the test labels choose the features, and
    reports ~0.70 accuracy on data with no signal whatsoever.

  IMPORTANCE: WHICH TO TRUST

    impurity      fast, biased to high-cardinality, computed on TRAIN
    permutation   slower, honest, computed on HELD-OUT data
    SHAP          per-row attributions, consistent, slowest

Key Takeaways

• Extra features cost accuracy through variance, instability, and the curse of dimensionality, and they cost operationally through drift and monitoring.
• Filter methods are fast and evaluate features one at a time, so they discard features that only matter jointly.
• Embedded methods such as L1 penalties and tree importances are the best default because selection is consistent with the model being fitted.
• Impurity importance favours high-cardinality features; prefer permutation importance measured on held-out data.
• Selection is part of the model and must run inside cross-validation, or you can report high accuracy on pure noise.

🧪 Practice

1. Build an XOR dataset, show that mutual information for each feature is near zero, and that a model using both predicts almost perfectly.
2. Add a pure-noise continuous column to a dataset and compare its impurity importance rank with its permutation importance.
3. Reproduce the noise experiment: select the ten best of two thousand random columns before and inside cross-validation, and report both accuracies.
4. Interview: A stakeholder asks you to cut a 400-feature model to 20 without losing accuracy. What is your plan and what do you measure? (Hint: which family of methods gives you both a ranking and a stopping rule in one fit?)

Handling High-Cardinality Features

High cardinality means a categorical variable with very many distinct values: postcodes, product SKUs, user IDs, IP addresses, merchant codes, URLs. They are common, they often carry real signal, and every standard encoding handles them badly. One-hot gives you an unusable number of columns. Ordinal invents an order over ten thousand levels. Target encoding works and leaks if you are careless.

The first question is diagnostic and frequently decides everything: is the identifier itself predictive, or is it a proxy for attributes you could use directly? A postcode predicts default risk because it proxies income, density, and local economy. If you can obtain those attributes, use them: they generalise to postcodes you have never seen, they are stable, and they are explainable. The identifier only earns its place when the entity-specific effect is real and irreducible.

The second question is about the tail. Cardinality distributions are almost always extremely skewed: a few values cover most rows and thousands appear once or twice. Those singletons cannot support an estimate of anything. The standard moves:

• Collapse the tail. Keep values covering, say, 95% of rows or appearing at least k times, and map the rest to a single "rare" bucket. One threshold, fitted on training data, and most of the problem goes away.
• Group semantically. A postcode's first two characters are a region. A URL has a domain. A SKU has a category prefix. Hierarchical structure is usually present in the identifier itself and is far more stable than the leaf.
• Encode with a statistic. Frequency and smoothed out-of-fold target encoding, as covered above.
• Hash. Fixed width, no vocabulary, accepts collisions. The right choice when new values arrive continuously and you cannot refit a vocabulary.
• Learn an embedding. A dense vector per level, learned jointly with the model. This is what neural networks do with entity embeddings and it works well given enough rows per level. It is also how recommender systems handle millions of item IDs.
• Use a model with native support. LightGBM and CatBoost handle categorical columns internally, with CatBoost's ordered target statistics designed specifically for this problem.

The operational trap to plan for: the category set changes. Tomorrow brings a new SKU. Your encoder must have a defined behaviour for it, the fallback value must be sensible, and you should monitor the rate of unseen values, because a spike in that rate is an early warning of upstream change that no accuracy metric will show you until much later.

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import OneHotEncoder
from sklearn.feature_extraction import FeatureHasher
from sklearn.ensemble import HistGradientBoostingClassifier

rng = np.random.default_rng(13)
n = 6000
# A realistic skewed identifier: 3000 levels, most of them nearly unused.
weights = rng.dirichlet(np.full(3000, 0.25))
sku = rng.choice([f"SKU-{c:02d}-{i:04d}" for c in range(8)
                  for i in range(375)], n, p=weights)
df = pd.DataFrame({"sku": sku})
df["category"] = df["sku"].str.slice(4, 6)          # the hierarchy is IN the id
y = pd.Series((rng.random(n) < df["category"].map(
    dict(zip([f"{i:02d}" for i in range(8)], rng.beta(2, 6, 8))))).astype(int))

vc = df["sku"].value_counts()
print(f"{len(vc)} distinct SKUs in {n} rows")
print(f"rows covered by the top 100 SKUs: {vc.head(100).sum() / n:.1%}")
print(f"SKUs appearing once: {(vc == 1).sum()} "
      f"({(vc == 1).sum() / len(vc):.0%} of levels, "
      f"{vc[vc == 1].sum() / n:.1%} of rows)")

# --- THE ONE-HOT COST ---
print(f"\none-hot would produce {len(vc)} columns for {n} rows:"
      f" {len(vc) / n:.2f} columns per row")

# --- MOVE 1: COLLAPSE THE TAIL. One fitted threshold. ---
def fit_tail_collapse(s: pd.Series, min_count: int = 20):
    keep = set(s.value_counts()[lambda v: v >= min_count].index)
    return keep

def apply_tail_collapse(s: pd.Series, keep: set, other: str = "__other__"):
    return s.where(s.isin(keep), other)

tr_idx, te_idx = train_test_split(df.index, test_size=0.3, random_state=0)
keep = fit_tail_collapse(df.loc[tr_idx, "sku"], min_count=20)
df["sku_collapsed"] = apply_tail_collapse(df["sku"], keep)
other_share = (df["sku_collapsed"] == "__other__").mean()
print(f"\ncollapsed to {df['sku_collapsed'].nunique()} levels;"
      f" '__other__' covers {other_share:.1%} of rows")
# `keep` is a fitted parameter. Recomputing it on test data leaks and also
# makes train and test disagree about which levels exist.

# --- MOVE 2: SEMANTIC GROUPING. The identifier contains a hierarchy. ---
print(f"\nhierarchy from the id: {df['category'].nunique()} categories")
grouped_rate = y.groupby(df["category"]).agg(["mean", "size"])
print(grouped_rate.round(3).to_string())
print("this generalises to SKUs never seen, because the category is known")

# --- MOVE 3: STATISTICS. Frequency plus out-of-fold target encoding. ---
df["sku_freq"] = df["sku"].map(df["sku"].value_counts(normalize=True))

from sklearn.model_selection import KFold
def oof_te(keys, target, m=30.0, splits=5):
    out = pd.Series(np.nan, index=keys.index, dtype=float)
    for fit_i, app_i in KFold(splits, shuffle=True, random_state=0).split(keys):
        t, k = target.iloc[fit_i], keys.iloc[fit_i]
        prior = t.mean()
        st = t.groupby(k).agg(["sum", "count"])
        mp = (st["sum"] + m * prior) / (st["count"] + m)
        out.iloc[app_i] = keys.iloc[app_i].map(mp).fillna(prior).values
    return out

df["sku_te"] = oof_te(df["sku"], y)
print(f"\nout-of-fold target encoding: 1 column,"
      f" correlation {df['sku_te'].corr(y):+.3f}")

# --- MOVE 4: HASHING. Fixed width, no vocabulary, tolerates new values. ---
for width in (64, 256, 1024):
    H = FeatureHasher(n_features=width, input_type="string").transform(
        [[v] for v in df["sku"]])
    distinct = len(np.unique(H.toarray(), axis=0))
    print(f"hash to {width:5d} cols -> {distinct:5d} distinct patterns"
          f" of {len(vc)} levels ({1 - distinct/len(vc):.0%} collided)")
# Collisions merge unrelated levels. Wider is safer and costs more.

# --- MOVE 5: NATIVE CATEGORICAL SUPPORT ---
native = df[["sku_collapsed", "category"]].astype("category")
model = HistGradientBoostingClassifier(categorical_features="from_dtype",
                                       random_state=0, max_iter=120)
print("\nnative categorical handling, CV accuracy:",
      round(cross_val_score(model, native, y, cv=5).mean(), 3))

# --- COMPARE THE REPRESENTATIONS HONESTLY ---
import warnings; warnings.filterwarnings("ignore")
options = {
    "category only (8 levels)": pd.get_dummies(df["category"]),
    "frequency only": df[["sku_freq"]],
    "collapsed one-hot": pd.get_dummies(df["sku_collapsed"]),
    "category + frequency": pd.concat(
        [pd.get_dummies(df["category"]), df[["sku_freq"]]], axis=1),
}
print("\nCV accuracy and width:")
for name, rep in options.items():
    s = cross_val_score(HistGradientBoostingClassifier(random_state=0, max_iter=80),
                        rep, y, cv=5).mean()
    print(f"  {name:26s} {rep.shape[1]:5d} cols   {s:.3f}")
# The 8-level hierarchy does most of the work. That is the common outcome:
# the identifier was a proxy, and the attribute it proxied is better.

# --- OPERATIONAL: monitor the unseen rate. ---
train_levels = set(df.loc[tr_idx, "sku"])
unseen_rate = (~df.loc[te_idx, "sku"].isin(train_levels)).mean()
print(f"\nunseen SKU rate on held-out data: {unseen_rate:.1%}")
print("track this in production; a jump means the upstream catalogue changed")
print("and it will move before any accuracy metric does.")

  THE SKEW YOU WILL ALWAYS FIND

  rows
   |#
   |#
   |##
   |###
   |#####___________________________________________
   +--------------------------------------------------> levels, sorted
    top 100 levels: ~60% of rows     |  thousands of singletons: ~10% of rows
    estimable                        |  not estimable, ever

  STRATEGY BY SITUATION

  situation                              approach
  ---------------------------------------------------------------------------
  the id proxies known attributes        use the attributes, drop the id
  the id has structure (prefix, domain)  extract the hierarchy
  moderate cardinality, enough rows      collapse tail + one-hot
  high cardinality, enough rows/level    out-of-fold target encoding
  new levels arrive constantly           hashing
  very many levels, deep model, big n    learned embedding
  gradient boosting available            native categorical support

  WIDTH COMPARISON, 10,000 LEVELS

    one-hot       10,000 columns
    hashing(512)      512 columns, with collisions
    binary             14 columns, uninterpretable
    frequency           1 column
    target              1 column, needs out-of-fold discipline
    embedding(16)      16 columns, learned

  ALWAYS: a defined fallback for unseen levels, and a monitor on their rate.

Key Takeaways

• Ask first whether the identifier is predictive in itself or a proxy for attributes you could use directly, because the attributes generalise to unseen levels.
• Cardinality distributions are heavily skewed, so collapsing the tail at a fitted count threshold removes most of the problem.
• Identifiers usually contain a hierarchy such as a region, domain, or category prefix, and the coarser level is far more stable.
• Hashing gives fixed width with no vocabulary and is the right answer when new levels arrive continuously.
• Define the fallback for unseen levels explicitly and monitor the unseen rate, because it moves before accuracy does.

🧪 Practice

1. For a high-cardinality column, report how many rows the top hundred levels cover and what share of levels appear only once.
2. Collapse the tail at a threshold fitted on training data, then apply it to a test split and confirm no new level appears.
3. Compare cross-validated accuracy for the raw identifier, an extracted hierarchy, and a frequency encoding, with the column count for each.
4. Interview: A fraud model keys on merchant ID and performs well offline. What is your main concern about deploying it? (Hint: what fraction of tomorrow's merchants existed in your training window?)

Leakage Prevention

Data leakage is when information reaches the model that will not be available at prediction time. It is the most damaging error in applied machine learning, because its symptom is success. The model scores wonderfully, everyone is pleased, it ships, and production performance is a fraction of what was promised. Nothing about a leaked model looks wrong until it is live.

Leakage takes three forms, and it is worth being able to name them.

Target leakage happens when a feature contains the answer, often indirectly. A hospital readmission model that uses "discharge medication count" is fine; one that uses "days in the follow-up program" is not, because enrolment happens after readmission. The classic tell is a feature that is implausibly predictive on its own. If one column gets you 0.99 AUC, you have not found a great feature, you have found the label.

Train-test contamination happens when the preprocessing that produced the features saw the test rows. Every fitted transform is a channel: imputation medians, scaler means, bin edges, target encodings, PCA components, vectorizer vocabularies, feature selection. Fitting any of them on the full dataset before splitting moves information across the boundary. This is by far the most common form, and the fix is entirely mechanical: put every transform in a pipeline and let cross-validation refit it per fold.

Temporal leakage happens when training data includes information from after the prediction moment. A random split of time-ordered data lets the model learn from a future it will not have. Aggregates computed over the full history do the same thing more subtly. The fix is a time-based split plus a stated as-of timestamp for every feature.

There is a fourth, easy to miss: group leakage. If one patient contributes several visits, or one user several sessions, a random split puts the same entity in both training and test. The model recognises the entity rather than learning the pattern, and the score measures memorisation. Split by group.

A practical defence checklist, in the order that catches the most:

1. Split first. Do everything else after.
2. Wrap every fitted transform in a pipeline.
3. For each feature, write the timestamp at which its value is known. If you cannot, that is your answer.
4. Split by time for temporal problems and by group for repeated entities.
5. Be suspicious of any single feature with implausible standalone power.
6. Check feature-target correlations and investigate anything extreme.
7. Compare a model's validation score against a plausible ceiling. If a fraud model is at 0.999 AUC, something is wrong.
8. Keep a final holdout used once, at the end, to confirm the story.

The reason to internalise this list rather than trusting your instincts is that leakage feels like competence while it is happening.

import numpy as np
import pandas as pd
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import (cross_val_score, train_test_split,
                                     GroupKFold, TimeSeriesSplit)
from sklearn.metrics import roc_auc_score

rng = np.random.default_rng(21)
n = 1500

# --- FORM 1: TARGET LEAKAGE. A column that exists because of the label. ---
churned = rng.random(n) < 0.25
data = pd.DataFrame({
    "tenure": rng.integers(1, 60, n),
    "monthly_spend": rng.gamma(3, 20, n).round(2),
    # This field is populated by the retention team AFTER churn is known.
    "retention_call_made": np.where(churned, rng.random(n) < 0.9,
                                    rng.random(n) < 0.02).astype(int),
})
y = pd.Series(churned.astype(int))

honest_cols = ["tenure", "monthly_spend"]
print("AUC with an honest feature set:",
      round(cross_val_score(LogisticRegression(max_iter=1000),
                            data[honest_cols], y, cv=5, scoring="roc_auc").mean(), 3))
print("AUC including the post-outcome field:",
      round(cross_val_score(LogisticRegression(max_iter=1000),
                            data, y, cv=5, scoring="roc_auc").mean(), 3))

# THE DETECTOR: score each feature alone. An implausible single-feature AUC
# is the signature of a label in disguise.
print("\nsingle-feature AUC:")
for col in data.columns:
    print(f"  {col:22s} {roc_auc_score(y, data[col]):.3f}")
print("  anything above ~0.90 alone deserves an explanation, not celebration")

# --- FORM 2: TRAIN-TEST CONTAMINATION via a fitted transform. ---
noise = pd.DataFrame(rng.normal(size=(300, 3000)))
rand_y = pd.Series(rng.integers(0, 2, 300))

sel = SelectKBest(f_classif, k=15).fit(noise, rand_y)       # sees all labels
leaky = cross_val_score(LogisticRegression(max_iter=1000),
                        noise.iloc[:, sel.get_support(indices=True)],
                        rand_y, cv=5, scoring="roc_auc").mean()
honest = cross_val_score(
    Pipeline([("sel", SelectKBest(f_classif, k=15)),
              ("m", LogisticRegression(max_iter=1000))]),
    noise, rand_y, cv=5, scoring="roc_auc").mean()
print(f"\nPURE NOISE (true AUC = 0.500)")
print(f"  selection before CV: {leaky:.3f}")
print(f"  selection inside CV: {honest:.3f}")

# THE FIX IS STRUCTURAL: one pipeline holding every fitted step.
num = ["tenure", "monthly_spend"]
safe_pipeline = Pipeline([
    ("prep", ColumnTransformer([
        ("num", Pipeline([("imp", SimpleImputer(strategy="median")),
                          ("sc", StandardScaler())]), num),
    ])),
    ("select", SelectKBest(f_classif, k=2)),
    ("model", LogisticRegression(max_iter=1000)),
])
print("\npipelined AUC on honest features:",
      round(cross_val_score(safe_pipeline, data[num], y, cv=5,
                            scoring="roc_auc").mean(), 3))
print("every fitted step above is refitted per fold, automatically")

# --- FORM 3: TEMPORAL LEAKAGE. A random split on ordered data. ---
dates = pd.date_range("2023-01-01", periods=n, freq="6h")
drift = np.linspace(0, 2.5, n)                     # the relationship changes
ts = pd.DataFrame({"x": rng.normal(size=n) + drift, "date": dates})
ts_y = pd.Series((rng.random(n) < 1 / (1 + np.exp(-(ts["x"] - drift)))).astype(int))

random_cv = cross_val_score(LogisticRegression(max_iter=1000), ts[["x"]],
                            ts_y, cv=5, scoring="roc_auc").mean()
time_cv = cross_val_score(LogisticRegression(max_iter=1000), ts[["x"]], ts_y,
                          cv=TimeSeriesSplit(5), scoring="roc_auc").mean()
print(f"\nrandom 5-fold CV     {random_cv:.3f}   <- optimistic")
print(f"time-series split CV {time_cv:.3f}   <- what production looks like")
print("every random fold contains rows from after its own test rows")

# --- FORM 4: GROUP LEAKAGE. The same entity on both sides. ---
groups = rng.integers(0, 120, n)                  # 120 users, many rows each
user_effect = rng.normal(0, 2, 120)[groups]       # a strong per-user effect
gX = pd.DataFrame({"x": rng.normal(size=n) + user_effect})
gy = pd.Series((rng.random(n) < 1 / (1 + np.exp(-user_effect))).astype(int))

naive = cross_val_score(LogisticRegression(max_iter=1000), gX, gy, cv=5,
                        scoring="roc_auc").mean()
grouped = cross_val_score(LogisticRegression(max_iter=1000), gX, gy,
                          cv=GroupKFold(5), groups=groups,
                          scoring="roc_auc").mean()
print(f"\nrandom CV  {naive:.3f}   grouped CV {grouped:.3f}")
print("the gap is memorised users, not learned behaviour")

# --- THE AS-OF AUDIT: write the availability time of every feature. ---
FEATURE_TIMING = pd.DataFrame([
    ("tenure",              "signup date",        "before prediction", "safe"),
    ("monthly_spend",       "billing, prior mo.",  "before prediction", "safe"),
    ("retention_call_made", "retention workflow",  "AFTER outcome",     "LEAK"),
    ("support_tickets_30d", "ticket system",       "before prediction", "safe"),
    ("cancellation_reason", "cancellation form",   "AFTER outcome",     "LEAK"),
], columns=["feature", "source", "available", "verdict"])
print("\n", FEATURE_TIMING.to_string(index=False))
# Producing this table is the single highest-value hour in any project.
# Every row you cannot fill in confidently is a row to investigate.

def audit_leakage(X: pd.DataFrame, y: pd.Series, auc_alarm: float = 0.90):
    """Flag any single feature whose standalone power is implausible."""
    flags = []
    for col in X.select_dtypes("number").columns:
        try:
            auc = roc_auc_score(y, X[col])
        except ValueError:
            continue
        auc = max(auc, 1 - auc)                    # direction does not matter
        if auc >= auc_alarm:
            flags.append((col, round(auc, 3)))
    return flags

print("\nautomated flag:", audit_leakage(data, y))

  THE FOUR FORMS

  form              information crossing            fix
  ---------------------------------------------------------------------------
  target            a feature caused by the label   audit feature provenance
  contamination     a transform fitted on all data  Pipeline + CV
  temporal          rows from after the prediction  TimeSeriesSplit, as-of
  group             same entity on both sides       GroupKFold

  EVERY FITTED TRANSFORM IS A CHANNEL

    imputation median      scaler mean/sd        bin edges
    target encoding map    PCA components        vectorizer vocabulary
    feature selection      outlier clip bounds   resampling (SMOTE)

    If it has a .fit(), it must be inside the pipeline.

  THE CORRECT ORDER

    raw data
       |
    SPLIT  <----------------- first, always
       |
       +-- train ---> fit transforms ---> fit model
       |                   |                  |
       +-- test  ----> apply the SAME --> score
                          fitted objects

  SMELL TEST

    one feature alone gives AUC > 0.95        -> it is the label
    validation >> test, gap grows with        -> leaky target encoding
      cardinality
    great offline, poor live, no code change  -> temporal or group leakage
    accuracy above the known human ceiling    -> stop and investigate

Key Takeaways

• Leakage's symptom is an excellent score, which is why it survives review and fails in production.
• Any object with a fit step is a leakage channel; putting every one inside a pipeline closes the whole class mechanically.
• Split by time for temporal problems and by group for repeated entities, because a random split breaks both.
• Write down when each feature's value becomes available relative to the prediction moment; a feature you cannot place is a suspect.
• Treat an implausibly strong single feature or a score above the known ceiling as evidence of leakage rather than of success.

🧪 Practice

1. Construct a dataset with a feature populated after the outcome, then show the cross-validated AUC with and without it.
2. Reproduce the noise experiment for feature selection and report both scores, then explain to a non-specialist why the higher one is wrong.
3. Build panel data with a strong per-entity effect and compare random k-fold against grouped k-fold scores.
4. Write the feature-availability audit table for a model you have built or read about, marking each feature safe or leaking.
5. Interview: A colleague's model reaches 0.99 AUC on fraud detection. Which three questions do you ask first? (Hint: when was each feature written, how was the data split, and what did the preprocessing see?)

<a id="64-handling-imbalanced-and-skewed-data"></a>

6.4 Handling Imbalanced and Skewed Data

Real target variables are rarely balanced and real numeric variables are rarely symmetric; this section covers how to diagnose both conditions, the techniques that genuinely help, and the much larger set of techniques that only appear to.

Class Imbalance Diagnosis

Fraud is 0.2% of transactions. Rare disease is 1% of patients. Click-through is 3% of impressions. Imbalance is the normal state of any interesting binary problem, and the first thing to understand is that imbalance itself is not automatically a problem. Diagnosing whether it is a problem, and what kind, comes before reaching for any technique.

The headline failure is that accuracy stops meaning anything. A classifier that predicts "not fraud" for every transaction is 99.8% accurate and worth nothing. This is easy to state and easy to fall for, because 99.8% reads as excellent to anyone not looking closely.

The metrics that survive imbalance:

Precision is, of the cases you flagged, how many were real. It answers "how much of my review team's time is wasted?" Recall is, of the real cases, how many you caught. It answers "how much fraud got through?" These trade off against each other and you cannot maximise both.

F1 is their harmonic mean, useful as a single number and opaque about which side it came from. Report precision and recall alongside it.

ROC AUC measures ranking quality across all thresholds. Its known weakness under heavy imbalance is that the false positive rate has a huge denominator, so thousands of false positives barely move it. A model can look strong by ROC AUC and be unusable in practice.

PR AUC, the area under the precision-recall curve, uses precision instead and therefore reflects the cost that actually hurts. Under heavy imbalance it is the more informative summary, and note that its baseline is the positive rate rather than 0.5.

The second insight is that the real problem is often not the ratio but the absolute count. A thousand positives in a hundred thousand rows is a 1% imbalance and enough data to learn from. Twenty positives in two thousand rows is the same ratio and hopeless, because twenty examples cannot define a decision boundary in any interesting feature space. Before resampling anything, count the minority class in absolute terms. If it is small, you have a data collection problem, and no technique will substitute.

The third insight is about thresholds. Most classifiers produce a probability, and the 0.5 cutoff is a default, not a decision. Changing it moves you along the precision-recall curve at no training cost. A great deal of what people attempt with resampling is better achieved by picking a threshold from the costs of the two error types, which is a business input, not a modelling one.

Finally, decide whether you need calibration. If a downstream process uses the probability itself, for example to compute expected loss, the probabilities must mean what they say. Resampling deliberately breaks calibration by changing the base rate, which is one of its least-discussed costs.

import numpy as np
import pandas as pd
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split, StratifiedKFold
from sklearn.linear_model import LogisticRegression
from sklearn.dummy import DummyClassifier
from sklearn.metrics import (accuracy_score, precision_score, recall_score,
                             f1_score, roc_auc_score, average_precision_score,
                             confusion_matrix, precision_recall_curve)

X, y = make_classification(n_samples=20000, n_features=20, n_informative=6,
                           weights=[0.995, 0.005], flip_y=0.01, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0,
                                      stratify=y)   # stratify, always
print(f"positives: {y.sum()} of {len(y)} ({y.mean():.3%})")
print(f"test set positives: {yte.sum()}")

# --- ACCURACY IS USELESS HERE ---
dummy = DummyClassifier(strategy="most_frequent").fit(Xtr, ytr)
print(f"\n'always negative' accuracy: {accuracy_score(yte, dummy.predict(Xte)):.4f}")
print(f"'always negative' recall:   {recall_score(yte, dummy.predict(Xte)):.4f}")

model = LogisticRegression(max_iter=2000).fit(Xtr, ytr)
proba = model.predict_proba(Xte)[:, 1]
pred = (proba >= 0.5).astype(int)

cm = confusion_matrix(yte, pred)
print("\nconfusion matrix at threshold 0.5:")
print(pd.DataFrame(cm, index=["actual neg", "actual pos"],
                   columns=["pred neg", "pred pos"]).to_string())

print("\nmetrics at 0.5:")
print(f"  accuracy   {accuracy_score(yte, pred):.4f}   <- looks great, means little")
print(f"  precision  {precision_score(yte, pred, zero_division=0):.4f}")
print(f"  recall     {recall_score(yte, pred, zero_division=0):.4f}")
print(f"  F1         {f1_score(yte, pred, zero_division=0):.4f}")
print(f"  ROC AUC    {roc_auc_score(yte, proba):.4f}   <- threshold-free ranking")
print(f"  PR AUC     {average_precision_score(yte, proba):.4f}"
      f"   <- baseline is {yte.mean():.4f}, not 0.5")

# --- WHY ROC AUC FLATTERS: the FPR denominator is enormous. ---
tn, fp, fn, tp = cm.ravel()
print(f"\n{fp} false positives out of {tn + fp} negatives"
      f" -> FPR {fp / (tn + fp):.4f}")
print("a reviewer sees those false positives; the ROC curve barely notices them")

# --- RATIO VS ABSOLUTE COUNT: the same imbalance, different problems. ---
print("\nsame 1% imbalance, different feasibility:")
for total in (2_000, 20_000, 200_000):
    print(f"  {total:>7,} rows -> {int(total * 0.01):>5,} positives"
          f"   {'too few to learn a boundary' if total * 0.01 < 100 else 'workable'}")

# --- THRESHOLD SELECTION: free movement along the curve. ---
prec, rec, thr = precision_recall_curve(yte, proba)
print("\nthreshold sweep:")
for t in (0.02, 0.05, 0.10, 0.30, 0.50):
    p = (proba >= t).astype(int)
    print(f"  t={t:<5} precision {precision_score(yte, p, zero_division=0):.3f}"
          f"  recall {recall_score(yte, p, zero_division=0):.3f}"
          f"  flagged {p.sum():>5}")
# No retraining happened. The whole curve was already available.

# --- COST-OPTIMAL THRESHOLD: a business input, not a default. ---
COST_FN, COST_FP = 500.0, 10.0       # missed fraud vs a wasted review
grid = np.linspace(0.001, 0.6, 300)
costs = [(t, COST_FN * ((yte == 1) & (proba < t)).sum()
             + COST_FP * ((yte == 0) & (proba >= t)).sum()) for t in grid]
best_t, best_c = min(costs, key=lambda kv: kv[1])
cost_at_half = dict(costs).get(min(grid, key=lambda t: abs(t - 0.5)))
print(f"\ncost-minimising threshold {best_t:.3f} at cost {best_c:,.0f}")
print(f"cost at the default 0.5:  {cost_at_half:,.0f}")
print("the ratio of the two error costs sets the threshold, nothing else")

# --- F-BETA when one error matters more, without full cost numbers. ---
from sklearn.metrics import fbeta_score
for beta, meaning in [(0.5, "precision matters more"), (2.0, "recall matters more")]:
    scores = [(t, fbeta_score(yte, (proba >= t).astype(int), beta=beta,
                              zero_division=0)) for t in grid]
    t_best = max(scores, key=lambda kv: kv[1])[0]
    print(f"F{beta} ({meaning:24s}) peaks at threshold {t_best:.3f}")

# --- STRATIFY EVERY SPLIT, or folds may contain no positives at all. ---
skf = StratifiedKFold(5, shuffle=True, random_state=0)
counts = [ytr[te_i].sum() for _, te_i in skf.split(Xtr, ytr)]
print(f"\npositives per stratified fold: {counts}")
from sklearn.model_selection import KFold
plain = [ytr[te_i].sum() for _, te_i in KFold(5, shuffle=True,
                                              random_state=3).split(Xtr)]
print(f"positives per unstratified fold: {plain}  <- unstable, sometimes zero")

# --- CALIBRATION: does the probability mean what it says? ---
from sklearn.calibration import calibration_curve
frac_pos, mean_pred = calibration_curve(yte, proba, n_bins=5, strategy="quantile")
print("\ncalibration (predicted vs observed):")
for mp, fp_ in zip(mean_pred, frac_pos):
    print(f"  predicted {mp:.4f}  observed {fp_:.4f}")
# If a downstream process multiplies this probability by a loss amount, these
# two columns must agree. Resampling will make them disagree.

  WHY ACCURACY FAILS

  20,000 rows, 100 positives (0.5%)

  always-negative classifier
                pred neg   pred pos
  actual neg    19,900     0            accuracy  = 19,900/20,000 = 99.5%
  actual pos       100     0            recall    = 0
                                        precision = undefined
  A model with zero value scores 99.5%.

  ROC AUC VS PR AUC UNDER IMBALANCE

    FPR = FP / (FP + TN)     denominator ~ 19,900   -> 200 FPs move it by 0.01
    precision = TP / (TP+FP) denominator ~ flagged  -> 200 FPs halve it

    ROC AUC answers "does it rank well?"
    PR AUC answers  "is the flagged list usable?"

  DIAGNOSE BEFORE TREATING

    minority count < ~100          -> data collection problem, not a modelling one
    metric is accuracy             -> change the metric, not the data
    threshold is still 0.5         -> tune the threshold first, it is free
    probabilities feed a decision  -> preserve calibration, avoid resampling
    ranking is all you need        -> imbalance may not need any treatment

Key Takeaways

• Accuracy is meaningless under imbalance; use precision, recall, and the area under the precision-recall curve, whose baseline is the positive rate.
• ROC AUC flatters imbalanced models because the false positive rate has an enormous denominator.
• The absolute count of minority examples matters more than the ratio; too few positives is a data problem no technique solves.
• Threshold tuning moves you along the precision-recall curve for free and often replaces the need for resampling.
• Stratify every split, and remember that resampling breaks calibration if a downstream process uses the probabilities.

🧪 Practice

1. Build a dataset with a 1% positive rate, fit a majority-class classifier, and report its accuracy, precision, and recall.
2. Sweep the decision threshold and tabulate precision, recall, and the number of flagged cases at five thresholds.
3. Given a cost of 500 for a false negative and 10 for a false positive, compute the cost-minimising threshold and compare its cost with the default.
4. Interview: A model reports 0.93 ROC AUC on a 0.3% positive rate and the review team says the alerts are useless. Reconcile the two. (Hint: compute precision at the operating threshold and count the false positives a reviewer actually sees.)

Resampling Techniques

Resampling changes the class ratio in the training data, either by discarding majority examples or by duplicating minority ones. It is the first thing most people try and it is less often the right answer than its popularity suggests, so it is worth being precise about what it does and does not fix.

Random undersampling discards majority rows until the ratio improves. It is fast, it shrinks the dataset so training gets cheaper, and it throws away real data. With a large majority class that redundancy may be genuinely wasteful and undersampling costs little; with a small dataset it discards information you cannot afford to lose.

Random oversampling duplicates minority rows. No information is lost and no information is gained. Exact duplicates give the model repeated identical points, which invites memorising them, and for a model that can fit arbitrarily complex boundaries this is a direct route to overfitting.

Informed undersampling is more interesting than the random kind. Tomek links remove majority points that form a nearest-neighbour pair with a minority point, cleaning the boundary. Edited nearest neighbours removes majority points whose neighbours disagree with them. Cluster centroids replaces the majority class with a smaller set of cluster centres, keeping coverage of the feature space while reducing count. These aim to remove the least useful majority points rather than arbitrary ones.

Combined approaches oversample the minority and then clean the boundary, which is where SMOTE plus Tomek links or SMOTE plus edited nearest neighbours come from.

Three rules govern correct use, and violating any of them invalidates your evaluation.

Resample only the training data. The test set must keep the real class distribution, because that is the distribution production has. Resampling before splitting puts duplicated or synthesised minority points on both sides and produces a meaningless score. This is the most common mistake.

Resample inside the cross-validation loop. The resampler is a fitted step, so each fold must resample only its own training portion. Use a pipeline that understands resampling, such as imblearn's, because a plain scikit-learn pipeline applies its steps to the test fold as well.

Expect probabilities to become miscalibrated. By changing the base rate you have changed what the model's output means. If you need calibrated probabilities, either recalibrate on unresampled data afterwards or prefer class weighting, which leaves the data alone.

A last point that deserves more attention than it gets: before resampling, try class weights and threshold tuning. They are cheaper, they do not distort the data, they preserve calibration better, and in a great many cases they perform as well or better. Resampling is worth reaching for when the minority class is genuinely tiny, when the model has no class-weight option, or when an experiment shows it wins.

import numpy as np
import pandas as pd
from collections import Counter
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split, StratifiedKFold, cross_validate
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import average_precision_score, roc_auc_score, brier_score_loss
from imblearn.under_sampling import (RandomUnderSampler, TomekLinks,
                                     EditedNearestNeighbours, ClusterCentroids)
from imblearn.over_sampling import RandomOverSampler
from imblearn.pipeline import Pipeline as ImbPipeline

X, y = make_classification(n_samples=6000, n_features=15, n_informative=5,
                           weights=[0.97, 0.03], flip_y=0.02, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0,
                                      stratify=y)
print("train distribution:", Counter(ytr), " test distribution:", Counter(yte))

# --- WHAT EACH RESAMPLER DOES TO THE TRAINING SET ---
samplers = {
    "random under": RandomUnderSampler(random_state=0),
    "random over": RandomOverSampler(random_state=0),
    "Tomek links": TomekLinks(),
    "edited NN": EditedNearestNeighbours(),
    "cluster centroids": ClusterCentroids(random_state=0),
}
rows = []
for name, s in samplers.items():
    Xr, yr = s.fit_resample(Xtr, ytr)
    c = Counter(yr)
    rows.append({"sampler": name, "rows": len(yr), "negatives": c[0],
                 "positives": c[1], "ratio": round(c[1] / c[0], 3)})
print("\n", pd.DataFrame(rows).to_string(index=False))
# Note random under discards thousands of real rows; random over adds no new
# information; Tomek and ENN barely change the ratio and clean the boundary.

# Random oversampling creates exact duplicates, which is its whole problem.
Xo, yo = RandomOverSampler(random_state=0).fit_resample(Xtr, ytr)
minority = Xo[yo == 1]
uniq = len(np.unique(minority, axis=0))
print(f"\noversampled minority: {len(minority)} rows, {uniq} distinct"
      f" -> each real point repeated ~{len(minority)/uniq:.1f} times")

# --- THE FATAL MISTAKE: resampling before the split. ---
Xall, yall = RandomOverSampler(random_state=0).fit_resample(X, y)
Xa_tr, Xa_te, ya_tr, ya_te = train_test_split(Xall, yall, test_size=0.3,
                                              random_state=0, stratify=yall)
wrong = RandomForestClassifier(n_estimators=120, random_state=0).fit(Xa_tr, ya_tr)
print(f"\nresample-then-split PR AUC: "
      f"{average_precision_score(ya_te, wrong.predict_proba(Xa_te)[:, 1]):.3f}"
      "   <- meaningless")

right = ImbPipeline([("res", RandomOverSampler(random_state=0)),
                     ("m", RandomForestClassifier(n_estimators=120, random_state=0))
                     ]).fit(Xtr, ytr)
print(f"split-then-resample PR AUC: "
      f"{average_precision_score(yte, right.predict_proba(Xte)[:, 1]):.3f}"
      "   <- honest")
print("the first number is inflated because duplicated minority rows appear")
print("in BOTH the training and the test halves.")

# --- RESAMPLE INSIDE THE CV LOOP, never outside it. ---
cv = StratifiedKFold(5, shuffle=True, random_state=0)
configs = {
    "baseline (nothing)": ImbPipeline([
        ("m", LogisticRegression(max_iter=2000))]),
    "class_weight balanced": ImbPipeline([
        ("m", LogisticRegression(max_iter=2000, class_weight="balanced"))]),
    "random undersample": ImbPipeline([
        ("res", RandomUnderSampler(random_state=0)),
        ("m", LogisticRegression(max_iter=2000))]),
    "random oversample": ImbPipeline([
        ("res", RandomOverSampler(random_state=0)),
        ("m", LogisticRegression(max_iter=2000))]),
    "Tomek cleaning": ImbPipeline([
        ("res", TomekLinks()),
        ("m", LogisticRegression(max_iter=2000))]),
}
print("\ncross-validated on the TRAINING set (test distribution untouched):")
for name, pipe in configs.items():
    res = cross_validate(pipe, Xtr, ytr, cv=cv,
                         scoring=["average_precision", "roc_auc"])
    print(f"  {name:24s} PR AUC {res['test_average_precision'].mean():.3f}"
          f"   ROC AUC {res['test_roc_auc'].mean():.3f}")
# The baseline is often competitive. Measure before assuming resampling helps.

# --- CALIBRATION COST: resampling changes what the output means. ---
plain = LogisticRegression(max_iter=2000).fit(Xtr, ytr)
over = ImbPipeline([("res", RandomOverSampler(random_state=0)),
                    ("m", LogisticRegression(max_iter=2000))]).fit(Xtr, ytr)
print(f"\ntrue positive rate in test: {yte.mean():.4f}")
for name, m in [("no resampling", plain), ("oversampled", over)]:
    p = m.predict_proba(Xte)[:, 1]
    print(f"  {name:15s} mean predicted {p.mean():.4f}"
          f"   Brier {brier_score_loss(yte, p):.5f}")
# The oversampled model's mean prediction tracks the RESAMPLED base rate,
# not the real one. Any expected-value calculation using it is wrong.

# Recalibrating on unresampled data repairs it, at the cost of a holdout.
from sklearn.calibration import CalibratedClassifierCV
recal = CalibratedClassifierCV(over, method="isotonic", cv=3).fit(Xtr, ytr)
pr = recal.predict_proba(Xte)[:, 1]
print(f"  recalibrated    mean predicted {pr.mean():.4f}"
      f"   Brier {brier_score_loss(yte, pr):.5f}")

# --- UNDERSAMPLING WITH AN ENSEMBLE: keep the data, not the imbalance. ---
from imblearn.ensemble import BalancedBaggingClassifier
bag = BalancedBaggingClassifier(
    estimator=RandomForestClassifier(n_estimators=30, random_state=0),
    n_estimators=8, sampling_strategy="auto", random_state=0).fit(Xtr, ytr)
print(f"\nbalanced bagging PR AUC: "
      f"{average_precision_score(yte, bag.predict_proba(Xte)[:, 1]):.3f}")
print("each base learner sees a different balanced subsample, so across the")
print("ensemble no majority row is permanently discarded.")

  WHAT RESAMPLING ACTUALLY DOES

  original          undersampled         oversampled
  N N N N N N N N   N N N               N N N N N N N N
  N N N N N N N N   N N                 N N N N N N N N
  P                 P                   P P P P P P P P     <- P duplicated
                                        P P P P P P P P

  information lost: none | majority rows | none
  information gained:      none          | NONE
  overfitting risk:        lower (less data) | HIGHER (repeated points)
  calibration:             preserved      broken (base rate changed)

  THE ORDER THAT MATTERS

    WRONG                              RIGHT
    resample                           split
       |                                  |
    split  <- duplicates land on both   resample TRAIN only
       |       sides of the split          |
    score (inflated)                   score on untouched test

  TRY IN THIS ORDER

    1. change the metric        (free, and often the whole problem)
    2. tune the threshold       (free, no retraining)
    3. class weights            (one parameter, keeps all data)
    4. resample                 (distorts data, breaks calibration)
    5. collect more positives   (best, slowest)

Key Takeaways

• Undersampling discards real data and oversampling adds no information, so neither creates signal that was not there.
• Resample only the training portion, inside the cross-validation loop, using a pipeline that leaves test folds untouched.
• Resampling before splitting puts duplicated minority rows on both sides and produces a meaningless score.
• Informed undersampling such as Tomek links or edited nearest neighbours cleans the boundary instead of discarding rows arbitrarily.
• Try metric choice, threshold tuning, and class weights before resampling, and expect to recalibrate if you resample anyway.

🧪 Practice

1. Resample an imbalanced dataset before and after splitting, and report the gap between the two held-out precision-recall scores.
2. Compare cross-validated PR AUC for a baseline, class weighting, random undersampling, and random oversampling on the same data.
3. Show that oversampling shifts the mean predicted probability away from the true positive rate, then repair it with calibration.
4. Interview: Your teammate reports that oversampling lifted accuracy from 0.97 to 0.99. What do you check? (Hint: where in the workflow did the resampling happen, and which metric would have revealed nothing?)

SMOTE and Synthetic Data

Random oversampling duplicates existing minority points, which gives a model the same information repeated. SMOTE, the Synthetic Minority Over-sampling Technique, tries to do better by manufacturing new points instead. The idea is simple and geometric: take a minority point, pick one of its k nearest minority neighbours, and place a new point at a random position on the line segment between them.

The intuition is that the space between two examples of the same class is probably also that class, so filling it in gives the model a denser region to learn a boundary around rather than a handful of isolated points to memorise.

That intuition is also exactly where SMOTE fails, and the failure modes are worth knowing precisely because SMOTE is applied reflexively.

It interpolates across gaps. If the minority class forms two separate clusters, the segment between a point in one cluster and a point in the other passes through majority territory, and SMOTE will happily put synthetic minority points there. The model then learns a boundary around a region that contains no real minority examples at all.

It amplifies noise. A mislabelled minority point deep inside the majority region becomes the seed for a line of synthetic points reaching further into that region. One bad label becomes many.

It struggles with categorical features. The midpoint of "red" and "blue" is not a colour. SMOTE-NC handles this by interpolating numeric features and taking the most common category among neighbours, which is a workable compromise rather than a solution.

It degrades in high dimensions. Nearest neighbours become unreliable as dimensionality grows, so "the line between a point and its neighbour" becomes a less meaningful place to sample from. In very sparse high-dimensional spaces, such as text features, SMOTE is usually not appropriate.

The variants address specific weaknesses. Borderline-SMOTE only synthesises from minority points near the boundary, where the extra density is most useful. ADASYN generates more points for minority examples in harder neighbourhoods. SMOTE-Tomek and SMOTE-ENN synthesise and then clean up the resulting boundary. SVM-SMOTE uses support vectors to choose seeds.

The honest empirical picture is that SMOTE helps sometimes, does nothing sometimes, and hurts sometimes, and which of those happens is not reliably predictable from the dataset's description. Several careful studies have found that on tabular problems with strong modern classifiers, SMOTE rarely beats simple class weighting. So the guidance is: treat it as one candidate among several, evaluate it honestly on a test set with the real class distribution, and do not assume it because the class ratio is skewed.

Beyond SMOTE there is a broader category of synthetic data generation: generative adversarial networks, variational autoencoders, and copula or Bayesian-network models that learn the joint distribution and sample from it. These matter more for privacy-preserving data sharing than for imbalance, and they carry the same fundamental caveat, in a stronger form. A generator can only reproduce structure it learned from real data. It cannot add information, and it can leak the training examples it memorised.

import numpy as np
import pandas as pd
from collections import Counter
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split, StratifiedKFold, cross_validate
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import NearestNeighbors
from sklearn.metrics import average_precision_score
from imblearn.over_sampling import SMOTE, BorderlineSMOTE, ADASYN, SMOTENC
from imblearn.combine import SMOTETomek, SMOTEENN
from imblearn.pipeline import Pipeline as ImbPipeline

# --- SMOTE FROM SCRATCH: the whole algorithm is five lines. ---
def smote_manual(Xmin: np.ndarray, n_new: int, k: int = 5, seed: int = 0):
    rng = np.random.default_rng(seed)
    nn = NearestNeighbors(n_neighbors=min(k + 1, len(Xmin))).fit(Xmin)
    _, idx = nn.kneighbors(Xmin)
    out = np.empty((n_new, Xmin.shape[1]))
    for i in range(n_new):
        a = rng.integers(len(Xmin))                    # a real minority point
        b = idx[a][rng.integers(1, idx.shape[1])]      # one of its neighbours
        lam = rng.random()                             # position on the segment
        out[i] = Xmin[a] + lam * (Xmin[b] - Xmin[a])
    return out

real = np.array([[1.0, 1.0], [1.0, 2.0], [2.0, 1.5], [8.0, 8.0]])
fake = smote_manual(real, 6, k=3)
print("real minority points:\n", real)
print("\nsynthesised points:\n", np.round(fake, 2))
print("\nsome lie on the segment toward the isolated point at (8, 8),")
print("which may be majority territory or a mislabelled row.")

# --- FAILURE MODE 1: interpolating across a gap between clusters. ---
rng = np.random.default_rng(0)
cluster_a = rng.normal([0, 0], 0.4, (20, 2))
cluster_b = rng.normal([10, 10], 0.4, (20, 2))
two_clusters = np.vstack([cluster_a, cluster_b])
synth = smote_manual(two_clusters, 300, k=5, seed=1)
midzone = ((synth[:, 0] > 3) & (synth[:, 0] < 7)).sum()
print(f"\nsynthetic points landing in the empty corridor between clusters:"
      f" {midzone} of 300 ({midzone/300:.0%})")
print("no real minority example lives there, and the majority class does")

# --- FAILURE MODE 2: one mislabelled point becomes a line of them. ---
mostly_tight = np.vstack([rng.normal([0, 0], 0.3, (30, 2)), [[9.0, 9.0]]])
noisy = smote_manual(mostly_tight, 400, k=5, seed=2)
far = (np.linalg.norm(noisy, axis=1) > 3).sum()
print(f"\nwith ONE mislabelled point, {far} of 400 synthetic points"
      f" ({far/400:.0%}) reach toward it")

# --- FAILURE MODE 3: categorical features have no midpoint. ---
mixed = pd.DataFrame({"age": [25.0, 60.0], "colour_code": [0.0, 2.0]})
print("\nnaive interpolation of a categorical code:")
print(f"  midpoint of red(0) and blue(2) is {mixed['colour_code'].mean():.1f}"
      " = green, a category neither point had")
# SMOTENC interpolates the numeric columns and votes on the categorical ones.

# --- COMPARE THE VARIANTS HONESTLY, on real class distributions. ---
X, y = make_classification(n_samples=5000, n_features=18, n_informative=6,
                           n_clusters_per_class=3, weights=[0.96, 0.04],
                           flip_y=0.03, random_state=1)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0,
                                      stratify=y)
print(f"\ntrain {Counter(ytr)}   test {Counter(yte)} (untouched)")

candidates = {
    "baseline": None,
    "class_weight balanced": "weight",
    "SMOTE": SMOTE(random_state=0, k_neighbors=5),
    "Borderline-SMOTE": BorderlineSMOTE(random_state=0),
    "ADASYN": ADASYN(random_state=0),
    "SMOTE + Tomek": SMOTETomek(random_state=0),
    "SMOTE + ENN": SMOTEENN(random_state=0),
}
cv = StratifiedKFold(5, shuffle=True, random_state=0)
print("\ncross-validated PR AUC (resampling inside each fold):")
for name, sampler in candidates.items():
    if sampler == "weight":
        pipe = ImbPipeline([("m", RandomForestClassifier(
            n_estimators=150, class_weight="balanced_subsample", random_state=0))])
    elif sampler is None:
        pipe = ImbPipeline([("m", RandomForestClassifier(
            n_estimators=150, random_state=0))])
    else:
        pipe = ImbPipeline([("res", sampler),
                            ("m", RandomForestClassifier(n_estimators=150,
                                                         random_state=0))])
    s = cross_validate(pipe, Xtr, ytr, cv=cv, scoring="average_precision")
    print(f"  {name:24s} {s['test_average_precision'].mean():.3f}"
          f"  +/- {s['test_average_precision'].std():.3f}")
# Read the standard deviations. Differences smaller than the spread are not
# results. Class weighting frequently matches or beats every SMOTE variant.

# --- MIXED TYPES: SMOTENC ---
Xm = pd.DataFrame(X[:, :4], columns=["n1", "n2", "n3", "n4"])
Xm["cat"] = pd.Series(np.where(X[:, 4] > 0, "A", "B"))
nc = SMOTENC(categorical_features=[4], random_state=0)
Xr, yr = nc.fit_resample(Xm, y)
print(f"\nSMOTENC output categories: {sorted(Xr['cat'].unique())}"
      "  <- no invented category")

# --- THE PIPELINE POSITION IS NOT OPTIONAL ---
print("\nWRONG: SMOTE(...).fit_resample(X, y) then train_test_split")
print("RIGHT: ImbPipeline([('res', SMOTE()), ('m', model)]) inside CV")
wrong_X, wrong_y = SMOTE(random_state=0).fit_resample(X, y)
a, b, c, d = train_test_split(wrong_X, wrong_y, test_size=0.3, random_state=0,
                              stratify=wrong_y)
m = RandomForestClassifier(n_estimators=120, random_state=0).fit(a, c)
print(f"  resample-then-split PR AUC "
      f"{average_precision_score(d, m.predict_proba(b)[:, 1]):.3f} (inflated)")
# Synthetic points derived from a training row can land in the test half,
# which is a near-duplicate of something the model has already seen.

  HOW SMOTE PLACES A POINT

    minority point a ----------o--------- neighbour b
                          new point at a + lam*(b-a), lam ~ U(0,1)

  WHEN THAT IS SENSIBLE           WHEN IT IS NOT

    P P P                          P P                    P P
    P P P   dense cluster,         P P   <-- corridor --> P P
    P P P   interior is            N N N N N N N N N N N N
            plausibly P            the segment crosses majority territory

  NOISE AMPLIFICATION

    before:  P P P P               N N N [P] N N N     one mislabelled P
    after:   P P P P  p p p p p p p p p [P] N N N      a bridge of synthetic
                                                       points into the N region

  VARIANT SELECTION

  variant            seeds from                     use when
  --------------------------------------------------------------------------
  SMOTE              all minority points            default attempt
  Borderline-SMOTE   minority points near the edge   boundary is the problem
  ADASYN             harder neighbourhoods           density varies a lot
  SVM-SMOTE          support vectors                 SVM-shaped problems
  SMOTE-Tomek/ENN    synthesise, then clean          noisy labels
  SMOTENC            mixed numeric + categorical     any real tabular data

  THE BOTTOM LINE
    SMOTE cannot create information. It redistributes the information in the
    minority points you already have. Compare it against class weighting and
    let the held-out score decide.

Key Takeaways

• SMOTE interpolates between a minority point and one of its neighbours, which densifies the minority region instead of duplicating points.
• When the minority class has separate clusters, interpolation places synthetic points in majority territory.
• A single mislabelled minority point seeds a whole line of synthetic points reaching into the majority region.
• Categorical features have no midpoint and high dimensions make nearest neighbours unreliable, so SMOTE-NC and caution are required.
• SMOTE creates no new information and frequently fails to beat class weighting, so benchmark it rather than assuming it.

🧪 Practice

1. Implement SMOTE in under ten lines and use it on a minority class made of two well-separated clusters, then count how many synthetic points fall in the empty corridor.
2. Add one mislabelled minority point deep in majority territory and measure how many synthetic points reach toward it.
3. Benchmark baseline, class weighting, SMOTE, and Borderline-SMOTE with cross-validated PR AUC, reporting standard deviations, and state whether any difference exceeds the noise.
4. Interview: When would you actively avoid SMOTE? (Hint: think about the feature types, the dimensionality, and whether anything downstream consumes the predicted probability.)

Class Weighting

Class weighting attacks imbalance from the loss function rather than the data. Instead of changing how many minority examples exist, you change how much each one costs to get wrong. Misclassifying a minority example incurs a larger penalty, so the optimiser stops finding it acceptable to ignore the whole class.

This is usually the right first thing to try, and the reasons are worth listing because they are all practical rather than theoretical. You keep every row of data. You add no duplicated or synthesised points, so nothing is invented. Training time does not increase. Calibration is distorted less than by resampling. And it is one parameter, which means it is trivial to include in a hyperparameter search alongside everything else.

The mechanics are straightforward. The loss for each example is multiplied by its class's weight. The common default, "balanced", sets each class's weight inversely proportional to its frequency, so that the total weight of each class is equal. In formula terms the weight for class c is n / (k * n_c), where n is the total row count, k the number of classes, and n_c the count in that class.

There is an important conceptual relationship: class weighting is approximately equivalent to oversampling by the same factor, without materialising the extra rows. A weight of 30 on the minority class has a similar effect on the loss as duplicating each minority row thirty times. The weighted version is cheaper, does not create repeated identical points, and keeps the dataset honest.

Two refinements matter in practice.

First, "balanced" is not necessarily optimal. It equalises the classes, which is a specific and arbitrary target. What you actually want is the weight ratio that reflects the real cost ratio of the two error types, or failing that, the ratio that maximises your chosen metric. Tune it.

Second, sample weights are more general than class weights. Most implementations accept a per-row weight, which lets you encode a per-example cost: a fraud of ten thousand matters more than one of fifty. This is cost-sensitive learning proper, and it is frequently the most direct formulation of the actual business problem, which the class-balance framing only approximates.

Support for weighting is near-universal: scikit-learn's class_weight parameter on most classifiers, scale_pos_weight in XGBoost and LightGBM, sample_weight in fit for nearly everything, and weighted loss functions in every deep learning framework. Focal loss, common in object detection, goes further by down-weighting easy examples regardless of class, which addresses a related problem.

One caution: weighting changes the model's output scale, so predicted probabilities are no longer calibrated to the true base rate. If downstream logic consumes probabilities rather than rankings or labels, recalibrate on unweighted data.

import numpy as np
import pandas as pd
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split, StratifiedKFold, cross_validate
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.utils.class_weight import compute_class_weight, compute_sample_weight
from sklearn.metrics import (average_precision_score, recall_score,
                             precision_score, brier_score_loss, confusion_matrix)

X, y = make_classification(n_samples=8000, n_features=16, n_informative=6,
                           weights=[0.96, 0.04], flip_y=0.02, random_state=2)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0,
                                      stratify=y)

# --- WHAT "balanced" COMPUTES ---
classes = np.unique(ytr)
w = compute_class_weight("balanced", classes=classes, y=ytr)
print(dict(zip(classes.tolist(), w.round(2).tolist())))
n, k = len(ytr), len(classes)
for c in classes:
    print(f"  class {c}: n/(k*n_c) = {n}/({k}*{(ytr == c).sum()})"
          f" = {n / (k * (ytr == c).sum()):.2f}")
print(f"imbalance ratio: {(ytr == 0).sum() / (ytr == 1).sum():.1f}:1")

# --- THE EFFECT ON PREDICTIONS ---
plain = LogisticRegression(max_iter=3000).fit(Xtr, ytr)
weighted = LogisticRegression(max_iter=3000, class_weight="balanced").fit(Xtr, ytr)

print("\nat the default threshold of 0.5:")
for name, m in [("unweighted", plain), ("balanced", weighted)]:
    p = m.predict(Xte)
    tn, fp, fn, tp = confusion_matrix(yte, p).ravel()
    print(f"  {name:11s} flagged {p.sum():>4}  TP {tp:>3}  FP {fp:>4}"
          f"  recall {recall_score(yte, p):.3f}"
          f"  precision {precision_score(yte, p, zero_division=0):.3f}")
# Weighting moves the model's operating point toward recall without any
# change to the data. It is doing at training time what a lower threshold
# does at scoring time, which is why the two are alternatives.

# --- RANKING QUALITY BARELY CHANGES; THE OPERATING POINT DOES. ---
print("\nPR AUC, which is threshold-free:")
for name, m in [("unweighted", plain), ("balanced", weighted)]:
    ap = average_precision_score(yte, m.predict_proba(Xte)[:, 1])
    print(f"  {name:11s} {ap:.3f}")
print("weighting mainly relocates the decision, not the ranking")

# --- "balanced" IS ONE CHOICE. TUNE THE RATIO. ---
cv = StratifiedKFold(5, shuffle=True, random_state=0)
print("\nsweeping the minority weight:")
for mult in (1, 2, 5, 10, 24, 50, 100):
    pipe = LogisticRegression(max_iter=3000, class_weight={0: 1, 1: mult})
    s = cross_validate(pipe, Xtr, ytr, cv=cv,
                       scoring=["average_precision", "recall"])
    print(f"  weight {mult:>3}  PR AUC {s['test_average_precision'].mean():.3f}"
          f"   recall {s['test_recall'].mean():.3f}")
print(f"  ('balanced' corresponds to weight ~{w[1]:.0f})")
# PR AUC is often flat across this whole sweep while recall climbs steadily.
# That is the signal that you are choosing an operating point, not improving
# the model, and a threshold would have done the same thing for free.

# --- WEIGHTING VS OVERSAMPLING: two routes to the same place. ---
from imblearn.over_sampling import RandomOverSampler
from imblearn.pipeline import Pipeline as ImbPipeline
print("\nequivalent formulations, compared:")
for name, pipe in {
    "class_weight balanced": LogisticRegression(max_iter=3000,
                                                class_weight="balanced"),
    "random oversample to 1:1": ImbPipeline([
        ("res", RandomOverSampler(random_state=0)),
        ("m", LogisticRegression(max_iter=3000))]),
}.items():
    s = cross_validate(pipe, Xtr, ytr, cv=cv, scoring="average_precision")
    print(f"  {name:26s} {s['test_average_precision'].mean():.3f}")
print("similar scores; weighting needs no extra rows and no synthetic data")

# --- SAMPLE WEIGHTS: per-row cost, which is the real business problem. ---
rng = np.random.default_rng(0)
amount = rng.gamma(2, 400, len(ytr)).round(2)      # value at risk per row
# A missed large fraud costs more than a missed small one. Encode that.
cost_weight = np.where(ytr == 1, amount / amount.mean(), 1.0)

rf = RandomForestClassifier(n_estimators=200, random_state=0)
rf.fit(Xtr, ytr, sample_weight=cost_weight)
proba = rf.predict_proba(Xte)[:, 1]
print(f"\ncost-weighted random forest PR AUC: "
      f"{average_precision_score(yte, proba):.3f}")
print("this optimises expected loss rather than the count of errors, which")
print("is usually what the business actually asked for.")

# compute_sample_weight builds class-balanced row weights for you.
auto = compute_sample_weight("balanced", ytr)
print(f"\nauto sample weights: minority {auto[ytr == 1][0]:.2f},"
      f" majority {auto[ytr == 0][0]:.2f}")

# --- CALIBRATION COST ---
print(f"\ntrue positive rate {yte.mean():.4f}")
for name, m in [("unweighted", plain), ("balanced", weighted)]:
    p = m.predict_proba(Xte)[:, 1]
    print(f"  {name:11s} mean predicted {p.mean():.4f}"
          f"   Brier {brier_score_loss(yte, p):.5f}")
# The weighted model's probabilities are inflated toward the minority class.
# Fine if you need a ranking or a label. Wrong if you multiply by a loss.

# --- LIBRARY-SPECIFIC PARAMETERS ---
print("""
where to set it:
  scikit-learn      class_weight="balanced" | {0: 1, 1: 24} | sample_weight in fit
  XGBoost           scale_pos_weight=neg/pos, or sample_weight
  LightGBM          class_weight / is_unbalance / scale_pos_weight
  CatBoost          class_weights / auto_class_weights="Balanced"
  PyTorch           nn.BCEWithLogitsLoss(pos_weight=...)
                    nn.CrossEntropyLoss(weight=...)
  TensorFlow/Keras  model.fit(..., class_weight={0: 1.0, 1: 24.0})
""")

  WEIGHTING VS RESAMPLING

                        class weighting        oversampling
  data kept             all                    all
  rows added            none                   many duplicates
  training time         unchanged              grows with the new rows
  memory                unchanged              grows
  duplicate points      none                   many, invites memorisation
  calibration           distorted              distorted, usually more
  tunable               one parameter          sampling strategy
  equivalent to         oversampling by w      weighting by the duplication

  WHAT WEIGHTING ACTUALLY MOVES

    loss = sum over rows of  w[class(row)] * per_row_loss(row)

    raising w on the minority makes ignoring it expensive, so the fitted
    decision surface shifts toward the minority. The RANKING changes little;
    the OPERATING POINT changes a lot. A threshold does the same thing.

  ESCALATION ORDER, RESTATED

    metric -> threshold -> class weight -> sample weight (real costs)
           -> resampling -> synthetic data -> collect more positives

Key Takeaways

• Class weighting changes the loss rather than the data, keeping every row and adding no synthetic points.
• A weight of w on the minority class is approximately oversampling by w without materialising the rows.
• The "balanced" setting equalises class totals, which is arbitrary; tune the ratio against your metric or the real cost ratio.
• Sample weights let you encode a per-row cost, which is usually a closer match to the business problem than class balance.
• Weighting inflates predicted probabilities toward the minority class, so recalibrate if anything downstream consumes probabilities.

🧪 Practice

1. Compute balanced class weights by hand from the class counts and confirm they match the library's output.
2. Sweep the minority class weight across a wide range and report both precision-recall area and recall, then say which one the sweep actually changed.
3. Fit a model with per-row weights proportional to a monetary amount and compare its flagged set against a class-balanced model's.
4. Interview: Class weighting and threshold tuning both shift the precision-recall trade-off. When would you prefer each? (Hint: one changes what the model learns, the other changes only how you read it.)

Log and Power Transformations

Skew is the other kind of imbalance: not in the label but in a numeric variable's own distribution. Income, house prices, page views, transaction amounts, city populations, and file sizes all share a shape. Most values are modest, a few are enormous, and the mean sits well above the median. This is right skew, and it breaks a surprising number of default assumptions.

What it breaks specifically: the mean stops describing a typical value; linear regression's residuals become non-constant in variance, which invalidates its standard errors; distance-based methods are dominated by the tail; visualisations become unreadable because 99% of the data crushes into the first tick; and a squared-error loss chases the extremes because a large error on a large value dominates the objective.

A log transform fixes much of this at once. It converts multiplicative structure into additive structure, which is the deeper reason it works so well. Many real quantities vary by ratios rather than by absolute amounts: a salary rise is a percentage, a population grows by a rate, a file is twice as large. For such quantities the log scale is the natural one, and differences on it are ratios on the original scale.

Two practical details. log is undefined at zero and negative values, so use log1p, which computes log(1 + x) and is exact at zero. For data with genuine negative values, the inverse hyperbolic sine transform behaves like a log for large magnitudes and is defined everywhere, including negatives.

The other transforms:

Square root is a gentler compression than log and is the natural choice for count data, where variance grows with the mean.

Box-Cox is a family parameterised by lambda, which is estimated from the data to maximise normality. Lambda of 0 is the log, 0.5 is the square root, 1 is no transform. It requires strictly positive input.

Yeo-Johnson is the same idea extended to handle zero and negative values, which makes it the usable default when you want an automatically chosen power transform.

Quantile transformation maps to ranks and then to a normal or uniform distribution. It will normalise essentially anything, and it discards the original spacing entirely, which is a large price.

Two cautions that matter most.

Interpretation changes. A coefficient on log(income) is not an effect per dollar. In a log-linear model it is approximately the proportional change in the outcome per proportional change in income, which is an elasticity. If you transform the target, back-transforming a prediction with exp gives the predicted median, not the mean, because the expectation does not pass through a non-linear function. Reporting it as the mean understates systematically, and the correction requires adding a variance term.

Trees do not need it. A monotonic transform does not change the ordering of values, and a decision tree only uses ordering. Log-transforming a feature before gradient boosting changes nothing. It can still help if you transform the target, because that changes the loss function's geometry, which is a different thing entirely.

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.preprocessing import PowerTransformer, QuantileTransformer
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.metrics import mean_absolute_error

rng = np.random.default_rng(6)
income = pd.Series(rng.lognormal(10.5, 0.9, 3000).round(0), name="income")

print(f"mean   {income.mean():>12,.0f}")
print(f"median {income.median():>12,.0f}")
print(f"max    {income.max():>12,.0f}")
print(f"skew   {income.skew():>12.2f}   kurtosis {income.kurtosis():.2f}")
print(f"the mean sits at the {stats.percentileofscore(income, income.mean()):.0f}th"
      " percentile: most people earn less than 'average'")

# --- COMPARE THE TRANSFORMS ---
def describe(s: pd.Series) -> dict:
    return {"skew": round(s.skew(), 2), "kurtosis": round(s.kurtosis(), 2),
            "shapiro_p": round(stats.shapiro(s.sample(500, random_state=0))[1], 4)}

bc = PowerTransformer(method="box-cox").fit(income.to_frame())
yj = PowerTransformer(method="yeo-johnson").fit(income.to_frame())
qt = QuantileTransformer(output_distribution="normal", n_quantiles=500,
                        random_state=0).fit(income.to_frame())

variants = {
    "raw": income,
    "log1p": np.log1p(income),
    "sqrt": np.sqrt(income),
    "box-cox": pd.Series(bc.transform(income.to_frame()).ravel()),
    "yeo-johnson": pd.Series(yj.transform(income.to_frame()).ravel()),
    "quantile-normal": pd.Series(qt.transform(income.to_frame()).ravel()),
}
print("\n", pd.DataFrame({k: describe(v) for k, v in variants.items()}).T.to_string())
print(f"\nBox-Cox chose lambda = {bc.lambdas_[0]:.3f}"
      "   (0 = log, 0.5 = sqrt, 1 = no transform)")

# --- ZERO AND NEGATIVE VALUES ---
with_zero = pd.Series([0.0, 1.0, 10.0, 100.0])
print(f"\nlog(0) = {np.log(with_zero[0])}   log1p(0) = {np.log1p(with_zero[0])}")
with_neg = pd.Series([-500.0, -1.0, 0.0, 1.0, 500.0])
asinh = np.arcsinh(with_neg)                  # log-like, defined everywhere
print("\narcsinh handles negatives:")
print(pd.DataFrame({"x": with_neg, "arcsinh": asinh.round(3)}).to_string(index=False))

# --- MULTIPLICATIVE BECOMES ADDITIVE, which is the real reason it works. ---
base = np.array([30_000.0, 60_000.0, 120_000.0])       # each double the last
print("\nequal RATIOS become equal DIFFERENCES on the log scale:")
print(pd.DataFrame({"income": base, "ratio_to_prev": [np.nan, 2.0, 2.0],
                    "log": np.log(base).round(3),
                    "log_diff": [np.nan, *np.diff(np.log(base)).round(3)]}
                   ).to_string(index=False))

# --- HETEROSCEDASTICITY: the reason regression assumptions break. ---
x = rng.uniform(1, 10, 1500)
y_mult = np.exp(1.0 + 0.3 * x + rng.normal(0, 0.35, 1500))   # multiplicative noise

raw_fit = LinearRegression().fit(x.reshape(-1, 1), y_mult)
raw_resid = y_mult - raw_fit.predict(x.reshape(-1, 1))
log_fit = LinearRegression().fit(x.reshape(-1, 1), np.log(y_mult))
log_resid = np.log(y_mult) - log_fit.predict(x.reshape(-1, 1))

lo, hi = x < np.median(x), x >= np.median(x)
print("\nresidual spread in the low vs high half of x:")
print(f"  raw target: {raw_resid[lo].std():8.2f} vs {raw_resid[hi].std():8.2f}"
      f"   ratio {raw_resid[hi].std()/raw_resid[lo].std():5.2f}  <- fans out")
print(f"  log target: {log_resid[lo].std():8.3f} vs {log_resid[hi].std():8.3f}"
      f"   ratio {log_resid[hi].std()/log_resid[lo].std():5.2f}  <- constant")
print("constant residual variance is what the standard errors assume")

# --- BACK-TRANSFORMING GIVES THE MEDIAN, NOT THE MEAN. ---
resid_var = log_resid.var(ddof=1)
pred_log = log_fit.predict(x.reshape(-1, 1))
naive = np.exp(pred_log)                                  # predicted MEDIAN
corrected = np.exp(pred_log + resid_var / 2)              # predicted MEAN
print(f"\nactual mean of y            {y_mult.mean():12,.1f}")
print(f"mean of exp(prediction)     {naive.mean():12,.1f}   <- biased low")
print(f"with the variance correction {corrected.mean():12,.1f}   <- unbiased")
print("E[exp(Z)] != exp(E[Z]); the gap grows with the residual variance")

# --- TRANSFORMING A FEATURE: helps linear models, not trees. ---
Xf = pd.DataFrame({"income": income})
target = (0.8 * np.log(income) + rng.normal(0, 0.4, len(income))).values

for model_name, model in [("linear regression", LinearRegression()),
                          ("gradient boosting", GradientBoostingRegressor(
                              random_state=0, n_estimators=120))]:
    raw_s = cross_val_score(model, Xf, target, cv=5, scoring="r2").mean()
    log_s = cross_val_score(model, np.log1p(Xf), target, cv=5, scoring="r2").mean()
    print(f"\n{model_name}: raw R2 {raw_s:.3f}   log feature R2 {log_s:.3f}")
print("\nthe tree model is unchanged: a monotonic transform cannot alter the")
print("ordering of values, and a split only uses the ordering.")

# --- TRANSFORMING THE TARGET changes the loss, so trees DO respond. ---
ytr_skew = np.exp(2 + 0.5 * x + rng.normal(0, 0.4, len(x)))
Xs = x.reshape(-1, 1)
a, b, c, d = train_test_split(Xs, ytr_skew, test_size=0.3, random_state=0)
gb_raw = GradientBoostingRegressor(random_state=0).fit(a, c)
gb_log = GradientBoostingRegressor(random_state=0).fit(a, np.log1p(c))
print(f"\nMAE, raw target        {mean_absolute_error(d, gb_raw.predict(b)):10,.2f}")
print(f"MAE, log target, back  "
      f"{mean_absolute_error(d, np.expm1(gb_log.predict(b))):10,.2f}")
print("squared error on a raw skewed target chases the largest values;")
print("on the log scale each row contributes on a proportional footing.")

# --- INTERPRETATION CHEAT SHEET ---
print("""
  model form              coefficient b means
  --------------------------------------------------------------------------
  y   ~ x                 +1 unit of x -> +b units of y
  y   ~ log(x)            +1% in x     -> +b/100 units of y
  log(y) ~ x              +1 unit of x -> +100*b percent in y (approx.)
  log(y) ~ log(x)         +1% in x     -> +b percent in y     (an elasticity)
""")

  RIGHT SKEW AND WHAT LOG DOES TO IT

  raw income
   |#
   |##
   |####
   |#######_______________________________________________
   +--------------------------------------------------------> value
   mean sits here ^ (above ~70% of the data)

  log(income)
   |        ####
   |      ########
   |    ############
   |  ################
   +--------------------------------------------------------> log value
   roughly symmetric; mean and median coincide

  CHOOSING A TRANSFORM

  transform      handles 0 / negatives   strength    interpretable
  --------------------------------------------------------------------
  log1p          0 yes, negatives no     strong      yes (ratios)
  sqrt           0 yes, negatives no     moderate    partly (counts)
  arcsinh        both yes                strong      roughly like log
  Box-Cox        strictly positive only  tuned       via lambda
  Yeo-Johnson    both yes                tuned       barely
  quantile       both yes                total       no

  WHEN IT MATTERS

  situation                            transform the feature? the target?
  ---------------------------------------------------------------------------
  linear regression, skewed predictor  YES                    -
  residual variance fans out           -                      YES
  distance-based model                 YES                    -
  tree ensemble                        no (no effect)         often yes
  reporting a typical value            report the median instead of the mean

  THE BACK-TRANSFORM TRAP
    exp(mean of logs) = geometric mean = predicted MEDIAN.
    For the predicted MEAN add sigma^2/2 before exponentiating.

Key Takeaways

• Right skew makes the mean unrepresentative, inflates residual variance, and lets the tail dominate distance and squared-error objectives.
• A log transform converts multiplicative structure into additive structure, which is why it suits quantities that vary by ratios.
• Use log1p for data containing zero and the inverse hyperbolic sine for data containing negatives; Yeo-Johnson is the general automatic option.
• Exponentiating a prediction from a log-target model gives the median, and recovering the mean requires a variance correction.
• Transforming a feature does nothing for tree ensembles, while transforming the target changes the loss and often helps them.

🧪 Practice

1. Report skew, kurtosis, and a normality test statistic for a log-normal column raw, log-transformed, and Yeo-Johnson transformed.
2. Fit a regression on a multiplicatively noisy target with and without a log transform and compare residual spread across the low and high halves of the predictor.
3. Back-transform predictions from a log-target model with and without the variance correction, and compare both against the actual mean.
4. Show that log-transforming a feature leaves a gradient boosting model's cross-validated score unchanged while changing a linear model's.
5. Interview: A stakeholder asks why you report median order value rather than the mean. What do you say in two sentences? (Hint: describe where the mean falls in the distribution and what a single large order does to it.)

---

<a id="7-exploratory-data-analysis-and-visualization"></a>