Peer-group robust z-scores, isolation forests, and a geographic lens on provider-level Medicaid claims
Author
Paulina Del Mundo Del Fierro, MD, MPH
Published
June 18, 2026
TipAbstract
Medicaid is the largest source of health coverage for low-income Americans, and the federal government releases monthly snapshots of who got paid for what. The interesting question isn’t “who bills the most” (huge hospital systems dominate by default) but “who bills atypically compared to other providers doing the same thing?” That signal is what flags candidates for fraud review or atypical practice patterns worth a second look. This notebook works through 11 GB of provider-level billing data on a laptop using a small streaming-database trick, then computes, for every provider on every billing code, how far their per-patient spending sits from the typical spending in their peer group.
Two key findings come out. First, the spending data is extremely lopsided — payments per patient span four orders of magnitude — so the textbook “z-score” method fails completely. A robustness fix that’s resistant to extreme values (the median absolute deviation) works instead. Second, when you apply that to thousands of providers, raw cutoffs produce huge numbers of false alarms; pairing the robust z-score with a false-discovery-rate correction keeps the alert list to providers actually worth investigating. The notebook also builds a county-level interactive map of where Medicaid dollars flow by category of care. Two caveats are load-bearing: the peer group as currently defined is too broad (“everyone billing this code” mixes specialties and practice sizes), and there’s no comparison against a list of known excluded providers — so the output is best read as “exploratory surfacing for human review,” not validated fraud detection.
1 Why this analysis, and what you’ll learn from it
Medicaid is the largest source of health insurance for low-income Americans, and the federal Department of Health and Human Services releases a public file every month listing how much each healthcare provider was paid for each procedure or product they billed. That file is around 11 gigabytes of comma-separated text — far too large to open in Excel, but well within reach of a laptop if you process it the right way.
The question we want to ask is not “who got paid the most” — that’s always going to be huge hospital systems and dialysis chains by sheer size. The interesting question is “who got paid atypically, given what their peers were paid for the same procedure?” That’s the signal that fraud-and-abuse review teams, network-integrity auditors, and health-services researchers all want to surface from this kind of data.
By the end of this notebook, you should understand:
How to process a dataset much larger than your computer’s memory — the streaming-database approach in prep.py reads 11 GB of text without ever loading it all at once.
Why the textbook z-score fails on real-world spending data, and how the robust z-score (built on median and median absolute deviation instead of mean and standard deviation) fixes it.
What the false discovery rate (FDR) is, and why ranking provider rows by an FDR-adjusted q-value is honest about the fact that you’re running thousands of statistical tests at once.
How an isolation forest — an unsupervised, tree-based anomaly detector — produces a second opinion on the same data using the joint distribution of multiple features, not just one variable at a time.
How to draw inequality conclusions about geographic spending patterns using the Gini coefficient and the Lorenz curve.
What real research-grade caveats look like when you can only see “who billed what” but not “what diagnosis the patient had” or “whether the provider was later sanctioned.”
NoteA note on terminology — the acronyms used in this notebook
Healthcare claims data carries a lot of jargon. Here are the terms we’ll use repeatedly:
Term
What it means
HCPCS
Healthcare Common Procedure Coding System — a code system for billable procedures, drugs, services, and equipment. CPT codes are a subset of HCPCS.
NPI
National Provider Identifier — the unique 10-digit number every US healthcare provider has.
NPPES
National Plan & Provider Enumeration System — the CMS database that links every NPI to a provider’s address, specialty, and credentials.
NUCC taxonomy
A code system maintained by the National Uniform Claim Committee that maps providers to specialties (e.g., “207RC0000X” = Cardiology).
ZCTA
ZIP Code Tabulation Area — the Census Bureau’s polygon approximations of ZIP code boundaries. ZIPs themselves are postal-route definitions, not geographic shapes.
FIPS
Federal Information Processing Standard — codes for US states (2 digits) and counties (5 digits).
MAD
Median Absolute Deviation — for a set of numbers, take each value’s distance from the median, then take the median of those distances. A robust analog of standard deviation.
Robust z-score
\((x - \text{median}) / (1.4826 \times \text{MAD})\). Reads on the same scale as a standard z-score under normality, but doesn’t blow up when a few extreme values are present.
BH-FDR
Benjamini–Hochberg False Discovery Rate — a procedure for adjusting p-values when you’re running many tests at once, so that no more than a fixed fraction of your “significant” flags are false alarms in expectation.
Isolation forest
An unsupervised anomaly-detection algorithm that builds many random binary trees and scores points by how easily they get isolated (anomalies are isolated quickly).
Gini coefficient
A 0-to-1 number measuring how unequally a quantity is distributed across units. 0 = perfect equality, 1 = one unit has everything.
Lorenz curve
The graphical companion to Gini: plots the cumulative share held by the bottom X% of units.
DuckDB
An in-process analytical database — like SQLite but optimised for analytics queries over large datasets. Reads CSVs and Parquets directly without loading them into memory.
Parquet
A column-oriented binary file format, much smaller and faster to scan than CSV for analytic queries.
Concretely, the work has three phases:
Phase 1 — prep step. Stream the 11 GB Medicaid Provider Spending CSV through DuckDB, aggregate to one row per (provider, billing code), join in NPPES (for specialty and state) and the Census ZCTA file (for geographic coordinates), then save the result as committed Parquet files. This phase is automated by prep.py; the notebook itself reads only the smaller, processed Parquets.
Phase 2 — outlier surfacing. For every (provider, billing code) row, compute the robust z-score relative to peer providers billing the same code. Rank the providers by FDR-adjusted q-value. Cross-check the ranking against an isolation forest fit on a joint feature space.
Phase 3 — geographic atlas. Bucket every billing code into one of 12 spending categories using a rule-based classifier, then map county-level paid-per-beneficiary as a 12-panel interactive choropleth. Pair the maps with Gini concentration numbers and state-by-state hotspot tables.
NoteData sources
On the pathway → 01 · Measurement: data provenance and ontology-aware claims, knowing what an HCPCS code does and does not capture.
The raw files stay local on disk (11 GB Medicaid + 45 GB NPPES, far too large to commit). The prep.py script runs one DuckDB pipeline that produces seven committed Parquet artifacts under data/processed/; this notebook reads only those.
2 Setup
Show code
suppressPackageStartupMessages({library(arrow)library(dplyr)library(ggplot2)library(tidyr)library(stringr)library(forcats)library(scales)library(sf)library(leaflet)library(reticulate)})reticulate::py_require(c("pandas", "numpy", "scikit-learn", "matplotlib"))theme_set(theme_minimal(base_size =12) +theme(plot.background =element_rect(fill ="#fbfaf7", colour =NA),panel.background =element_rect(fill ="#fbfaf7", colour =NA),panel.grid.major =element_line(colour ="#e8e5dd"),panel.grid.minor =element_line(colour ="#f3f1ec"),text =element_text(colour ="#14181f"),axis.text =element_text(colour ="#4a5160"),plot.title =element_text(colour ="#14181f", face ="bold"),plot.subtitle =element_text(colour ="#4a5160") ))# Method packages — installed on first render if absent.for (pkg inc("kableExtra")) {if (!requireNamespace(pkg, quietly =TRUE)) {install.packages(pkg, repos ="https://cloud.r-project.org") }}# Lightweight wrapper for hero tables — bootstrap styling, no monospaced# greyscale dump. Plain knitr::kable() still works elsewhere for inline /# diagnostic tables where the styled version is overkill.pretty_kable <-function(x, caption =NULL, ...) { knitr::kable(x, caption = caption, ...) |> kableExtra::kable_styling(bootstrap_options =c("striped", "hover", "condensed", "responsive"),full_width =FALSE,position ="left",font_size =13 )}data_dir <- here::here("traces", "applied", "02-medicaid-outliers", "data", "processed")outliers <- arrow::read_parquet(file.path(data_dir, "outliers.parquet"))hcpcs_sum <- arrow::read_parquet(file.path(data_dir, "hcpcs_summary.parquet"))state_sum <- arrow::read_parquet(file.path(data_dir, "state_summary.parquet"))specialty_sum<- arrow::read_parquet(file.path(data_dir, "specialty_summary.parquet"))taxonomy <- arrow::read_parquet(file.path(data_dir, "taxonomy.parquet"))format_money <-function(x) scales::label_dollar(accuracy =1, scale_cut = scales::cut_short_scale())(x)
Show code
import pandas as pdimport numpy as npimport matplotlib.pyplot as pltfrom sklearn.ensemble import IsolationForestfrom sklearn.preprocessing import StandardScalerplt.rcParams.update({"figure.dpi": 110,"figure.facecolor": "#fbfaf7","axes.facecolor": "#fbfaf7","savefig.facecolor": "#fbfaf7","axes.edgecolor": "#4a5160","axes.labelcolor": "#14181f","axes.titlecolor": "#14181f","text.color": "#14181f","xtick.color": "#4a5160","ytick.color": "#4a5160","grid.color": "#e8e5dd",})
3 Step 1 — what does the dataset actually contain?
Before any analysis, look at the raw shape of the data: how many providers, how many billing codes, how much money in total, how many candidate “outlier” rows the prep step already flagged for closer inspection. Doing this first catches obvious data-quality problems (missing files, truncated reads, accidentally-empty columns) before they propagate through the rest of the analysis.
This is the diagnostic phase. We’re asking four questions before trusting any downstream model:
Are the seven Parquet files the size and shape we expect? If a join failed or a file got truncated in the prep step, we want to know now.
Is the spending distribution actually heavy-tailed? This is the key empirical claim that motivates using the robust z-score instead of a regular z-score. If the data were close to a bell curve, the regular z-score would be fine and the extra machinery would be unnecessary.
How complete is the geocoding join? If most provider addresses didn’t resolve to a county centroid, a county-level map will be hollow even where the data has signal.
Where does the observed robust-z distribution sit relative to the conventional z ≥ 5 flag threshold? (The histogram below shows the empirical distribution of robust-z — actual data, including any real outliers. It is not the “null distribution,” which would be the hypothetical distribution if every provider were exactly typical. The shape and tail of the observed distribution tells us how aggressive our flag threshold is.)
Show code
artifact_summary <-function(name, df) { tibble::tibble(artifact = name,rows =nrow(df),cols =ncol(df),n_complete =sum(stats::complete.cases(df)),pct_complete =100* n_complete /pmax(nrow(df), 1) )}dplyr::bind_rows(artifact_summary("hcpcs_summary", hcpcs_sum),artifact_summary("state_summary", state_sum),artifact_summary("specialty_summary",specialty_sum),artifact_summary("outliers", outliers),artifact_summary("taxonomy", taxonomy)) |> dplyr::mutate(rows =format(rows, big.mark =","),cols =format(cols, big.mark =","),n_complete =format(n_complete, big.mark =","),pct_complete =sprintf("%.1f%%", pct_complete) ) |> knitr::kable(col.names =c("Parquet artifact", "Rows", "Cols","Complete-case rows", "% complete-case"),caption ="Row × column shape and complete-case rate of each prep artifact. Aggregates (hcpcs / state / specialty summary) are dense by construction; the outliers detail table carries some NA in lat/lon (provider ZIP didn't resolve to a Census centroid)." )
Row × column shape and complete-case rate of each prep artifact. Aggregates (hcpcs / state / specialty summary) are dense by construction; the outliers detail table carries some NA in lat/lon (provider ZIP didn’t resolve to a Census centroid).
Parquet artifact
Rows
Cols
Complete-case rows
% complete-case
hcpcs_summary
4,114
6
4,114
100.0%
state_summary
84
7
84
100.0%
specialty_summary
768
9
768
100.0%
outliers
1,089,059
14
1,061,399
97.5%
taxonomy
879
4
639
72.7%
Show code
paid_per_bene_plot <-ggplot( outliers |> dplyr::filter(paid_per_bene >0, !is.na(paid_per_bene)),aes(paid_per_bene)) +geom_histogram(fill ="#1f4f99", color ="white", bins =60) +scale_x_log10(labels = scales::label_dollar()) +labs(x ="Paid per beneficiary (log10)", y ="Rows", subtitle ="Per-row paid/bene")bene_plot <-ggplot( outliers |> dplyr::filter(beneficiaries >0, beneficiaries <=5000),aes(beneficiaries)) +geom_histogram(fill ="#1f4f99", color ="white", bins =60) +labs(x ="Beneficiaries (≤ 5,000)", y ="Rows", subtitle ="Per-row beneficiary count")peer_plot <-ggplot(hcpcs_sum, aes(peer_n)) +geom_histogram(fill ="#1f4f99", color ="white", bins =60) +scale_x_log10(labels = scales::label_comma()) +labs(x ="Peer providers per HCPCS (log10)", y ="HCPCS codes",subtitle ="Peer-group sizes")z_plot <-ggplot( outliers |> dplyr::filter(!is.na(robust_z), robust_z >=-10, robust_z <=50),aes(robust_z)) +geom_histogram(fill ="#1f4f99", color ="white", bins =60) +geom_vline(xintercept =5, linetype ="dashed",color ="#be123c", linewidth =0.6) +labs(x ="Robust z (clipped to [-10, 50])", y ="Rows",subtitle ="Robust-z distribution")# Use patchwork if available, then gridExtra, then a plain stacked print.if (requireNamespace("patchwork", quietly =TRUE)) { patchwork::wrap_plots(paid_per_bene_plot, bene_plot, peer_plot, z_plot, ncol =2)} elseif (requireNamespace("gridExtra", quietly =TRUE)) { gridExtra::grid.arrange(paid_per_bene_plot, bene_plot, peer_plot, z_plot, ncol =2)} else {for (p inlist(paid_per_bene_plot, bene_plot, peer_plot, z_plot)) print(p)}
Three distributional sanity checks. Top-left — paid per beneficiary across all (NPI × HCPCS) rows on a log10 scale; the long upper tail is exactly why MAD-based robust z’s outperform a mean/SD z-score here. Top-right — beneficiary count per row, capped at 5,000 for visibility. Bottom-left — peer-group sizes per HCPCS (the denominators behind every robust z); HCPCS with peer_n < 25 are filtered out in prep. Bottom-right — robust-z values (those with a finite z) showing the long right tail; the orange dashed line marks the program-integrity ≥ 5 cutoff.
Outlier-table rows total : 1,089,059
with NPPES state : 1,076,259 (98.8%)
with NPPES specialty : 1,076,255 (98.8%)
with ZCTA → lat/lon centroid : 1,061,403 (97.5%)
The geo-coverage line is the join-quality check that any geographic claim downstream depends on: a county-level atlas built from a 60% NPPES-state coverage would mean the map is showing 60% of the data, not 100% — and the difference matters for any conclusion about state-level Medicaid policy.
4.1 Overdispersion and zero-inflation diagnostics
Some of the columns in our data are counts — number of beneficiaries served, number of claims billed, number of months active. If we ever wanted to model these counts (predict them as a function of provider characteristics), we’d need to pick a probability distribution to assume. The simplest count model, the Poisson, has the curious property that its variance equals its mean — which is rarely true in real data. Two diagnostics check whether Poisson is even worth considering:
The variance-to-mean ratio (VMR) computes \(\text{Var}(X)/E[X]\). Under Poisson this ratio is 1; a VMR much greater than 1 means the data are overdispersed and a negative binomial (NB) distribution would fit better.
The fraction of zeros tells us whether the count distribution has more zeros than even NB can produce — in which case a zero-inflated NB model would be appropriate.
We’re not actually fitting a count model here (the outlier-flagging is a different task), but documenting these diagnostics is good practice: it tells any future modeller exactly what they’re walking into.
Show code
vmr <-function(x) { x <- x[!is.na(x) & x >=0]if (length(x) ==0||mean(x) ==0) return(c(mean =NA, var =NA, vmr =NA))c(mean =mean(x), var = stats::var(x), vmr = stats::var(x) /mean(x))}zero_pct <-function(x) { x <- x[!is.na(x)]if (!length(x)) return(NA_real_)100*mean(x ==0)}vmr_tbl <-rbind(paid =vmr(outliers$paid),beneficiaries =vmr(outliers$beneficiaries),claims =vmr(outliers$claims),months_active =vmr(outliers$months_active))family_call <-function(v) {if (is.na(v)) return("—")if (v <1.2) "Poisson"elseif (v <5) "Quasi-Poisson / NB"else"Negative binomial / ZINB"}pretty_kable(data.frame(Variable =rownames(vmr_tbl),Mean =sprintf("%.1f", vmr_tbl[, "mean"]),Variance =sprintf("%.1f", vmr_tbl[, "var"]),VMR =sprintf("%.1f", vmr_tbl[, "vmr"]),`% zeros`=sprintf("%.2f%%",vapply(list(outliers$paid, outliers$beneficiaries, outliers$claims, outliers$months_active), zero_pct, numeric(1))),`Family`=vapply(vmr_tbl[, "vmr"], family_call, character(1)),check.names =FALSE ),caption ="Variance-to-mean ratio (VMR), zero-inflation, and the implied regression family for each count-like feature. Caveat: the outliers parquet is already filtered to robust_z >= 2 OR top-20 spend codes AND beneficiaries >= 11 in prep, so the % zeros reads near 0 by construction. The VMR is still informative because it characterizes the conditional dispersion that any future modeling step will need to absorb.")
Variance-to-mean ratio (VMR), zero-inflation, and the implied regression family for each count-like feature. Caveat: the outliers parquet is already filtered to robust_z >= 2 OR top-20 spend codes AND beneficiaries >= 11 in prep, so the % zeros reads near 0 by construction. The VMR is still informative because it characterizes the conditional dispersion that any future modeling step will need to absorb.
Variable
Mean
Variance
VMR
% zeros
Family
paid
paid
610931.6
115566678372174.3
189164684.9
1.62%
Negative binomial / ZINB
beneficiaries
beneficiaries
2740.6
253168488.7
92375.8
0.00%
Negative binomial / ZINB
claims
claims
5845.1
8580957013.7
1468053.1
0.00%
Negative binomial / ZINB
months_active
months_active
26.5
753.9
28.4
0.00%
Negative binomial / ZINB
NoteGoing deeper — why \(\text{Var} = E\) for Poisson, and what NB / ZINB / hurdle each fix
The Poisson distribution with rate \(\lambda\) has the curious property that its variance equals its mean: \(\Pr(X = k) = \lambda^k e^{-\lambda}/k!\), and from the moment-generating function \(E[X] = \lambda\) and \(\text{Var}(X) = \lambda\). So VMR = variance/mean = 1 in expectation; a sample with VMR much above 1 falsifies the Poisson assumption — the data are overdispersed. Negative binomial (NB) is the standard fix, built as a Poisson-Gamma mixture: \(X \mid \mu \sim \text{Poisson}(\mu)\) where \(\mu\) itself varies as \(\text{Gamma}(\theta)\). The marginal has \(\text{Var}(X) = \mu_0 + \mu_0^2/\theta\) — the extra term absorbs unobserved heterogeneity across observations. Zero-inflated NB layers a second mechanism on top: with probability \(\pi\) the observation is a structural zero (forced to \(X = 0\)), and with probability \(1 - \pi\) it’s an ordinary NB draw. Use ZINB when zeros are categorically different from small counts (a doctor who never prescribes opioids vs. one who prescribes them rarely). A hurdle model is similar but cleaner: model “any vs. none” with one distribution, then conditional on positive, model “how much” with a truncated NB — useful when enrollment and intensity are governed by different processes. As a concrete example, 1,000 providers each filing a beneficiary count with \(E[X] = 40\) and \(\text{Var}(X) = 380\) has VMR = 9.5, strongly overdispersed; the NB fit solves \(380 = 40 + 1600/\theta\) for \(\theta \approx 4.7\).
4.2 Within-flagged-set selection bias
The outliers parquet is itself a non-random sample of the full DuckDB output (kept rows are those with robust_z >= 2 OR billed under the top-20 spending HCPCS codes). We can’t directly compare against the dropped rows without re-running the 11 GB DuckDB pipeline, but we can check whether the high-z subset (robust_z >= 5, the program-integrity threshold) systematically differs from the moderate-z subset (2 <= robust_z < 5) on the variables we’d otherwise model — a “what makes the flagged-flagged set different from the merely flagged” check. The seven tests below are descriptive: independence assumptions are mildly violated (some NPIs appear in multiple rows on different HCPCS), and we do not apply a multiplicity correction across the seven, so the p-values should be read as evidence of direction rather than as calibrated significance tests.
knitr::kable( bias_df,caption ="Within-flagged-set bias check: do high-z (>= 5) rows differ from moderate-z (2 - 5) rows on the inputs that would feed any peer-comparison model? The 'high-z is just larger providers' null can be ruled out by checking the beneficiary and claim distributions; specialty mix differences expose which clinical areas drive the high-z tail.")
Within-flagged-set bias check: do high-z (>= 5) rows differ from moderate-z (2 - 5) rows on the inputs that would feed any peer-comparison model? The ‘high-z is just larger providers’ null can be ruled out by checking the beneficiary and claim distributions; specialty mix differences expose which clinical areas drive the high-z tail.
Variable
Test
Statistic
p-value
beneficiaries
beneficiaries
Mann-Whitney U
4.622618e+10
2.2e-60
claims
claims
Mann-Whitney U
4.574995e+10
7.3e-119
paid
paid
Mann-Whitney U
4.271721e+10
<1e-300
paid_per_bene
paid_per_bene
Mann-Whitney U
4.170076e+10
<1e-300
months_active
months_active
Mann-Whitney U
4.640249e+10
6.1e-44
state
state
Pearson chi-squared
1.421810e+04
<1e-300
specialty_class
specialty_class
Pearson chi-squared
3.494460e+04
<1e-300
Show code
hcpcs_totals <- hcpcs_sum |>arrange(desc(total_paid)) |>mutate(rank_pct =row_number() /n(),cum_share =cumsum(total_paid) /sum(total_paid) )ggplot(hcpcs_totals, aes(rank_pct, cum_share)) +geom_line(color ="#1f4f99", linewidth =1) +geom_hline(yintercept =0.8, linetype ="dashed", color ="#be123c") +scale_x_continuous(labels =percent_format()) +scale_y_continuous(labels =percent_format()) +labs(x ="Fraction of HCPCS codes (sorted from top to bottom spender)",y ="Cumulative share of total paid")
Cumulative share of total Medicaid paid dollars, ordered by HCPCS code from highest- to lowest-spending. A steep early curve means a small set of codes drives most public outlay.
On the pathway → 02 · Model: robust statistics for heavy tails and multiplicity control, so a handful of extreme providers do not dominate the result.
This is the analytical heart of the notebook. Our goal is to identify (provider, billing-code) pairs where the payment per patient is much higher than what other providers billing the same code received — controlling for the fact that we’re running this comparison for hundreds of thousands of rows, so some will look extreme just by chance.
The recipe is:
For each billing code (HCPCS), compute the median paid-per-beneficiary across all providers billing that code. This is the peer-group reference value.
For each provider-billing-code row, compute a robust z-score: how many normal-scale “standard deviations” above the peer median is this row, using the median absolute deviation (MAD) as the scale.
Convert each robust z-score into a tail probability using the standard normal distribution as a reference.
Apply the Benjamini–Hochberg correction to control the expected false-discovery proportion across all the rows being tested.
In one sentence. Compute a MAD-based robust z per (NPI × HCPCS) row, then apply Benjamini–Hochberg FDR ≤ 0.01 instead of a hard z-cutoff so the flag list comes with a multiplicity guarantee, not just an arbitrary threshold.
WarningThis is surfacing, not detection — read the flag list as exploratory
Two reasons the framing matters for a portfolio reader:
The peer group is too broad. “Everyone billing this HCPCS” mixes specialties, panel sizes, and care settings into one comparison. A correct peer group is specialty × panel-size band × setting. This is the foundational caveat in the README and it isn’t resolved in the current pipeline — every flagged outlier could be explained by “wrong peer group” before “anomalous billing behaviour.” A publication-grade or deployable version of this analysis would require a re-run with that peer-group adjustment as the first step.
No ground-truth validation. The flag list has not been compared against the HHS-OIG List of Excluded Individuals/Entities (LEIE) or any other ground-truth set of providers known to have been sanctioned. The robust-z + BH-FDR output is therefore a prioritisation tool — “look at these rows first” — not a detector whose precision and recall have been measured.
Read every “flagged” count below in that frame: surfaced for review, not confirmed as anomalous. The methodology is correct; the validation step that turns it into a deployable detector is the missing piece.
NoteWhy a robust (MAD-based) z and not a normal z?
The classical z-score is (x − mean) / SD. It works beautifully when the data are roughly normal. It fails completely on heavy-tailed data because a few extreme outliers blow up the SD, which then deflates every other observation’s z toward zero — the outliers hide themselves by inflating the noise scale.
The robust z replaces the two non-robust pieces:
Median instead of mean (robust to outliers; doesn’t move when you change extreme values)
MAD (median absolute deviation) instead of SD (also robust). Multiplied by 1.4826 so that under a normal distribution, MAD × 1.4826 estimates the SD — making the resulting z-score interpretable on the familiar scale (“z = 5 means about 5 SDs from the median”).
Use robust z whenever your data have heavy tails (claims data, gene expression, financial returns, anything log-skewed). Use the classical z when you’ve genuinely confirmed normality (rare in real-world health data). (Hampel 1974; Rousseeuw & Croux 1993)
NoteGoing deeper — where the constant 1.4826 comes from
The MAD by itself is a small number, not the same scale as a standard deviation; 1.4826 is the rescaling factor that makes MAD estimate \(\sigma\)if your data is normally distributed. The derivation: for \(X \sim N(\mu, \sigma^2)\), the absolute deviation \(|X - \mu|\) has a half-normal distribution. The median of \(|X - \mu|\) satisfies \(\Pr(|X - \mu| \le \text{MAD}) = \tfrac{1}{2}\). Since the standardised \(|Z| = |X - \mu|/\sigma\) has CDF \(2\Phi(z) - 1\), setting this equal to \(\tfrac{1}{2}\) gives \(\Phi(\text{MAD}/\sigma) = \tfrac{3}{4}\), so \(\text{MAD}/\sigma = \Phi^{-1}(0.75) \approx 0.6745\), and \(\sigma \approx \text{MAD}/0.6745 = 1.4826 \cdot \text{MAD}\). Exact for normal data; approximate otherwise. The corresponding “robust z” \(= (x - \text{median})/(1.4826 \cdot \text{MAD})\) reads on the same scale as a regular z-score but is built from quantities with high breakdown point — the fraction of arbitrarily bad observations an estimator tolerates before producing nonsense. The mean has breakdown 0 (one bad point moves it arbitrarily); both median and MAD have breakdown 0.5. Example: for \(\{10, 11, 12, \ldots, 18, 1000\}\), mean = 112.6, SD ≈ 313, so the classical z for 1000 is just 2.83 (barely flagged — the outlier inflated its own scale). The median is 14.5, MAD = 2.5, so robust z = 266 — the outlier shows up as obviously extreme because the scale didn’t bend to accommodate it.
NoteWhy Benjamini–Hochberg FDR — and what does “FDR ≤ 0.01” actually mean?
When you test thousands of hypotheses, some will look significant just by chance. At p < 0.05 across 100,000 tests, you’d expect 5,000 false positives even if no real signal exists. A hard z ≥ 5 cutoff has the same problem: at z = 5 in a normal null, p ≈ 6 × 10⁻⁷, but across 100,000 tests you’d still expect ~0.06 false positives per HCPCS code — and the data have hundreds of HCPCS codes.
The Benjamini–Hochberg (BH) procedure controls the false discovery rate (FDR) — the expected fraction of false positives among the rows you flag. “FDR ≤ 0.01” means: of the rows the rule flags, no more than 1% are expected to be false positives. That’s a much more honest claim than “p < 0.05” on each individual test.
How BH works: sort the p-values from smallest to largest. For the i-th smallest, compare it to (i / m) × α, where m is the total number of tests and α is your target FDR. The largest i where the test passes defines the cutoff; flag everything at or below it. Built into base R as p.adjust(method = "BH").
Use BH whenever you’re testing many hypotheses simultaneously — genomics, neuroimaging, fraud detection, A/B-testing dashboards. Bonferroni (divide α by m) is more conservative; use Bonferroni when you need a strict family-wise error rate guarantee, BH when controlling the false discovery proportion is enough. (Benjamini & Hochberg 1995; Benjamini & Yekutieli 2001)
NoteGoing deeper — why the BH algorithm actually controls FDR
Two facts power BH. First, under the null, p-values are uniformly distributed on \([0,1]\), so among \(m_0\) true nulls, the expected number with \(p \leq t\) is exactly \(m_0 t\). Second, the BH cutoff is set so that if we reject the smallest \(i\) p-values, the expected count of false rejections among them is at most \(\alpha i\). The algorithm: sort p-values \(p_{(1)} \leq p_{(2)} \leq \cdots \leq p_{(m)}\), find the largest \(i^*\) with \(p_{(i^*)} \leq (i^*/m) \cdot \alpha\), and reject everything up to \(i^*\). Equivalently, define the q-value \(q_{(i)} = \min_{j \geq i} m \cdot p_{(j)}/j\) and reject all tests with \(q \leq \alpha\). Under independence (or “positive regression dependency”), Benjamini & Hochberg proved \(\text{FDR} = E[V / \max(R, 1)] \leq (m_0/m) \cdot \alpha \leq \alpha\), where \(V\) is false discoveries and \(R\) is total discoveries — a martingale argument on the ordered p-values. Contrast with Bonferroni: Bonferroni rejects \(p \leq \alpha/m\) for every test, controlling the family-wise error rate (probability of any false rejection); much stricter, much less powerful. As a 10-test example with \(\alpha = 0.05\): p-values \((0.001, 0.008, 0.012, 0.04, 0.05, 0.10, \ldots)\) compared against the cutoffs \((0.005, 0.010, 0.015, 0.020, 0.025, \ldots)\) — the first three pass the linear BH cutoff while the fourth (\(p = 0.04 > 0.020\)) does not. BH rejects 3 hypotheses; Bonferroni at the same \(\alpha\) would reject only the very first (\(p < 0.005\)). BH finds 2 extra real discoveries while keeping the expected false-discovery proportion at 5%.
The DuckDB prep step computed a robust z-score per HCPCS using the median absolute deviation (MAD) scaled by 1.4826. Under normality, a robust z of 5 corresponds to roughly five standard deviations above the peer median — a common heuristic threshold for program-integrity screens. On heavy-tailed spending data the correspondence is approximate (extreme values are intrinsically more common than a normal would predict), so a fixed z-cutoff at any level is a heuristic rather than a calibrated test.
To guard against false alarms from running hundreds of thousands of comparisons, we also compute Benjamini–Hochberg-adjusted q-values and flag at q ≤ 0.01. Two caveats are essential to read these q-values correctly. First, the p-values feeding BH come from a standard-normal reference applied to robust z’s whose true null distribution is heavier-tailed than normal — so the p-values systematically understate tail probabilities, and the strict FDR-control guarantee doesn’t hold at the nominal level. Second, the outliers parquet was pre-filtered upstream to keep rows with robust_z ≥ 2 or membership in the top-20 spend codes; applying BH within this pre-filtered family controls the false-discovery proportion given the filter, not over the full population of providers. The honest reading: the q-values rank rows by extremity in a multiplicity-aware way, and the resulting list is a defensible prioritisation for human review — not a calibrated detector.
Show code
flagged <- outliers |> dplyr::filter(!is.na(robust_z), beneficiaries >=30) |> dplyr::mutate(# One-sided upper-tail p-value from a standard-normal reference,# consistent with our one-sided flag (only over-billing is of# interest). pnorm(-z) = Pr(Z > z) for standard normal Z.# Caveat: the true null distribution of robust z on heavy-tailed# spending data has fatter tails than normal, so this approximation# UNDERSTATES the true tail probability for large z. The reported# q-values are therefore not strictly FDR-controlling at the nominal# level on the global universe of all NPI x HCPCS rows; they are# an empirical prioritisation score, not calibrated p-values.p_value = stats::pnorm(-robust_z),q_value = stats::p.adjust(p_value, method ="BH"),flag_z5 = robust_z >=5,flag_fdr = q_value <=0.01 )flagging_summary <-data.frame(Rule =c("Raw z >= 5","BH-FDR q <= 0.01 (one-sided)","Both rules agree (intersection)"),Flagged =c(sum(flagged$flag_z5),sum(flagged$flag_fdr),sum(flagged$flag_z5 & flagged$flag_fdr)))flagging_summary$`% of tested`<-sprintf("%.2f%%",100* flagging_summary$Flagged /nrow(flagged))pretty_kable( flagging_summary,caption =sprintf("Flagging counts for %s NPI x HCPCS rows tested (after the beneficiaries >= 30 minimum-sample filter and the upstream prep-step filter of robust_z >= 2 OR top-20 spend codes). The raw z >= 5 cutoff is a common program-integrity heuristic. The BH-FDR rule controls the expected false-discovery proportion *within this pre-filtered family*; because the family was selected on the test statistic, the q-values are not globally calibrated, and we treat them as a prioritisation score rather than calibrated probabilities. The intersection is the highest-confidence subset for downstream review.",format(nrow(flagged), big.mark =",") ))
Flagging counts for 907,705 NPI x HCPCS rows tested (after the beneficiaries >= 30 minimum-sample filter and the upstream prep-step filter of robust_z >= 2 OR top-20 spend codes). The raw z >= 5 cutoff is a common program-integrity heuristic. The BH-FDR rule controls the expected false-discovery proportion *within this pre-filtered family*; because the family was selected on the test statistic, the q-values are not globally calibrated, and we treat them as a prioritisation score rather than calibrated probabilities. The intersection is the highest-confidence subset for downstream review.
Rule
Flagged
% of tested
Raw z >= 5
214160
23.59%
BH-FDR q <= 0.01 (one-sided)
391206
43.10%
Both rules agree (intersection)
214160
23.59%
Show code
top_outliers <- flagged |>filter(flag_fdr) |>arrange(q_value, desc(robust_z)) |>slice_head(n =25) |>transmute( npi, state, specialty_name, hcpcs, beneficiaries, paid,paid_per_bene = paid / beneficiaries, hcpcs_median_ppb,robust_z =round(robust_z, 1),q_value =format.pval(q_value, eps =1e-300, digits =2) )pretty_kable( top_outliers |>mutate(paid =format_money(paid),paid_per_bene =format_money(paid_per_bene),hcpcs_median_ppb =format_money(hcpcs_median_ppb) ),caption ="Top 25 peer-group outliers ranked by BH-FDR q-value (FDR ≤ 0.01, ≥ 30 beneficiaries). The q-value column is the false-discovery-adjusted p — these are the rows the multiplicity-controlled rule keeps.")
Top 25 peer-group outliers ranked by BH-FDR q-value (FDR ≤ 0.01, ≥ 30 beneficiaries). The q-value column is the false-discovery-adjusted p — these are the rows the multiplicity-controlled rule keeps.
Per-beneficiary spend by specialty (top 15 by flagged-outlier count). Box spans p25–p75 with median; whiskers are p05 / p95 from the full peer distribution computed in the DuckDB prep step. Orange dots are the individual robust-z ≥ 5 outliers layered on top.
6 Step 4 — a second opinion: unsupervised anomaly detection
On the pathway → ∗ · Defend it: a second opinion that tries to corroborate the flags a different way, rather than confirm them.
NoteWhy an isolation forest as a second opinion?
The robust z asks one question: is this row’s paid-per-beneficiary unusually far from the peer median? It looks at one variable at a time. Isolation forest asks a richer question: is this row unusual considering the joint distribution of all features at once? A row could look normal on every individual feature but still be an outlier in the joint space (high beneficiaries + low claims + high paid is a different story from high beneficiaries + high claims + high paid).
How it works: the algorithm builds many random binary trees by splitting on random features at random thresholds. Outliers are easier to isolate — they get separated into their own leaf in fewer splits. The “anomaly score” for a row is roughly the average path length from root to its leaf across the trees; shorter path = more anomalous. No labels needed.
Use unsupervised methods (isolation forest, local outlier factor, one-class SVM, autoencoders) when you don’t have a labeled training set of “fraud” vs “not fraud” — which is almost always the case in healthcare. Treat the agreement between a rule-based method (robust z) and an unsupervised method as your highest-confidence flag set, and the disagreement as the interesting set worth a human review. (Liu, Ting & Zhou 2008, 2012)
NoteGoing deeper — isolation forest scoring and the \(c(n)\) normalisation
Each isolation tree is built by recursively splitting a random subsample of \(\psi\) observations on a random feature at a random cutoff. The split continues until every point is alone in its own leaf. “Normal” points require many splits to isolate (they have many neighbours, each requiring a cut); outliers — sitting in low-density regions — get isolated quickly with few splits. The path length\(h(x)\) from root to leaf is the diagnostic. To make path lengths comparable across sample sizes, divide by the expected path length of a binary search tree with \(\psi\) nodes: \[c(\psi) = 2\,H(\psi - 1) - \tfrac{2(\psi - 1)}{\psi}, \quad H(i) = \ln(i) + 0.5772 \text{ (Euler–Mascheroni)}.\] The anomaly score \(s(x, \psi) = 2^{-E[h(x)]/c(\psi)}\) then maps the unbounded path length onto \([0, 1]\) — values near 1 indicate strong anomalies (isolated near the root), values near 0 indicate likely inliers, and ~0.5 means average-depth (no signal). The contamination=0.01 parameter isn’t a property of the data; it’s an alert budget — flag the top 1% of scores for human review. Set higher when reviews are cheap, lower when false positives are expensive. Concrete numbers: with the default \(\psi = 256\), \(c(256) \approx 10.25\); a row with average path length 4 across 100 trees scores \(2^{-4/10.25} \approx 0.76\) (strong anomaly), while a row with path length 12 scores \(\approx 0.45\) (inlier-ish).
This step adds a second opinion: an isolation forest, an unsupervised tree-based anomaly detector that scores each row on how unusual it looks in the joint feature space. We fit it on the top 20 billing codes by total spend (where peer groups are large enough to be informative), standardise the features so they’re on comparable scales, and then look at the overlap between rows flagged by the isolation forest and rows flagged by the robust-z rule.
The two methods are doing different things: the robust-z is a deterministic, auditable per-code rule and the isolation forest is an unsupervised, joint-feature, probabilistic model. Where they agree, we have higher confidence; where they disagree, there’s usually something interesting going on (a row that’s extreme on one variable but typical on the others, for example).
6.1 Feature collinearity: a brief check before fitting
When two features are nearly perfect substitutes for each other (e.g., total claims and total beneficiaries — each claim usually covers about one beneficiary, so the two numbers move in lock-step), we say they are collinear. For a linear model, collinearity inflates the standard errors of the affected coefficients — the model can’t tell which feature is “doing the work.” Tree-based models like the isolation forest don’t have this problem in the same way: trees just pick whichever feature to split on, and substitutability doesn’t break inference. But heavy collinearity still has a subtler cost: it reduces tree diversity in the ensemble, because many trees will pick essentially the same split. Three diagnostics give a quick read on how bad the collinearity is:
Pearson correlation matrix — pairwise linear association, \(r \in [-1, 1]\). Values near \(\pm 1\) flag perfectly substitutable variables.
Condition number \(\kappa\) — the ratio of the largest to smallest singular value of the standardised design matrix. Conventional thresholds: \(\kappa < 10\) clean, \(10\)–\(30\) moderate, \(> 30\) problematic.
Variance inflation factor (VIF) per feature — for each feature, regress it on all the others; VIF \(= 1/(1-R^2)\). The same idea as the condition number, but localised to each feature.
What the numbers mean here: claims and beneficiaries are essentially the same construct measured two ways (each claim covers approximately one beneficiary), so r ≈ 0.95+ and a VIF in the double digits on both is the expected reading — not a bug. paid_per_bene is by construction orthogonal to beneficiaries (it’s paid divided by beneficiaries) and months_active carries an independent activity-duration axis. We keep all four features anyway: the isolation forest splits don’t over-count the redundant pair the way a linear regression would, and dropping claims would lose the small but real per-beneficiary multi-claim variation in the long tail. The condition-number print just makes that decision auditable rather than implicit.
NoteGoing deeper — singular values and the condition number \(\kappa\)
A design matrix \(\mathbf{X}\) stretches different directions of input space by different amounts. The singular value decomposition\(\mathbf{X} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\top\) writes those stretch factors as the singular values \(\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_p \geq 0\) on the diagonal of \(\mathbf{\Sigma}\). The condition number is the ratio of the most-stretched direction to the least: \[\kappa(\mathbf{X}) = \sigma_1 / \sigma_p.\] A large \(\kappa\) means the matrix is nearly singular — small changes in \(y\) produce huge changes in \(\hat\beta\), because there’s a near-zero-stretch direction the predictions barely depend on. Concretely, the relative error in \(\hat\beta\) is bounded by \(\kappa\) times the relative error in \(y\). Belsley, Kuh & Welsch’s conventional thresholds: \(\kappa < 10\) negligible multicollinearity, \(10\)–\(30\) moderate (coefficients interpretable but inflated SEs), \(\geq 30\) severe (coefficients unstable under small data perturbations). VIF is the per-feature counterpart of the same diagnostic — global vs. per-coefficient. The relationship: \(\kappa^2 \geq \max_j \text{VIF}_j\) for a standardised design, so a high VIF implies high \(\kappa\) but \(\kappa\) can also be inflated by interaction patterns VIF misses. Example: three standardised features whose covariance eigenvalues are 2.8, 0.95, 0.05 produce \(\kappa = \sqrt{2.8/0.05} \approx 7.5\) — moderate, no action. If the smallest eigenvalue were 0.01 instead, \(\kappa\) jumps to 17, and dropping a feature or applying ridge regression (which adds a small constant to the diagonal to bound \(\kappa\)) becomes appropriate.
6.2 Fitting the isolation forest
Show code
X = StandardScaler().fit_transform(X_raw.to_numpy())iso = IsolationForest( n_estimators=400, contamination=0.01, random_state=42, n_jobs=-1)sub["iforest_score"] =-iso.fit(X).score_samples(X) # higher = more anomaloussub["iforest_flag"] = iso.predict(X) ==-1print(f"Rows analysed: {len(sub):,}")
Paid-per-beneficiary vs. beneficiary count for the top 20 HCPCS codes by spend. Isolation-forest flags in orange overlap substantially with — but not identically to — the robust-z rule.
7 Step 5 — where does the dollar go? A geographic atlas
On the pathway → 03 · Estimate: reading the result on the absolute and spatial scale where decisions actually live.
In one sentence. A county-level interactive map across all 12 spending categories, shaded by where each county sits within the distribution for its category — colours are comparable inside a category but not across them.
Up to this point the analysis has answered “who bills atypically?” — useful for a Program Integrity team. A different audience (policy researchers, state Medicaid administrators, journalists) wants the inverted question: “which kinds of care cost the most, and where geographically does the money flow?” Answering that requires three additional pieces of data engineering that the prep.py script already produced:
A HCPCS → category classifier that buckets every billing code into one of 12 spending categories (Behavioral Health, Dental, Dialysis & ESRD, Prescription Drugs, and so on). The classifier is a deterministic rule-based map; the code lives below for transparency.
A Census 2020 ZIP → county crosswalk that translates every provider’s practice ZIP into the standard 5-digit county FIPS code, allowing county-level aggregation.
A county × category aggregate with within-category ranking so that the colour scale reads consistently across the 12 panels: darker shades always mean “high spend relative to other counties in this same category.”
Two important interpretive caveats up front: (1) The 12 categories are mostly service-type groupings that the HCPCS code system itself organises by — not disease categories. HCPCS codes don’t carry diagnoses. Reading the atlas as “where Medicaid spends on Behavioral Health services” is supportable; reading it as “where Medicaid spends on depression” is not. (2) The within-category colour scale is relative — the darkest navy in the Dental panel and the darkest navy in the Behavioral Health panel both mean “this county is among the highest in its category,” but they are not directly comparable dollar levels across categories.
NoteWhy county, not HSA or ZIP3?
For a Medicaid cost atlas, counties are the defensible middle path — clean TIGER polygons, HUD/Census maintained crosswalks, and the unit most federal-agency publications (MACPAC, KFF, Urban Institute) use. HSAs (Dartmouth Atlas) are the gold standard for Medicare market analysis but were built from Medicare patient flows, so they’re slightly miscalibrated for Medicaid. ZIP3 is tempting because it’s one LEFT(zip, 3) away, but ZIP3 boundaries are USPS logistics geography — a high Modifiable-Areal-Unit-Problem exposure for a spatial-stats story.
Show code
county_atlas <- arrow::read_parquet(file.path(data_dir, "county_group_atlas.parquet"))group_totals <- arrow::read_parquet(file.path(data_dir, "group_totals.parquet"))top_categories <- group_totals |>arrange(desc(total_paid)) |>pull(disease_group)# Recompute a within-category decile rank on paid-per-beneficiary. The parquet# carries quintiles from prep.py; we override with 10 bins here so shading# shows more granular within-category variation without re-running the# DuckDB prep. For categories with < 100 billing counties some deciles will# be tied — the lambda collapses those cleanly.county_atlas <- county_atlas |>group_by(disease_group) |>mutate(decile_ppb = dplyr::ntile(paid_per_bene, 10)) |>ungroup()
Total national Medicaid paid by HCPCS-derived category, 2024. Ranked descending by total paid. ‘Other Medical Services’ captures 99000-series codes the classifier doesn’t yet disambiguate.
