Missing values and messy data

What this chapter covers

Back in the threshold chapter I put a footnote on the Pima data and moved on: a few of those glucose readings are zero, which is impossible, so treat them with suspicion. This chapter is me paying that footnote off in full. It turns out the suspicion was the whole story. Nearly half the insulin column is zeros, and a zero insulin reading means the patient is dead, not that their pancreas is idle. Those aren't measurements. They're missing values wearing a measurement's clothes, and every model in this book will happily swallow them and quietly get worse.

So before any algorithm runs, someone has to find the holes in the data and decide what to do about them. That's this chapter. We build the tools from scratch in NumPy and pandas: a detector that catches missing values including the disguised kind, three ways to fill them in (mean, median, mode) and an argument for when each one is right, and two outlier rules that disagree on purpose. Then we let scikit-learn's SimpleImputer fill the same gaps and check we land on identical numbers. The dataset is the Pima set again, chosen because its dirt is famous and instructive, and the reveal is a bar chart that counts exactly how many impossible zeros are hiding in each column. It's more than you'd guess.

A bit of history

The formal study of missing data is younger than you'd think. For most of the 20th century the standard move was to throw away any row with a hole in it — "complete-case analysis" — and hope the rows you deleted looked like the rows you kept. Donald Rubin killed that hope in a 1976 paper, "Inference and Missing Data," which gave the field the vocabulary it still uses. Rubin's point was that you can't decide how to handle a gap until you understand why the value is missing, and he named three cases. Missing completely at random (MCAR): the gaps have nothing to do with anything, a dropped test tube. Missing at random (MAR): the gaps depend on data you do have, like older patients skipping a test more often. And missing not at random (MNAR): the gap depends on the missing value itself, like a scale that refuses to report weights above 300 pounds. The three cases need three different treatments, and the dangerous one is MNAR, because the absence is data and deleting it throws away signal. Hold that thought — the Pima insulin column is closer to MNAR than anyone likes to admit, and it bites us at the end of the chapter.

The other half of this chapter's lineage is John Tukey and exploratory data analysis. Tukey's 1977 book "Exploratory Data Analysis" made an argument that now sounds obvious and wasn't: before you fit anything, look at the data. Plot it, summarize it with statistics that a few crazy values can't wreck, find the points that don't belong. The boxplot is his, and so is the fence rule we'll use for outliers — the median and the interquartile range instead of the mean and standard deviation, precisely because the mean and standard deviation are the first things an outlier corrupts. Rubin tells you to respect the holes; Tukey tells you to look before you leap. This chapter is both of them, in code.

The intuition

Here's the trap in one sentence: a missing value doesn't have to be blank to be missing. In a clean world, the absence of a measurement is stored as NaN, a null, an empty cell — something an honest detector can find. In the real world, whoever built the table needed a number to put there, and zero was sitting right next to their hand. So the missing insulin readings became 0, the missing blood pressures became 0, and now your dataset is full of values that pass every type check, load without a warning, and are complete fiction.

The tell is domain knowledge, not code. A body mass index of zero is a person with no mass. A diastolic blood pressure of zero is a corpse. Once you know the column can't physically be zero, every zero in it is a confession that the measurement never happened. This is the single most common way real data lies to you, and it's invisible unless you go looking. Here's what going looking turns up in the Pima set — the count of impossible zeros in each column that has them:

The number on each bar is the raw count of impossible zeros; the axis is that as a percentage of all 768 patients. Insulin is 374 zeros — 48.7% of the column is fake. Skin thickness is 227, almost 30%. Blood pressure, BMI, and glucose have smaller but real holes. Add it up and 652 individual cells in this "clean, ready-to-model" benchmark dataset are missing values in disguise. If you trained on the raw table you'd be teaching a model that half your patients secretly produce no insulin, which is not a fact about diabetes. It's a fact about data entry.

The math

Three small pieces of math run this chapter, and none of them is heavy. First, the fills. To impute a column you replace every missing entry with one summary of the values that aren't missing. The mean is the average:

where the sum runs over the $n$ observed (non-missing) entries only. The median is the middle value of those same entries once you sort them:

with $x_{(k)}$ the $k$-th smallest observed value. The mean and the median agree on a symmetric column and disagree on a skewed one, and the size of that disagreement is the whole argument for which to use — I'll make it concrete on the insulin column in a moment.

Second, the two outlier rules, which are just two different definitions of "too far from the middle." The z-score rule measures distance in standard deviations from the mean:

$\sigma$ is the column's standard deviation, and the usual threshold is 3. Tukey's rule measures distance from the quartiles in units of the interquartile range instead:

where $Q_1$ and $Q_3$ are the 25th and 75th percentiles and $\text{IQR} = Q_3 - Q_1$ is the width of the middle half of the data. The two rules disagree, and the reason they disagree is the point: the z-score's own ingredients, $\bar{x}$ and $\sigma$, are inflated by the very outliers it's trying to catch, so on a skewed column it under-flags. The IQR rule is built from percentiles, which a handful of extreme values can't move, so it keeps working where the z-score goes soft.

What it's good at, what it isn't

Imputation's virtue is that it keeps your rows. Delete every patient with a missing insulin reading and you've thrown away 374 of 768 people — half your data, gone, to patch one column. Fill the gaps instead and every row survives with every other measurement intact. That's the case for imputing, and on a dataset this small it's decisive. The cost is that you're inventing numbers, and a filled value is a guess dressed as a fact. Fill with a constant like the median and you shrink the column's variance, you weaken its correlations, and you plant a suspicious spike of identical values that a later model can key on for the wrong reasons. Mean and median imputation are the cheap, honest-enough default; they are not free.

The judgment call is drop versus impute, and it comes down to how much is missing and why. A column that's 1% missing you impute without a second thought. A column that's 48% missing — insulin, here — is a real dilemma: half of it is invented after you fill it, and no summary statistic can conjure the signal that was never recorded. Sometimes the right answer for a column that empty is to drop the column, not the rows. And sometimes, when the missingness itself is informative (Rubin's MNAR), the best move is to impute a placeholder and add a flag marking which values were missing, so the model can learn from the absence instead of being fooled by the fill. That last trick is the one that saves us at the end.

The data

Same Pima Indians Diabetes set as week 1: 768 patients, eight measurements each, a binary outcome for a diabetes diagnosis. National Institute of Diabetes and Digestive and Kidney Diseases collected it; it's been the standard teaching set for missing-data problems for decades, precisely because of the mess we just charted. Five of the eight feature columns encode missingness as zero: glucose, blood pressure, skin thickness, insulin, and BMI. The other three don't — pregnancies of zero is a real and common value, and age and the diabetes-pedigree score are never zero — so a blanket "treat all zeros as missing" would corrupt the pregnancies column. The detector has to know which columns to suspect. That knowledge is domain knowledge, hardcoded, because there's no statistical test that can tell you a zero is impossible; only a doctor can.

We'll split once — 576 train, 192 test, stratified on the outcome, seed 0 — and that split is fixed for every accuracy on this page. Everything that follows fits its statistics on the training rows only. The insulin column is the one we'll watch closely, because it's the worst offender and because its shape makes the mean-versus-median argument impossible to miss.

Build it, one function at a time

The first job is detection, and it's the only function in the file that needs to know something about the world. Real NaNs are easy; the disguised zeros need a per-column decision about whether zero is even possible:

def missing_mask(X, cols_missing_zero):
    """Boolean mask, True where a value is missing — real NaN OR a disguised zero.

    X is (N, D). `cols_missing_zero` is a boolean length-D vector marking the
    columns where 0 secretly means "missing". For those columns a 0 counts as
    missing; for the rest, only an actual NaN does. The whole chapter turns on
    this one idea: missingness can hide in plain sight as a legal-looking number.
    """
    mask = np.isnan(X)
    mask[:, cols_missing_zero] |= (X[:, cols_missing_zero] == 0)
    return mask


def missing_counts(X, names, cols_missing_zero):
    """Per-column (count, percent) of missing values, using missing_mask."""
    mask = missing_mask(X, cols_missing_zero)
    n = X.shape[0]
    return [{"feature": names[j], "count": int(mask[:, j].sum()),
             "percent": round(100.0 * mask[:, j].sum() / n, 2)}
            for j in range(X.shape[1])]

The mask is the foundation everything else stands on. missing_mask starts from the honest nulls with np.isnan, then flips on the impossible zeros for exactly the columns we flagged as suspect, and leaves pregnancies alone. missing_counts just tallies it per column — that's the function behind the reveal chart above. Once you have the mask, the fills are almost trivial, and they all share one non-negotiable rule: compute the summary from the observed values only, never including the holes you're about to fill.

def impute_fit(X, mask, strategy="median"):
    """Learn the fill value for each column from the NON-missing rows only.

    For every column we ignore the missing entries and compute one summary of
    what's left: the mean, the median, or the mode (most frequent value). Median
    is the default on purpose — it sits in the bulk of a skewed column where the
    mean gets dragged into the empty tail. Returns a length-D vector of fills.
    """
    D = X.shape[1]
    fill = np.empty(D)
    for j in range(D):
        col = X[~mask[:, j], j]                       # observed values only
        if col.size == 0:                             # whole column missing
            fill[j] = 0.0
        elif strategy == "mean":
            fill[j] = col.mean()
        elif strategy == "median":
            fill[j] = np.median(col)
        elif strategy == "mode":
            vals, counts = np.unique(col, return_counts=True)
            fill[j] = vals[counts.argmax()]           # first max on ties, like sklearn
        else:
            raise ValueError(f"unknown strategy {strategy!r}")
    return fill


def impute_transform(X, mask, fill):
    """Replace every masked entry with its column's fitted fill value."""
    out = X.copy()
    rows, cols = np.where(mask)
    out[rows, cols] = fill[cols]
    return out

impute_fit reads one number per column off the non-missing entries — a mean, a median, or a mode — and impute_transform drops that number into every masked cell. Splitting fit from transform is the same discipline as the scaling chapter: you fit on train, you apply the frozen fill to test, and the test set's own distribution never leaks into the numbers you filled with. Median is the default argument on purpose, for a reason the animation is about to make visual. Now the piece that turns a liability into a feature — the was-missing indicator:

def missing_indicator(mask, cols_missing_zero):
    """Binary "was this value imputed?" columns for the features that had gaps.

    Imputing throws away one bit of truth: whether the number is real or a
    guess. Appending an indicator per imputed column hands that bit back to the
    model, which can then learn "a missing insulin reading is itself a signal".
    Returns an (N, k) 0/1 array, one column per feature in cols_missing_zero.
    """
    return mask[:, cols_missing_zero].astype(float)

This is one column of zeros and ones per gappy feature, marking which values were imputed. It costs almost nothing and it hands the model a bit of truth that imputation otherwise destroys: whether each number is real or invented. On a dataset where the missingness is informative, this is the difference between a clean-but-lobotomized table and a clean table that kept its signal. Last, the two outlier rules, straight from the math:

def zscore_outliers(col, thresh=3.0):
    """Flag points more than `thresh` standard deviations from the mean.

    z = (x - mean) / std. The classic rule calls anything past |z| = 3 an
    outlier. It assumes a roughly bell-shaped column, so it under-flags a
    skewed one (the mean and std are already inflated by the very tail you're
    hunting) — which is exactly why the IQR rule below is the safer default.
    """
    mu, sigma = col.mean(), col.std()
    if sigma == 0:
        return np.zeros(col.shape, dtype=bool)
    z = (col - mu) / sigma
    return np.abs(z) > thresh


def iqr_outliers(col, k=1.5):
    """Flag points outside [Q1 - k·IQR, Q3 + k·IQR] — Tukey's rule.

    IQR = Q3 - Q1 is the width of the middle half of the data. Tukey's fences
    sit k=1.5 IQRs beyond the quartiles. Built from percentiles, the rule
    doesn't care how long the tail is, so it keeps working on skewed columns
    where the z-score quietly fails.
    """
    q1, q3 = np.percentile(col, [25, 75])
    iqr = q3 - q1
    lo, hi = q1 - k * iqr, q3 + k * iqr
    return (col < lo) | (col > hi)

zscore_outliers flags anything past three standard deviations; iqr_outliers flags anything outside Tukey's fences. On the observed insulin values they disagree, and by a lot — that disagreement is the practical lesson of the two rules, and we read the exact counts off results.json below.

Watch it work

This is the argument for median over mean, made on the real insulin column, one frame at a time. The histogram is insulin readings across all 768 patients; the two vertical lines are the column's mean (red, dashed) and median (green, solid). Watch where they sit as we clean the column.

Frame one is the raw column, and it's a disaster you can see. That amber tower at zero is 374 patients — the disguised-missing — and it drags the mean all the way down to 80 μU/mL and the median down to 30, both of them meaningless because they're computed over a pile of fake zeros. Frame two throws the zeros out and shows the real distribution of the 394 patients who were actually measured, and now you see the shape that decides everything: a right-skewed hump with a long tail of high values. The median lands at 125 μU/mL, sitting in the fat part where most patients actually are. The mean is dragged out to 156 by the tail — 31 units higher, parked in a region where comparatively few real patients live.

Frames three and four are the two fills. Mean imputation drops every one of the 374 missing values at 156, building a purple spike out in the sparse tail, a crowd of invented patients standing where almost no real patient stood. Median imputation drops them at 125, inside the bulk, next to people who look like them. Neither fill is truth — both invent 374 numbers — but the median's invented numbers are far less absurd, and that's the entire rule: on a skewed column, mean imputation piles your guesses where the data isn't. Reach for the median. The more skewed the column, the wider the gap between the two lines, and insulin is skewed enough that you can read the argument off the chart without a single statistic.

The full implementation

The whole file, no library — detection, the three fills, the was-missing indicator, and the two outlier rules, exactly as the checkpoints and the animation used them:

"""Finding and fixing messy data, built from scratch.

The work before any model runs: spot the missing values (including the ones
disguised as zeros), fill them in a way that doesn't lie, and flag the outliers.
Everything here is pure NumPy (+ pandas for loading). The design mirrors the
scaler chapter — a `fit` reads a statistic off the TRAINING rows, a `transform`
applies it — so imputation never leaks test-set numbers into the fill.

Every function below appears in the chapter one step at a time (the
`# region:` markers are what the book's include directives pull in).
"""

import numpy as np
import pandas as pd

# Columns where a value of 0 is physically impossible, so 0 is a stand-in for
# "not recorded" — a disguised missing value. A glucose, blood pressure, skin
# fold, insulin, or BMI of zero would mean a dead patient. Pregnancies of 0 is
# real, Age is never 0, the pedigree score is never 0 — those are left alone.
ZERO_IS_MISSING = ("Glucose", "BloodPressure", "SkinThickness", "Insulin", "BMI")


# region: detect_missing
def missing_mask(X, cols_missing_zero):
    """Boolean mask, True where a value is missing — real NaN OR a disguised zero.

    X is (N, D). `cols_missing_zero` is a boolean length-D vector marking the
    columns where 0 secretly means "missing". For those columns a 0 counts as
    missing; for the rest, only an actual NaN does. The whole chapter turns on
    this one idea: missingness can hide in plain sight as a legal-looking number.
    """
    mask = np.isnan(X)
    mask[:, cols_missing_zero] |= (X[:, cols_missing_zero] == 0)
    return mask


def missing_counts(X, names, cols_missing_zero):
    """Per-column (count, percent) of missing values, using missing_mask."""
    mask = missing_mask(X, cols_missing_zero)
    n = X.shape[0]
    return [{"feature": names[j], "count": int(mask[:, j].sum()),
             "percent": round(100.0 * mask[:, j].sum() / n, 2)}
            for j in range(X.shape[1])]
# endregion


# region: impute
def impute_fit(X, mask, strategy="median"):
    """Learn the fill value for each column from the NON-missing rows only.

    For every column we ignore the missing entries and compute one summary of
    what's left: the mean, the median, or the mode (most frequent value). Median
    is the default on purpose — it sits in the bulk of a skewed column where the
    mean gets dragged into the empty tail. Returns a length-D vector of fills.
    """
    D = X.shape[1]
    fill = np.empty(D)
    for j in range(D):
        col = X[~mask[:, j], j]                       # observed values only
        if col.size == 0:                             # whole column missing
            fill[j] = 0.0
        elif strategy == "mean":
            fill[j] = col.mean()
        elif strategy == "median":
            fill[j] = np.median(col)
        elif strategy == "mode":
            vals, counts = np.unique(col, return_counts=True)
            fill[j] = vals[counts.argmax()]           # first max on ties, like sklearn
        else:
            raise ValueError(f"unknown strategy {strategy!r}")
    return fill


def impute_transform(X, mask, fill):
    """Replace every masked entry with its column's fitted fill value."""
    out = X.copy()
    rows, cols = np.where(mask)
    out[rows, cols] = fill[cols]
    return out
# endregion


# region: was_missing
def missing_indicator(mask, cols_missing_zero):
    """Binary "was this value imputed?" columns for the features that had gaps.

    Imputing throws away one bit of truth: whether the number is real or a
    guess. Appending an indicator per imputed column hands that bit back to the
    model, which can then learn "a missing insulin reading is itself a signal".
    Returns an (N, k) 0/1 array, one column per feature in cols_missing_zero.
    """
    return mask[:, cols_missing_zero].astype(float)
# endregion


# region: outliers
def zscore_outliers(col, thresh=3.0):
    """Flag points more than `thresh` standard deviations from the mean.

    z = (x - mean) / std. The classic rule calls anything past |z| = 3 an
    outlier. It assumes a roughly bell-shaped column, so it under-flags a
    skewed one (the mean and std are already inflated by the very tail you're
    hunting) — which is exactly why the IQR rule below is the safer default.
    """
    mu, sigma = col.mean(), col.std()
    if sigma == 0:
        return np.zeros(col.shape, dtype=bool)
    z = (col - mu) / sigma
    return np.abs(z) > thresh


def iqr_outliers(col, k=1.5):
    """Flag points outside [Q1 - k·IQR, Q3 + k·IQR] — Tukey's rule.

    IQR = Q3 - Q1 is the width of the middle half of the data. Tukey's fences
    sit k=1.5 IQRs beyond the quartiles. Built from percentiles, the rule
    doesn't care how long the tail is, so it keeps working on skewed columns
    where the z-score quietly fails.
    """
    q1, q3 = np.percentile(col, [25, 75])
    iqr = q3 - q1
    lo, hi = q1 - k * iqr, q3 + k * iqr
    return (col < lo) | (col > hi)
# endregion


def load_data(path="../data/diabetes.csv"):
    """The Pima diabetes set: 8 features, a binary Outcome, and hidden zeros.

    Returns (X, y, names, cols_missing_zero): X is (N, 8) float, y is (N,) int
    0/1, names lists the feature columns in order, and cols_missing_zero is the
    boolean vector marking columns where 0 means missing.
    """
    df = pd.read_csv(path)
    names = [c for c in df.columns if c != "Outcome"]
    X = df[names].to_numpy(float)
    y = df["Outcome"].to_numpy(int)
    cols_missing_zero = np.array([n in ZERO_IS_MISSING for n in names])
    return X, y, names, cols_missing_zero

The library version

You wouldn't hand-roll imputation in production, and scikit-learn's SimpleImputer is the thing you'd use instead. It's the same fit/transform object as everything else in the library, and it does exactly what our impute_fit/impute_transform pair does — with one setup wrinkle. SimpleImputer keys off NaN by default, so the disguised zeros have to be converted to real nulls first; it can't know that a zero is impossible any more than our detector could without being told.

def sklearn_impute(Xtr, Xte, strategy):
    """Fit a SimpleImputer on TRAIN, transform both splits.

    strategy is 'mean', 'median', or 'most_frequent' (our 'mode'). Missing
    entries must already be np.nan. Fitting on Xtr alone and transforming Xte
    with those learned fills is the no-leakage pattern — the imputer never sees
    a test row while learning what to fill with.
    """
    imp = SimpleImputer(strategy=strategy, missing_values=np.nan)
    imp.fit(Xtr)
    return imp.transform(Xtr), imp.transform(Xte), imp.statistics_

strategy="mean", "median", and "most_frequent" line up one-to-one with our mean, median, and mode, and the chapter's checkpoint asserts that our from-scratch fills match SimpleImputer's learned statistics and its transformed output to within 1e-9, for all three strategies, before any chart is drawn. If a formula drifts, the run fails loudly instead of shipping a wrong picture. The model that scores each cleaning strategy is held fixed so the only thing that changes is how the missing values were handled — a standardizer feeding logistic regression, fit on train alone:

def classify_accuracy(Xtr, ytr, Xte, yte, seed=0):
    """Test accuracy of one fixed classifier on whatever cleaning it's handed.

    A standardizer plus logistic regression, both frozen across every run and
    fit on the training split only. Any change in this number is the cleaning
    strategy's doing, not the model's — that isolation is the experiment.
    """
    scaler = StandardScaler().fit(Xtr)
    clf = LogisticRegression(max_iter=1000, random_state=seed)
    clf.fit(scaler.transform(Xtr), ytr)
    return float(clf.score(scaler.transform(Xte), yte))

Scratch versus library

The usual scratch-versus-library face-off is a formality here — the checkpoint already proved the two imputers produce identical numbers, so a bar chart of "our median versus sklearn's median" would be two bars of the same height. The comparison worth making is the one this whole chapter is about: what each cleaning strategy does to a downstream model. Same fixed classifier, same train/test split, four runs — the raw dirty zeros left untouched, then mean-imputed, then median-imputed, then median-imputed with the was-missing flag added:

Now read this honestly, because it does not say what a tidy chapter would want it to say. The raw dirty zeros score 0.781 test accuracy — the highest of the four — and the two naive imputations both drop to 0.750, five patients worse. Cleaning the data made the model worse. That's not a mistake in the code; it's the Pima insulin column being closer to MNAR than to random. Whether a patient's insulin was measured at all is itself correlated with their outcome, so the impossible zeros were accidentally smuggling a real signal into the model, and mean and median imputation erased it by making every missing value look like an ordinary one. This is exactly Rubin's warning from 1976, showing up as a number: when the absence is informative, deleting or overwriting it throws away data.

But you can't keep the zeros. A physiologically impossible value poisons every distance, every coefficient, every summary statistic downstream, and it's a lie sitting in your table waiting to mislead the next person who queries it. The fix is the last bar: median-impute so the numbers are honest, and add the was-missing flag so the model can still learn from the absence. That recovers the signal on clean data — 0.771, back within one patient of the dirty baseline — with none of the corruption. Mean and median tie here, which is why I made the mean-versus-median case on the distribution instead of on this accuracy; the two fills land the same downstream number on this split even though median is the demonstrably saner fill. The lesson isn't that any one bar wins. It's that "clean the data" is not a button, and if you impute without asking why the values were missing, you can scrub a real signal right out of your dataset.

Takeaways

Look at your data before you model it. Not a summary, not the first five rows — actually look, column by column, at the ranges and the zeros and the things that can't be true. I've watched pipelines train for hours on tables where half a column was a data-entry placeholder nobody had noticed, and no amount of model tuning fixes a feature that's 48% fiction. Garbage in, garbage out is the oldest line in computing and it still beats every algorithm choice you'll agonize over later. The Pima set is the standard "clean" benchmark and it has 652 disguised missing values in it; assume yours is dirtier than the docs claim, because it is.

The working rules are short. Disguised missingness is everywhere zero is a legal value for a column that can't physically be zero — hunt for it with domain knowledge, not a null check. When you fill, median-impute skewed columns; the animation showed why, mean imputation piles your guesses out in the empty tail where no real data lives. When a lot is missing, or when you suspect the missingness itself carries signal, add a was-missing indicator column so the model learns from the absence instead of being fooled by the fill — that one cheap column is what turned our cleaned data from a downgrade back into a tie with the dirty baseline. And never impute with statistics computed on your test set; fit the fill on train, apply it everywhere, same discipline as scaling. Do this part right and every model in the rest of the book gets better for free. Skip it and the fanciest algorithm in the world just learns your mistakes faster.