meds_summary_stats

Fast, privacy-safe summary statistics over MEDS datasets, for validating ETLs and catching regressions in them over time.

You ship a new version of an ETL. The dataset version has not changed, so the subject count should not have moved, the schema should be identical, and the code vocabulary should be the same set of strings. This tool turns that expectation into a file you can commit and a check you can run in CI.

[!WARNING] Under active initial development. The stats surface and comparison semantics are not yet stable; pin a version.

Two commands

# Reduce a dataset to a small, aggregated, non-sensitive fingerprint.
extract_MEDS_summary_stats /data/mimic-iv-meds -o baseline.json

# ... later, after an ETL change ...
extract_MEDS_summary_stats /data/mimic-iv-meds -o current.json

# Diff the two against a tolerance policy. Exit code 1 if it regressed.
compare_MEDS_summary_stats baseline.json current.json

Both are also python -m meds_summary_stats extract / ... compare, and both are importable (extract_summary_stats, compare_summary_stats).

A worked example

Three complete MEDS datasets are committed under examples/datasets/, so everything below runs against a fresh clone with no setup. Every output on this page is executed as a doctest, so it is what the tool actually prints.

>>> print_directory(EXAMPLES / "datasets" / "v1_baseline")
├── data
│   ├── held_out
│   │   └── 0.parquet
│   ├── train
│   │   └── 0.parquet
│   └── tuning
│       └── 0.parquet
└── metadata
    ├── codes.parquet
    ├── dataset.json
    └── subject_splits.parquet

160 subjects, 2258 measurements, a ten-code vocabulary — half of it mapped into LOINC, ICD10CM, and SNOMED via parent_codes, half unmapped. examples/build_datasets.py is the readable source of truth for what is in them.

1. Extract a baseline

extract_MEDS_summary_stats examples/datasets/v1_baseline -o examples/stats/v1_baseline.json

That fingerprint is committed at examples/stats/v1_baseline.json — 8 KB of JSON, safe to keep in version control next to the ETL that produced it. The headline block:

>>> from meds_summary_stats.extract import read_summary_stats
>>> baseline = read_summary_stats(EXAMPLES / "stats" / "v1_baseline.json")
>>> print(json.dumps(baseline["counts"], indent=2, sort_keys=True))
{
  "n_events": 2258,
  "n_measurements": 2258,
  "n_static_measurements": 160,
  "n_subjects": 160,
  "n_subjects_with_birth": 160,
  "n_subjects_with_death": 0,
  "n_subjects_with_static": 160,
  "n_unique_codes": 10,
  "numeric": {
    "mean": 103.9,
    "n_finite": 1037,
    "n_infinite": 0,
    "n_nan": 0,
    "n_negative": 0,
    "n_present": 1037,
    "n_zero": 0,
    "std": 26.3615
  },
  "quality": {
    "n_empty_code": 0,
    "n_null_code": 0,
    "n_null_subject_id": 0
  }
}

Splits survive the default suppression floor of 20 subjects, so they are reported by name:

>>> print(json.dumps(baseline["splits"]["by_split"], indent=2, sort_keys=True))
{
  "held_out": {
    "n_measurements": 315,
    "n_subjects": 24
  },
  "train": {
    "n_measurements": 1593,
    "n_subjects": 112
  },
  "tuning": {
    "n_measurements": 350,
    "n_subjects": 24
  }
}

Alongside the statistics, the fingerprint carries a checksum of the dataset itself — a sha256 over the bytes of every file the MEDS schema defines. That answers the prior question, "is this even the same dataset?", without reference to any statistic:

>>> baseline["dataset_digest"]["combined"]
'sha256:...'
>>> baseline["dataset_digest"]["n_files"], baseline["dataset_digest"]["data"]["n_files"]
(6, 3)

Files outside the schema are ignored, so task labels and caches in the dataset root cannot affect it, and metadata/dataset.json is digested separately because its created_at changes on every run. It is a digest of bytes, so it is also sensitive to re-encoding — see docs/stats-surface.md for what it is and is not invariant to, and for how to promote it from a warning to a hard error once you know your ETL is byte-reproducible.

Re-extracting the same dataset reproduces it exactly, which is the property the whole tool rests on:

>>> from meds_summary_stats.compare import compare_summary_stats
>>> from meds_summary_stats.report import render_text
>>> print(render_text(compare_summary_stats(baseline, baseline)))
PASS - 199 paths checked, 4 ignored, no findings.

2. A broken ETL, same source data

examples/datasets/v1_etl_regression is the same source data — dataset_version still 1.0 — put through an ETL carrying three independent bugs: TEMP is emitted in Celsius instead of Fahrenheit, LAB//SODIUM is silently renamed, and every twentieth subject is dropped.

extract_MEDS_summary_stats examples/datasets/v1_etl_regression -o current.json
compare_MEDS_summary_stats examples/stats/v1_baseline.json current.json # exit 1

>>> regressed = read_summary_stats(EXAMPLES / "stats" / "v1_etl_regression.json")
>>> report = compare_summary_stats(baseline, regressed)
>>> report.status, report.failed, report.counts
('fail', True, {'error': 44, 'info': 0, 'warning': 22})

All three bugs are caught, each on a path that names the problem:

>>> findings = {f.path: f for f in report.findings}

>>> print(findings["counts.n_subjects"].message)          # dropped subjects
160 -> 152 (-5.00%), over relative threshold 0.005

>>> findings["codes.digests.alphabetical"].kind           # renamed code
'exact'

>>> loinc = findings["code_metadata.ontology.by_vocabulary.LOINC.numeric.mean"]
>>> loinc.baseline, loinc.current                         # Fahrenheit -> Celsius
(103.9, 88.1595)

Note what the second one does not say. The vocabulary change is detected by a digest, so neither the old code nor the new one appears anywhere in the report:

>>> blob = json.dumps(report.to_dict())
>>> "SODIUM" in blob or "LAB//NA" in blob
False

And note what the third one does. The per-code moments digest says only "some code's distribution moved"; the per-vocabulary moments narrow it to the LOINC-mapped measurements. That is as specific as a report can get without naming a code.

3. A genuine data refresh

examples/datasets/v2_data_refresh bumps dataset_version to 2.0 and has 220 subjects instead of 160 — same schema, same vocabulary, more data. Strictly, that fails:

>>> refreshed = read_summary_stats(EXAMPLES / "stats" / "v2_data_refresh.json")
>>> compare_summary_stats(baseline, refreshed).failed
True

Which is the wrong answer, and why relax exists. It tells the comparison that a version bump excuses volume drift but nothing else:

compare_MEDS_summary_stats examples/stats/v1_baseline.json current.json \
	--on-dataset-version-change relax # exit 0

>>> from meds_summary_stats.policy import ComparePolicy
>>> policy = ComparePolicy.from_dict({"on_dataset_version_change": "relax"})
>>> relaxed = compare_summary_stats(baseline, refreshed, policy)
>>> relaxed.status, relaxed.failed, relaxed.counts
('warn', False, {'error': 0, 'info': 60, 'warning': 11})

The 60 volume findings drop to informational. What survives as a warning is the ordering of the vocabulary by frequency, which genuinely did shift as the cohort grew, plus the file checksums, which are of course different files:

>>> sorted({f.path.split(".")[0] for f in relaxed.findings if str(f.severity) == "warning"})
['codes', 'dataset_digest']
>>> sorted(f.path for f in relaxed.findings
...        if str(f.severity) == "warning" and f.path.startswith("codes."))
['codes.digests.by_measurement_frequency', 'codes.digests.by_subject_frequency',
 'codes.digests.numeric_moments.digest']

Crucially, relax does not blind the check. The set of codes is unchanged here, and if the refresh had also broken the vocabulary or the schema, that would still be a hard error:

>>> baseline["codes"]["digests"]["alphabetical"] == refreshed["codes"]["digests"]["alphabetical"]
True
>>> from copy import deepcopy
>>> also_broken = deepcopy(refreshed)
>>> also_broken["layout"]["data_columns"]["subject_id"] = "Int32"
>>> broken_report = compare_summary_stats(baseline, also_broken, policy)
>>> broken_report.failed
True
>>> [f.path for f in broken_report.findings if str(f.severity) == "error"]
['layout.data_columns.subject_id']

The two failure modes are distinguishable, which is the whole point: a data refresh moves volume, a broken ETL moves volume and the vocabulary and a distribution.

What it emits

A single canonical JSON document. Roughly:

docs/stats-surface.md documents every field, with live output.

The privacy contract

The output is designed to be committable to a public repository as a baseline. Three rules, enforced structurally rather than by convention:

1. No code string is ever emitted — at any threshold, under any flag. The vocabulary appears only as digests (sha256 over the sorted code list, over the frequency-ordered list, and over frequency-restricted subsets) and as aggregate distributions. That is enough to detect that the vocabulary changed, and never enough to reconstruct it. Nor is the vocabulary ever partitioned by the structure of the code strings: splitting on the common // convention is not done, both because MEDS does not require it and because for a flat vocabulary the prefix is the whole code, so the list of prefixes would be the list of codes.
2. No small cells. Any split, calendar-year, or ontology-vocabulary entry backed by fewer than min_cell_size distinct subjects (default 20) is dropped entirely — the key is never written, only a <suppressed> rollup of how many were dropped and their pooled volume.
3. No extremes. No min, no max, anywhere. Not on values, not on timestamps. The subject with the longest record is exactly where re-identification starts; quantiles replace extremes, and are themselves suppressed below min_numeric_cell_size observations.

The one exception is deliberate, and it is where the deeper vocabulary inspection lives: parent_codes entries name concepts in published external terminologies (LOINC/8867-4), which are properties of the terminology rather than of any patient — and, decisively, the grouping is declared by the dataset rather than guessed from string shape. So that block reports vocabularies by name, with per-vocabulary volume, subject counts, numeric moments, and key-set digests (still subject to suppression), plus dataset-level mapping coverage. It is the only axis that says where the vocabulary moved, and it catches a regression nothing else can: when a concept-mapping join silently stops resolving, every count over the data is unchanged and only ontology.coverage.frac_measurements_mapped falls.

How comparison works

Both documents flatten to path -> value. Each path resolves to exactly one rule — the most specific glob that matches it — and the rule says how to compare:

Each rule has a warn and an error tier. The built-in policy is already a reasonable ETL regression policy; a config file only states your overrides:

on_dataset_version_change: relax
fail_on: error
rules:
  - {path: counts.n_subjects, kind: relative, warn: 0.0, error: 0.002}
  - {path: codes.digests.alphabetical, kind: exact, on_change: error}
  - {path: time.by_year.*.n_measurements, kind: relative, error: 0.05}

See examples/policy.yaml for a fully commented one.

on_dataset_version_change: relax

The setting that makes the motivating workflow practical. The premise of the check is "same data, new ETL". When the data itself moved — dataset.dataset_version differs — volume drift is expected and flagging it is noise, but a schema or vocabulary regression is still a bug. relax draws that line: everything numeric drops to informational, exact rules stay hard.

>>> from copy import deepcopy
>>> from meds_summary_stats.policy import ComparePolicy
>>> baseline = {"stats_schema_version": 1, "config": {},
...             "dataset": {"dataset_version": "1.0"},
...             "layout": {"data_columns": {"code": "String"}},
...             "counts": {"n_subjects": 100_000}}
>>> refreshed = deepcopy(baseline)
>>> refreshed["dataset"]["dataset_version"] = "2.0"
>>> refreshed["counts"]["n_subjects"] = 140_000

Strictly, a 40% jump in subjects is a failure:

>>> compare_summary_stats(baseline, refreshed).status
'fail'

Under `relax`, it is expected -- the data changed:

>>> policy = ComparePolicy.from_dict({"on_dataset_version_change": "relax"})
>>> compare_summary_stats(baseline, refreshed, policy).status
'pass'

But a schema change is never excused:

>>> broken = deepcopy(refreshed)
>>> broken["layout"]["data_columns"]["code"] = "Int64"
>>> report = compare_summary_stats(baseline, broken, policy)
>>> report.status, [f.path for f in report.findings if str(f.severity) == "error"]
('fail', ['layout.data_columns.code'])

Reading a report

>>> from meds_summary_stats.report import render_text
>>> regressed = deepcopy(baseline)
>>> regressed["counts"]["n_subjects"] = 98_000
>>> print(render_text(compare_summary_stats(baseline, regressed)))
FAIL - 3 paths checked, 0 ignored; 1 error.
<BLANKLINE>
ERROR (1)
  counts.n_subjects
      100000 -> 98000 (-2.00%), over relative threshold 0.005

compare_MEDS_summary_stats also writes a machine-readable report (--output) and a GitHub-flavored markdown one (--markdown), and appends the markdown to $GITHUB_STEP_SUMMARY when it is running inside Actions.

In CI

A composite GitHub Action, usable from any repository:

  - name: Check for summary-stat regressions
    uses: Medical-Event-Data-Standard/meds_summary_stats@v0
    with:
      meds-dir: build/meds
      baseline: baselines/mimic-iv-3.1.json
      policy: .github/meds-stats-policy.yaml
      package-spec: meds_summary_stats==0.1.0

docs/github-action.md covers the reusable-workflow form, advisory rollout with fail-on: never, reacting to the outcome, label-gated baseline refresh, and the full input/output/exit-code tables.

Performance

One scan_parquet per shard, normalized and concatenated, then a handful of aggregates collected together on the polars streaming engine so the common scan is shared. The dataset is read a small fixed number of times regardless of its size, and the only frames materialized are per-subject and per-code — both orders of magnitude smaller than the data.

The default numeric-detail: moments needs no sorting. quantiles adds full-column quantiles over numeric_value and time, which does require materializing those columns; it is opt-in for that reason. Per-subject and per-code quantiles are always computed, since those frames are small.

Installation

uv add meds_summary_stats # or: uvx --from meds_summary_stats extract_MEDS_summary_stats --help

Requires Python 3.12+.

Contributing

See CONTRIBUTORS.md. Briefly: uv sync --group dev, then uv run pytest -v. Tests are predominantly doctests — the API-level behavior of essentially every function is pinned by examples in its docstring, with tests/ reserved for end-to-end runs and hypothesis properties (notably, that no code string ever escapes into the output for any dataset).

License

MIT