Holes, caps and spikes: missing values, censoring and outliers

Munging, and why a number needs a meaning

Lesson 35 followed the wrangling steps of the slides. Before the statistics, the deck adds one more word (p. 1102):

Both properties are warnings. Lossy: strip the punctuation and "U.S." and "US" become one string, with no way back; strip the HTML tags and the table that said which number was the price is gone. Vague: you munge while still finding out what the data needs, so nobody writes the steps down. Hence the habit: keep the raw file, munge a copy, record each step as a replayable script (the tools are in lesson 37).

A column without a name

The next slide is a single column of thirteen numbers under the header "???", with a large question mark beside it (p. 1103). Look at the left panel of Fig. 1 and decide which number is odd. The second build (p. 1104) reveals the header, "Age(Years)", draws an R boxplot and adds one sentence: "The value stands in the abnormal".

R's boxplot decides with a rule the slide does not print (background): a value more than 1.5 interquartile ranges beyond the box is drawn as a separate point. For the thirteen ages:

The whiskers stop at 45 and 90, the last values inside the fences; 110 lies beyond 102.5 and is drawn alone. That is arithmetic, and it holds under any header. The header adds the verdict. As an age, 110 sits at the edge of human life: a typo for 10 or 101, a misfiled record, or a rare real person. As a body weight in kilograms it would be an ordinary heavy adult. Hence the order of the two builds: first the meaning, then the expected, then the departure.

No expected value without the meaning of the column, and no outlier without an expected value.

Holes you can see, and holes that hide

The deck looks for missing data in two layers (pp. 1105–1106): holes you can see if you look, and damage that leaves the file looking complete.

Mind the last overt check: aggregates can be misleading. A monthly total of calls can be right while one week is missing and another was loaded twice; the errors cancel in the total and show only per week, region or source.

Hidden damage can pass every count. If a loader writes 0 for an absent value, a customer who made no calls and one whose minutes were lost look the same, and only metadata or domain knowledge can tell them apart. The cure is upstream: store a real empty value, never a default that looks like data.

Three ways a value goes missing

"Missing Values: Random" (p. 1107) is a bar chart of visitors per day: days 11, 12 and 16 have no bar, and beside it stand "System failures" and "Complete miss".

A failed collector removes whole records, and nothing about those days caused it: the holes fall at random. Leave them as holes. A bar of 0 would claim that nobody visited; an empty slot says that nobody knows. The mean per recorded day stays fair, while the monthly total is three days short.

"Missing Values: Wrong Ingestion" (p. 1108) is a "CSV to table / excel" accident. Three raw lines go into a six-column table:

Antony Mc James IV, 123 Untidy Cir, Suite# 234, Cumbersome, TX, 76849
Miller Johnson, 1102 Messy Data St, Middletow, OH
Betty Flyier 6483 Phew Lane, Apt A4, FixMe, CA, 91103

The loader splits at the commas and fills Customer, Address Line 1, Address Line 2, City, State and Zip from left to right. Line 1 has six values and lands correctly. Line 2 has four; line 3 has five, because the comma after the name is missing. The slide calls the defects "Missing fields" and "Merged fields" (Fig. 3).

Every value was in the file: the loader created the holes and put wrong values beside them ("OH" is now a city). A check that only counts empty cells may pass the file. The fix belongs in the parser — quote the fields, compare each line's number of values with the header — not in a statistical estimate.

"Missing Values: Inapplicability" (p. 1109) circles two blocks of a medical form, the women's health history and the men's (a PSA blood test for the prostate); every patient leaves one of them empty. The slide's two lines: "Partial data by nature", "Remember to leave empty slots".

Here an empty cell is a statement, not a hole. Filling it invents a PSA test for Anna; a default such as 1900-01-01 recreates the value-or-default confusion; dropping incomplete records deletes every patient. The schema must let the slot stay empty.

Three kinds of holes, three cures: estimate the random ones, fix the ingestion ones in the loader, leave the inapplicable ones empty.

Drop the record, or fill the hole?

Why not simply ignore every record with a hole? The slide "Imputing Values to Missing Data" (p. 1110) answers:

• "In federated data, between 30%-70% of the data points will have at least one missing attribute - data wastage if we ignore all records with a missing value"
• "Remaining data is seriously biased"
• "Lack of confidence in results"
• "Understanding pattern of missing data unearths data integrity issues"

Federated data comes from several sources, each contributing some of the attributes. Why do so many records have a hole when each attribute is missing only now and then? Because the chances multiply:

1. A record has \(p\) attributes. Each is missing with probability \(q\), so each is present with probability \(1 - q\).

2. The record is complete only if every attribute is present. If the holes fall independently, the probabilities multiply:

   \[ P(\text{complete}) = (1 - q)^p \]

3. Ten attributes with 5 % holes each:

   \[ 0.95^{10} = 0.599 \]

   40 % of the records have a hole although 95 % of every column is there. With twenty attributes, \(0.95^{20} = 0.36\): 64 % are lost. Both land inside the slide's 30–70 %.

Wastage is the visible cost; bias is the hidden one. The survivors are a fair sample only if the holes fall purely at random. If heavy callers come from a sloppier source and lose fields more often, the complete records are mostly light callers, and their average bill is too low. The widget shows both, and the two simplest fills.

The last line of p. 1110 turns the problem around: the pattern of the holes is itself a finding. Holes from one source only mean a broken feed, holes on one date a failed load; and the pattern tells you what not to fill at all: ingestion holes and inapplicable slots.

The imputation ladder

To impute is to fill a hole with an estimate. The deck climbs a ladder; each rung uses more of the data and costs more.

Rung 1 is "Standalone imputation" (p. 1111): fill the hole with "Mean, median, other point estimates" of its own column. The slide adds: "Assume: Distribution of the missing values is the same as the non-missing values."; "Does not take into account inter-relationships"; "Introduces bias"; "Convenient, easy to implement". The widget showed the bias: filled values sit at the centre, so the spread shrinks and correlations weaken.

Rung 2, "Better imputation - use attribute relationships" (p. 1112), needs a condition on the pattern of the holes: "Assume : all prior attributes are populated", "That is, monotonicity in missing values." The slide's table is its picture (Fig. 4).

Monotonicity lets you fill the table column by column, left to right, always predicting from columns already full. The slide names "Two techniques": "Regression (parametric)" and "Propensity score (nonparametric)". The first is the "Regression method" (p. 1113): "Use linear regression, sweep left-to-right",

"and so on", where "X3 in the second equation is estimated from the first equation if it is missing". The slide gives no numbers, so fit the first equation on its own table. It is a toy: three complete rows for three unknowns.

1. Rows 1–3 have \(X_3\); row 4 does not. Write the first equation for the three complete rows:

   \[ 3.5 = a + 20b + 1.0c \] \[ 4.0 = a + 18b + 1.1c \] \[ 2.2 = a + 22b + 1.9c \]

2. Subtract the first line from the other two, and \(a\) cancels:

   \[ 0.5 = -2b + 0.1c \] \[ -1.3 = 2b + 0.9c \]

3. Add these two, and \(b\) cancels: \(c = -0.8\). Back-substitute: \(b = -0.29\), then \(a = 3.5 + 5.8 + 0.8 = 10.1\).

4. Fill the hole of row 4, where \(X_1 = 0.9\) and \(X_2 = 15\):

   \[ \hat X_3 = 10.1 - 0.29 \cdot 15 - 0.8 \cdot 0.9 \] \[ \hat X_3 = 5.03 \]

5. Compare with rung 1: the mean of the observed \(X_3\) is \(9.7 / 3 = 3.23\), the median 3.5. Three methods, three answers; and 5.03 is an extrapolation, since \(X_2 = 15\) lies below the observed 18 to 22.

The second technique, "Propensity Scores (nonparametric)" (p. 1114), models the chance that \(X_j\) is missing instead of its value:

1. "Let Y_j=1 if X_j is missing, 0 otherwise"
2. "Estimate P(Y_j =1) based on X₁ through X_(j-1) using logistic regression"
3. "Group by propensity score P(Y_j =1)"
4. "Within each group, estimate missing X_js from known X_js using approximate Bayesian bootstrap."
5. "Repeat until all attributes are populated."

Records with about the same score, the same chance of a hole, are treated as comparable, and each hole gets a value resampled from the observed values of its group (the approximate Bayesian bootstrap). Hence "nonparametric": no formula for \(X_j\) is assumed, although the score comes from a parametric logistic regression.

Rung 3 handles an "Arbitrary missing pattern" with "Markov Chain Monte Carlo (MCMC)" (p. 1115). It assumes that the "data is multivariate Normal, with parameter Θ" (means and covariances) and alternates "(1) Simulate missing X, given Θ estimated from observed X ; (2) Re-compute Θ using filled in X", "Repeat until stable." It is "Expensive: Used most often to induce monotonicity": MCMC fills just enough cells to make a staircase, and the cheap sweep does the rest. The slide ends, in bold:

"Note that imputed values are useful in aggregates but can’t be trusted individually"

The widget measured it (task 4): regression filling missed the mean bill by 0.10 €, a single filled bill by 6.33 € on average. Errors cancel in an average; one guess stays a guess.

Caps and cuts: censored and truncated data

Some damage is done by a bound. The slide "Censoring and Truncation" (p. 1116) calls it "Well studied in Biostatistics, relevant to time dependent data e.g. duration" and separates two cases.

A censored call of 47 minutes is in the table as 20: you know only "at least 20". A truncated customer who called for one minute is not in the table at all, and nothing in it says that anyone is missing. Fig. 5 shows why p. 1106 says "check for spikes and dips".

The two bias the mean in opposite directions (background): capping long calls lowers the average duration, dropping light users makes the average customer heavier. Mnemonic: censored is capped, truncated is thrown out.

Intervals the study only half sees

Durations are where censoring is at home. The figure "Censored time intervals" (p. 1117) draws a study window between "Begin Study" and "End Study" and six intervals against it (Fig. 6).

An interval that began before the window is left censored: its start is unknown. One still running at the end is right censored: you know only that it lasted at least this long. Only "Complete" has a precise duration. Intervals 1 and 2 make the point: identical in the data, very different in reality. Read both as ending at End Study and both durations come out too short, the second by far. Biostatistics meets this daily: a patient still alive when a study ends has a right-censored survival time.

The bias can be mitigated on one condition, stated on p. 1118: "If censoring/truncation mechanism not known, analysis can be inaccurate and biased." "But if you know the mechanism, you can mitigate the bias from the analysis." Hence "Metadata should record the existence as well as the nature of censoring/truncation": capped at 20 minutes, say, or customers under 2 minutes excluded.

Spikes

The last picture (p. 1119) is a plot titled "Potential Data Problem": a noisy "Count" over "Current Accounts" with one spike above 600 and a smaller one above 200, and the caption "Spikes usually indicate censored time intervals caused by resetting of timestamps to defaults". When a system loses an account's real start date and writes a default instead (say, the date of a data migration), you know only that the account existed by then — its interval is censored — and every such account piles onto one value: defaults standing in for missing values.

Censored: the record stays and its value is capped. Truncated: the record is dropped. Both can show in a histogram, as a spike or a dip, and their bias can be mitigated only if metadata records them.

Departure from the expected: outliers

Back to the hook. The slide "Suspicious Data" (p. 1120) lists the eleven numbers and says "“92” is suspicious - an outlier". Then it refuses to convict: "Outliers are potentially legitimate"; "Often, they are data or model glitches"; "Or, they could be a data miner’s dream, e.g. highly profitable customers".

Three verdicts, all possible: a real but rare value; a glitch of the data (a typo for 9, a unit error, a default) or of the model (a wrong "expected" makes a normal value look odd); or your most valuable record, the big spender, the fraud worth finding. Deleting outliers automatically throws the third kind out with the second.

Outliers also move aggregates. Without the 92 the ten numbers have mean 5.6 and median 5.5; with it, mean 13.45 and median 6. One record in eleven more than doubled the mean; the median barely noticed.

Two yardsticks

The classic test divides the distance from the mean by the standard deviation. The robust test (background, not on the slide) uses ingredients an outlier cannot drag: the median and the median absolute deviation, MAD.

For the slide's numbers, \(\bar x = 13.45\) and \(s = 26.13\), so the 92 has \(z = 3.006\). The median is 6, the distances from it are 3, 2, 1, 2, 2, 3, 3, 1, 1, 0 and 86, their median is 2, so \(z_{\text{rob}} = 86 / 2.965 = 29.0\). One yardstick says "barely", the other "very far".

Why the classic score cannot pass 3.015

The masking is built into the formula (background, not on the slide): the outlier sits inside the mean and the standard deviation that judge it. Follow the most extreme case through:

1. Put \(n - 1\) values equal to \(c\) and one at \(c + d\). The mean moves only a little:

   \[ \bar x = c + \frac{d}{n} \]

2. The odd value lies \(d\,(n-1)/n\) above the mean, each of the others \(d/n\) below it. Add up the squared deviations:

   \[ (n-1)\,\frac{d^2}{n^2} + \frac{(n-1)^2\, d^2}{n^2} \] \[ = \frac{(n-1)\, d^2}{n} \]

3. Divide by \(n - 1\) for the sample variance: \(s^2 = d^2 / n\), so \(s = d / \sqrt{n}\). The yardstick grows in proportion to the outlier's own distance.

4. Divide the distance by the yardstick, and \(d\) cancels:

   \[ z = \frac{d\,(n-1)/n}{d / \sqrt{n}} = \frac{n-1}{\sqrt{n}} \]

   For \(n = 11\), \(10 / \sqrt{11} = 3.015\). Spreading the other values apart only lowers the score, so no value in a sample of eleven can score more.

The puzzle, resolved

The hook asked for the z-score of 920 and of 92 million among the ten small numbers. Both give 3.015, the ceiling \((n-1)/\sqrt{n}\) for eleven values: as the outlier moves away, it stretches the standard deviation as fast as its own distance grows. With eleven values a 3-sigma rule can flag a value only by a hair; among the slide's numbers it misses anything below about 72, and it misses both copies when the 92 comes twice. The robust yardstick does not move with the outlier: the 92 stands 29 robust units from the median, the 920 more than 300. The boxplot rule of Fig. 1 catches it too: for the eleven numbers the upper fence is \(7.5 + 1.5 \cdot 3.5 = 12.75\).

Many ways to define "expected"

The next slide (p. 1121) turns the example into a definition: "Outlier – “departure from the expected”", and "Types of outliers – defining “expected”". Every detection method states what you expected; the slide lists "Many approaches".

The geometric kind is familiar from clustering: the noise points of DBSCAN are points that no dense neighbourhood claims.

An outlier is a departure from the expected; every method is a way of defining "expected", and the verdict — legitimate, glitch or a data miner's dream — needs the meaning of the data.

Control charts: the expected band, for one variable and for many

Control charts come from "Quality control of production lots" (p. 1122): every hour a factory measures a small sample of its output and asks whether the process has changed. They are "Typically univariate: X-Bar, R, CUSUM": an X-bar chart plots the mean of each sample, an R chart its range (largest minus smallest value), a CUSUM chart the cumulative sum of the deviations from the target, which catches a small persistent shift. The "Main steps (based on statistical inference)":

1. "Define “expected” and “departure” e.g. Mean and standard error based on sampling distribution of sample mean (aggregate);"
2. "Compute aggregate each sample"
3. "Plot aggregates vs expected and error bounds"
4. "“Out of Control” if aggregates fall outside bounds"

For data quality a sample is a batch: one day's file, one source, one load (Fig. 7).

The bounds are usually three standard errors (background), \(\bar{\bar x} \pm 3\,\sigma / \sqrt{n}\) for samples of \(n\) items: standard errors, because the plotted point is an aggregate. This also explains "Distributional assumptions for charts not based on means e.g. R–charts": by the central limit theorem sample means are close to normal whatever the data, so an X-bar chart needs few assumptions; ranges enjoy no such theorem, so R-chart limits assume a normal population.

The kilogram example

"An Example" (p. 1123) is a real chart from the NIST handbook, "Shewhart control chart for kilogram calibrations": the "Corrections (mg)" of a reference kilogram calibrated again and again from the mid-1970s to the late 1980s. It prints "UCL = mean + 3s = -19.39 mg" and "LCL = mean - 3s = -19.57 mg", and the rest follows:

The limits are \(6s\) apart, hence the 6. Each point is a single calibration, not the mean of a sample, so its standard error is \(s\) itself and the bounds are \(\pm 3s\). A few points around 1987 reach the upper limit: those measurements get investigated. Under normality a point crosses a 3s limit by chance 0.27 % of the time, once in about 370 (background).

Two variables: the rectangle and the ellipse

Watch two correlated variables with one chart each and you get the rectangle of Fig. 8: a point is "in control" when it is inside both pairs of limits. "Multivariate Control Charts - 1" (p. 1124) says why that is not enough: bivariate charts are "based on bivariate Normal assumptions", and "component-wise limits lead to Type I, II errors", because the region of expected points is an ellipse tilted along the correlation (p. 1125).

Read the disagreements as errors of the rectangle, with "in control" as the null hypothesis. A point inside the rectangle but outside the ellipse is fine for each variable alone but implausible for the pair, such as a high X with a low Y when the two always rise together: the rectangle misses the alarm, a Type II error. A point at a tip of the ellipse, beyond one limit, is unusual for one variable but exactly what the correlation predicts: the rectangle raises a false alarm, a Type I error.

The ellipse is the set of points at the same Mahalanobis distance from the centre, the distance that measures in units of the covariance:

Hence the second idea of p. 1124, "Depth based control charts (nonparametric)": "map n-dimensional data to one dimension using depth e.g. Mahalanobis", "Build control charts for depth", "Compare against benchmark using depth e.g. Q-Q plots of depth of each data set". However many attributes a record has, its depth (how far inside the cloud it lies) is one number for an ordinary chart; if new data resemble the benchmark, the Q-Q plot of their depths follows the diagonal.

For signals, p. 1126 adds "Multiscale process control with wavelets": it "Detects abnormalities at multiple scales as large wavelet coefficients", is "Useful for data with heteroscedasticity" (a variance that changes over time) and was "Applied in chemical process control".

Models and baselines

"Model Fitting and Outliers" (p. 1127): "Models summarize general trends in data", "more complex than simple aggregates", like regressions that "focus on attribute relationships"; "Data points that do not conform to well fitting models are potential outliers" — potential, since the model may be wrong. Hence "Goodness of fit tests (DQ for analysis/mining)": "check suitableness of model to data", "verify validity of assumptions", "data rich enough to answer analysis/business question?". The residual analysis of regression, turned on the data.

"Set Comparison and Outlier Detection" (p. 1128) works on whole data sets: the "“Model” consists of partition based summaries" (counts, quantiles, histograms per day or region), and you "Perform nonparametric statistical tests for a rapid section-wise comparison of two or more massive data sets"; "If there exists a baseline “good’’ data set, this technique can detect potentially corrupt sections in the test data set". A control chart for files, rapid because summaries are small.

Every method of this lesson, from the boxplot of the ages to the baseline comparison, does the same two things: it states what is expected, and it measures the departure.