From 92ae58a2ec94674e2694d51ff4b1dcb4e55c7f9e Mon Sep 17 00:00:00 2001 From: kert Date: Tue, 22 Sep 2026 16:40:06 -0400 Subject: [PATCH] =?UTF-8?q?feat(pfs):=20event=20study=20around=20code-fami?= =?UTF-8?q?ly=20exposures=20=E2=80=94=20pfs.event=5Fstudy=20+=20notebooks/?= =?UTF-8?q?event=5Fstudy.py=20(refs=20#696)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit pfs.event_study: pre/post windows around an exposure year (the effective date is 1 January, so the exposure year is the first treated year), log ratio of window means, difference in log changes against the never-treated control set from pfs.exposure.control_codes (same RVU status and decile band, no exposure within ±2 years), 95% percentile-bootstrap interval over the controls with a fixed seed. A code with no pre-period of its own (a new code or a successor) is measured through its family's summed series (scope="family"), so G2058 → 99439 is one continuous line. pre_post() gives single-series window means for outcomes without a code-level control set (occupation wages, employment). utilization_panel() ties an exposure row, its pfs.utilization series, its controls and the estimate. notebooks/event_study.py: pick one of nine fixture exposures, see the treated services series with the exposure rule, the indexed control band (10th–90th pct), the contrast table with its interval and FR anchor, clinical-staff wages and health-care employment around the same year, and a table of every fixture contrast with anchors. A 2025 exposure has no post-period until the 2025 PUF ships — stated, never filled. Descriptive contrasts only; the notebook says so and the chat only links here. Market/industry series are deferred: no keyless quote source is reachable from the rack. --- notebooks/event_study.py | 555 +++++++++++++++++++++++++ src/pfs/event_study.py | 265 ++++++++++++ tests/notebooks/test_event_study_nb.py | 72 ++++ tests/pfs/test_event_study.py | 263 ++++++++++++ 4 files changed, 1155 insertions(+) create mode 100644 notebooks/event_study.py create mode 100644 src/pfs/event_study.py create mode 100644 tests/notebooks/test_event_study_nb.py create mode 100644 tests/pfs/test_event_study.py diff --git a/notebooks/event_study.py b/notebooks/event_study.py new file mode 100644 index 0000000..7fba744 --- /dev/null +++ b/notebooks/event_study.py @@ -0,0 +1,555 @@ +import marimo + +__generated_with = "0.23.13" +app = marimo.App(width="medium") + + +@app.cell(hide_code=True) +def _(): + import marimo as mo + + return (mo,) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(""" + # Event study: what happens around a code-family exposure + + The physician fee schedule is law, and `pfs.code_exposure` is its exposure calendar: + the dates on which a code became payable, was revalued, changed status, ended, or was + listed for telehealth, each anchored to the Federal Register paragraph that says so. + This notebook puts pre/post windows around one exposure at a time on three outcome + series — services billed (`pfs.utilization`), clinical-staff wages (`bls.oews`) and + health-care employment (`bls.ces`) — and, for utilization, compares the treated code + with a never-treated control set drawn from codes of the same RVU status and decile + band that had no exposure of their own in the window. A code with no pre-period of + its own — a new code, or a successor such as 99439 — is measured through its family's + summed series, so G2058 → 99439 reads as one continuous line rather than a birth. + + **What the numbers are.** The utilization estimate is a difference in log changes: + the treated code's log ratio of post-window to pre-window mean services, minus the mean + of the same quantity over its controls; the interval is a 95% percentile bootstrap over + the control set. It is a descriptive contrast around a dated legal event, reported with + its uncertainty. It is not a causal estimate, and the chat never states one from it. + """) + return + + +@app.cell(hide_code=True) +def _(): + # ── 0. Setup ── + import altair as alt + import polars as pl + + from conf import connect, path + + connect.theme() + NOTES = {} + _primary = path("db.aco") + REPLICA_PATH = _primary.parent / "replica" / f"{_primary.stem}.ro.duckdb" + + def _open_replica(): + try: + return connect.duckdb("aco", read_only=True) + except Exception as e: # noqa: BLE001 — degrade, never crash the page + NOTES["replica"] = f"replica unavailable: {e}" + return None + + def _open_bib(): + try: + return connect.bib() + except Exception as e: # noqa: BLE001 + NOTES["bib"] = f"bibliography unavailable: {e}" + return None + + con = _open_replica() + store = _open_bib() + + def fr_md(item_key, p_id): + if not item_key: + return "_unanchored_" + if store is None: + return f"{item_key} ¶{p_id}" + try: + from bib.frlink import md_link + + return md_link( + f"p-{p_id}", store=store, item_key=item_key, text=f"{item_key} ¶{p_id}" + ) + except Exception: # noqa: BLE001 + return f"{item_key} ¶{p_id}" + + def not_built(cmd): + return f"_Not built yet — run `{cmd}` and republish the replica._" + + return NOTES, REPLICA_PATH, alt, con, fr_md, not_built, pl, store + + +@app.cell(hide_code=True) +def _(con, mo, not_built): + # ── 1. Pick an exposure ── + from pfs.codetables import read_exposures as _read_exposures + + #: The fixture exposures from #693/#696 — the CCM 2015 and APCM 2025 panels the + #: nightly run must render, plus the other dated events the calendar was built on. + FIXTURES = [ + ("99490", 2015, "becomes-payable"), + ("99487", 2017, "becomes-payable"), + ("99490", 2022, "revalued"), + ("99439", 2021, "becomes-payable"), + ("G2058", 2021, "ends"), + ("99424", 2022, "becomes-payable"), + ("G2211", 2024, "becomes-payable"), + ("G0556", 2025, "becomes-payable"), + ("99441", 2025, "ends"), + ] + _rows = {} + if con is not None: + try: + for _c, _y, _k in FIXTURES: + for _r in _read_exposures(con, code=_c): + if _r.year == _y and _r.kind == _k: + _rows[(_c, _y, _k)] = _r + except Exception: # noqa: BLE001 — not built yet + _rows = {} + labels = { + f"{c} {k} {y}" + + (f" [{_rows[(c, y, k)].family}]" if (c, y, k) in _rows else ""): (c, y, k) + for c, y, k in FIXTURES + } + exposure_picker = mo.ui.dropdown( + options=list(labels), value=list(labels)[0], label="Exposure" + ) + _intro = mo.md( + "## 1. Pick an exposure\n\n" + "Eight dated events from the exposure calendar. Each is a row of " + "`pfs.code_exposure` with an effective date, the code's RVU status and decile " + "band that year, and the FR paragraph it is anchored to." + + ("" if _rows else "\n\n" + not_built("stack pfs exposure --write")) + ) + mo.vstack([_intro, exposure_picker]) + return FIXTURES, exposure_picker, labels + + +@app.cell(hide_code=True) +def _(exposure_picker, labels): + code, year, kind = labels[exposure_picker.value] + return code, kind, year + + +@app.cell(hide_code=True) +def _(alt, code, con, fr_md, kind, math_exp, mo, not_built, pl, year): + # ── 2. Utilization ── + from pfs.event_study import utilization_panel as _panel + + _p = None + if con is not None: + try: + _p = _panel(con, code, year, kind) + except Exception as _e: # noqa: BLE001 — a missing table is "not built yet" + _p = None + _head = ( + f"## 2. Utilization around {code} {kind} ({year})\n\n" + "Services billed per year for the treated code (solid) against its never-treated " + "controls, each control indexed to its own pre-window mean so the shapes are " + "comparable; the shaded band is the 10th–90th percentile of the indexed controls. " + "The table gives the log-change contrast and its bootstrap interval." + ) + if _p is None: + _view = mo.md(_head + "\n\n" + not_built("stack pfs utilization --write")) + else: + _e = _p.effect + _anchor = ( + fr_md(_p.exposure.item_key, _p.exposure.p_id) + if _p.exposure + else "_no exposure row_" + ) + _eff = _p.exposure.effective_date.isoformat() if _p.exposure else "—" + + def _fmt(v, pct=False): + if v is None: + return "—" + return f"{v:+.1%}" if pct else f"{v:,.0f}" + + _tbl = pl.DataFrame( + { + "quantity": [ + "effective date", + "anchor", + "treated series", + "pre window", + "post window", + "treated pre mean (services)", + "treated post mean (services)", + "treated change", + "control change (mean)", + "contrast (treated − control)", + "95% interval", + "controls", + ], + "value": [ + _eff, + _anchor, + f"{code} alone" + if _e.scope == "code" + else f"{_p.exposure.family} family total (the code had no pre-period of its own)", + ", ".join(map(str, _e.pre_years)) or "— (none on file)", + ", ".join(map(str, _e.post_years)) or "— (none on file)", + _fmt(_e.treated_pre), + _fmt(_e.treated_post), + _fmt(math_exp(_e.treated_change), True) + if _e.treated_change is not None + else "— (no pre-period: the code did not exist or was not billed)", + _fmt(math_exp(_e.control_change), True) + if _e.control_change is not None + else "—", + _fmt(_e.effect_pct, True), + f"{math_exp(_e.ci_low):+.1%} to {math_exp(_e.ci_high):+.1%}" + if _e.ci_low is not None + else "— (needs ≥2 controls and a treated change)", + str(_e.n_controls), + ], + } + ) + _panels = [ + mo.md(_head), + mo.ui.table(_tbl, label=f"pfs.event_study — {code} {year}"), + ] + if _p.treated or _p.controls: + _recs = [] + + def _indexed(series, pre_years): + base = [series[y] for y in pre_years if y in series] + b = sum(base) / len(base) if base else None + return {y: (v / b if b else None) for y, v in series.items()} + + for y, v in sorted(_p.treated.items()): + _recs.append( + { + "year": y, + "series": f"{code} (treated)", + "services": v, + "kind": "treated", + } + ) + _ctrl_idx = { + _c: _indexed(_s, _e.pre_years) for _c, _s in _p.controls.items() + } + _by_year = {} + for _c, _s in _ctrl_idx.items(): + for _y, _v in _s.items(): + if _v is not None: + _by_year.setdefault(_y, []).append(_v) + _band = [ + { + "year": y, + "p10": sorted(v)[int(0.1 * (len(v) - 1))], + "p90": sorted(v)[int(0.9 * (len(v) - 1))], + "median": sorted(v)[len(v) // 2], + } + for y, v in sorted(_by_year.items()) + ] + _t = ( + pl.DataFrame(_recs) + if _recs + else pl.DataFrame( + {"year": [], "series": [], "services": [], "kind": []} + ) + ) + _treated_chart = ( + alt.Chart(_t.to_pandas()) + .mark_line(point=True) + .encode( + x=alt.X("year:O", title="Calendar year"), + y=alt.Y("services:Q", title="Services (treated)"), + tooltip=["year", "services"], + ) + .properties(width=640, height=220, title=f"{code}: services billed") + ) + _rule = ( + alt.Chart(pl.DataFrame({"year": [year]}).to_pandas()) + .mark_rule(strokeDash=[4, 4], color="gray") + .encode(x="year:O") + ) + _panels.append(mo.ui.altair_chart(_treated_chart + _rule)) + if _band: + _b = pl.DataFrame(_band) + _band_chart = ( + alt.Chart(_b.to_pandas()) + .mark_area(opacity=0.25) + .encode( + x=alt.X("year:O", title="Calendar year"), + y=alt.Y("p10:Q", title="Indexed to pre-window mean"), + y2="p90:Q", + ) + .properties( + width=640, + height=220, + title=f"Controls ({_e.n_controls}): indexed services, 10th–90th pct and median", + ) + ) + _med = ( + alt.Chart(_b.to_pandas()) + .mark_line() + .encode(x="year:O", y="median:Q") + ) + _panels.append(mo.ui.altair_chart(_band_chart + _med + _rule)) + else: + _panels.append( + mo.md( + "_No utilization rows on file for this code or its controls yet — for a 2025 exposure that means the 2025 public use file has not been published._" + ) + ) + _view = mo.vstack(_panels) + _view + return + + +@app.cell(hide_code=True) +def _(): + import math as _math + + def math_exp(v): + return _math.exp(v) - 1 if v is not None else None + + return (math_exp,) + + +@app.cell(hide_code=True) +def _(alt, code, con, kind, mo, not_built, pl, year): + # ── 3. Wages and employment ── + from bls.ces import read_series as _read_ces + from bls.oews import read_occupations as _read_oews + from pfs.event_study import pre_post as _pre_post + + _occs = { + "31-9092": "Medical assistants", + "29-2061": "LPN/LVN", + "29-1141": "Registered nurses", + "29-1171": "Nurse practitioners", + } + _series_ids = { + "CES6562110001": "Offices of physicians, employment (k)", + "CES6562100001": "Ambulatory care, employment (k)", + "CES6562110003": "Offices of physicians, avg hourly earnings ($)", + } + _w, _c = [], [] + if con is not None: + try: + _w = _read_oews(con, list(_occs)) + except Exception: # noqa: BLE001 + _w = [] + try: + _c = _read_ces(con, list(_series_ids), annual=True) + except Exception: # noqa: BLE001 + _c = [] + _head = ( + f"## 3. Clinical-staff wages and employment around {year}\n\n" + "These series are national, not code-specific, so there is no control set: the " + "table reports the pre- and post-window means and the log change for each, " + "which is context for the utilization contrast above, not an estimate of the " + "exposure's effect." + ) + if not _w and not _c: + _view = mo.md( + _head + + "\n\n" + + not_built("stack bls oews --write` and `stack bls ces --write") + ) + else: + _rows = [] + for occ, name in _occs.items(): + s = { + r.year: r.a_mean + for r in _w + if r.occ_code == occ and r.a_mean is not None + } + if s: + p = _pre_post(s, year) + _rows.append( + { + "series": f"{occ} {name}: mean annual wage ($)", + "pre": ", ".join(map(str, p.pre_years)), + "pre mean": p.pre_mean, + "post": ", ".join(map(str, p.post_years)), + "post mean": p.post_mean, + "log change": p.change, + } + ) + for sid, name in _series_ids.items(): + s = { + r.year: r.value + for r in _c + if r.series_id == sid and r.value is not None + } + if s: + p = _pre_post(s, year) + _rows.append( + { + "series": f"{sid} {name}", + "pre": ", ".join(map(str, p.pre_years)), + "pre mean": p.pre_mean, + "post": ", ".join(map(str, p.post_years)), + "post mean": p.post_mean, + "log change": p.change, + } + ) + _tbl = pl.DataFrame(_rows) if _rows else pl.DataFrame({"series": []}) + _wdf = pl.DataFrame( + { + "occupation": [f"{r.occ_code} {_occs[r.occ_code]}" for r in _w], + "year": [r.year for r in _w], + "mean annual wage": [r.a_mean for r in _w], + } + ) + _rule = ( + alt.Chart(pl.DataFrame({"year": [year]}).to_pandas()) + .mark_rule(strokeDash=[4, 4], color="gray") + .encode(x="year:O") + ) + _panels = [mo.md(_head), mo.ui.table(_tbl, label="pre/post window means")] + if len(_wdf): + _chart = ( + alt.Chart(_wdf.to_pandas()) + .mark_line(point=True) + .encode( + x=alt.X("year:O", title="OEWS release year"), + y=alt.Y("mean annual wage:Q", title="Mean annual wage ($)"), + color=alt.Color("occupation:N", title="Occupation"), + tooltip=["occupation", "year", "mean annual wage"], + ) + .properties(width=640, height=260, title="Clinical-staff wages (OEWS)") + ) + _panels.append(mo.ui.altair_chart(_chart + _rule)) + _view = mo.vstack(_panels) + _view + return + + +@app.cell(hide_code=True) +def _(FIXTURES, con, fr_md, mo, not_built, pl): + # ── 4. All fixture exposures ── + from pfs.event_study import utilization_panel as _panel_all + + _rows = [] + if con is not None: + for _c, _y, _k in FIXTURES: + try: + _pn = _panel_all(con, _c, _y, _k) + except Exception: # noqa: BLE001 + continue + _pe = _pn.effect + _rows.append( + { + "exposure": f"{_c} {_k} {_y}", + "family": _pn.exposure.family if _pn.exposure else "", + "treated": _pe.scope, + "anchor": fr_md(_pn.exposure.item_key, _pn.exposure.p_id) + if _pn.exposure + else "—", + "pre": ", ".join(map(str, _pe.pre_years)), + "post": ", ".join(map(str, _pe.post_years)), + "contrast": f"{_pe.effect_pct:+.1%}" + if _pe.effect_pct is not None + else "—", + "95% interval": f"{__import__('math').exp(_pe.ci_low) - 1:+.1%} to {__import__('math').exp(_pe.ci_high) - 1:+.1%}" + if _pe.ci_low is not None + else "—", + "controls": _pe.n_controls, + "note": "" + if _pe.effect is not None + else ( + "no pre-period (new code)" + if _pe.treated_change is None + and _pe.pre_years + and not _pe.treated_pre + else ( + "no post-period on file yet" + if not _pe.post_years or _pe.treated_post is None + else "insufficient controls" + ) + ), + } + ) + _head = ( + "## 4. All fixture exposures\n\n" + "The same contrast for every fixture event, with the FR anchor that dates it. A new " + "code has no pre-period by construction — its first years are the outcome, not a " + "change — and a 2025 exposure has no post-period until CMS publishes the 2025 file; " + "both are stated rather than filled in." + ) + _view = ( + mo.vstack( + [ + mo.md(_head), + mo.ui.table(pl.DataFrame(_rows), label="event-study contrasts"), + ] + ) + if _rows + else mo.md( + _head + + "\n\n" + + not_built( + "stack pfs exposure --write` and `stack pfs utilization --write" + ) + ) + ) + _view + return + + +@app.cell(hide_code=True) +def _(NOTES, REPLICA_PATH, con, mo, pl): + # ── 5. Provenance ── + import datetime as _dt + + if REPLICA_PATH.exists(): + _mtime = _dt.datetime.fromtimestamp(REPLICA_PATH.stat().st_mtime).isoformat( + timespec="seconds" + ) + _line = f"`{REPLICA_PATH}` — last published {_mtime}" + else: + _line = f"`{REPLICA_PATH}` — not found" + + def _q(sql): + if con is None: + return pl.DataFrame() + try: + return con.execute(sql).pl() + except Exception as e: # noqa: BLE001 + NOTES[sql[:40]] = str(e) + return pl.DataFrame() + + _counts = _q( + "SELECT 'pfs.code_exposure' t, count(*) n FROM pfs.code_exposure " + "UNION ALL SELECT 'pfs.utilization', count(*) FROM pfs.utilization " + "UNION ALL SELECT 'bls.oews', count(*) FROM bls.oews " + "UNION ALL SELECT 'bls.ces', count(*) FROM bls.ces" + ) + _log = _q( + "SELECT ingested_at, module, table_name, rows, pincite_key FROM cms.ingest_log " + "WHERE module IN ('pfs.utilization', 'bls.oews', 'bls.ces') ORDER BY ingested_at DESC LIMIT 30" + ) + _notes = "\n".join(f"- {k}: {v}" for k, v in NOTES.items()) + mo.vstack( + [ + mo.md( + "## 5. Provenance\n\n" + f"Replica: {_line}. Every outcome row carries the bib key of the dataset " + "release it came from (`item_key`), and each release is a Source item " + "tagged with its table; `cms.ingest_log` records the ingest with the source " + "file's sha256. Market and industry series (sector ETFs, telehealth and " + "care-management tickers) are not on file: no keyless quote source is " + "reachable from the rack, so that panel waits on a data source (#696 follow-up)." + + (f"\n\n{_notes}" if _notes else "") + ), + mo.ui.table(_counts, label="row counts") + if len(_counts) + else mo.md("_row counts unavailable_"), + mo.ui.table(_log, label="cms.ingest_log — outcome series") + if len(_log) + else mo.md("_no ingest log rows_"), + ] + ) + return diff --git a/src/pfs/event_study.py b/src/pfs/event_study.py new file mode 100644 index 0000000..44d1fa9 --- /dev/null +++ b/src/pfs/event_study.py @@ -0,0 +1,265 @@ +"""Event study around code-family exposures (#696): pre/post windows on an +outcome series with never-treated control sets. + +The fee schedule is law; ``pfs.code_exposure`` is the exposure calendar +(P50). For an exposure in *year*, the window is ``pre`` file years before +it and ``post`` years from it on (the effective date is 1 January, so +the exposure year is the first treated year). The estimate is a +difference in log changes:: + + effect = log(post_mean/pre_mean)[treated] − mean_c log(post_mean/pre_mean)[control c] + +with a 95% percentile-bootstrap interval from resampling the control set +(``_BOOT`` draws, fixed seed — the same inputs always give the same +interval). It is a descriptive contrast, not a causal estimate: the +notebook says so and the chat only links here. + + windows(year, available, pre, post) → (pre_years, post_years) + log_change(series, pre, post) → log ratio of window means + did_effect(treated, controls, …) → Effect + pre_post(series, year, …) → PrePost (single-series, no controls) + utilization_panel(con, code, year, kind) → Panel over pfs.utilization + family_series(con, family) → member codes summed per year +""" + +from __future__ import annotations + +import math +import random +from dataclasses import dataclass +from typing import Any, Iterable, Mapping, Sequence + +_BOOT = 500 +_SEED = 20260922 + + +@dataclass(frozen=True) +class Effect: + code: str + year: int + kind: str + outcome: str + pre_years: tuple[int, ...] + post_years: tuple[int, ...] + treated_pre: float | None + treated_post: float | None + treated_change: float | None + control_change: float | None + effect: float | None + ci_low: float | None + ci_high: float | None + n_controls: int + #: ``"code"`` when the treated series is the code's own utilization; + #: ``"family"`` when the code had no pre-period of its own (a new + #: code, or a successor) and the family's summed series stood in — + #: G2058 → 99439 is then one continuous series, not a birth. + scope: str = "code" + + @property + def effect_pct(self) -> float | None: + return math.exp(self.effect) - 1 if self.effect is not None else None + + +@dataclass(frozen=True) +class PrePost: + year: int + pre_years: tuple[int, ...] + post_years: tuple[int, ...] + pre_mean: float | None + post_mean: float | None + change: float | None + + +@dataclass(frozen=True) +class Panel: + exposure: Any # ExposureRow | None + effect: Effect + treated: dict[int, float] + controls: dict[str, dict[int, float]] + + +# ── pure pieces ─────────────────────────────────────────────────────── + + +def windows( + year: int, *, available: Iterable[int], pre: int = 2, post: int = 2 +) -> tuple[tuple[int, ...], tuple[int, ...]]: + have = set(available) + return ( + tuple(y for y in range(year - pre, year) if y in have), + tuple(y for y in range(year, year + post) if y in have), + ) + + +def _mean(vals: Sequence[float]) -> float | None: + return sum(vals) / len(vals) if vals else None + + +def _window_means( + series: Mapping[int, float], pre: Sequence[int], post: Sequence[int] +) -> tuple[float | None, float | None]: + pv = [series[y] for y in pre if series.get(y) is not None] + qv = [series[y] for y in post if series.get(y) is not None] + return _mean(pv), _mean(qv) + + +def log_change( + series: Mapping[int, float], pre: Sequence[int], post: Sequence[int] +) -> float | None: + a, b = _window_means(series, pre, post) + if a is None or b is None or a <= 0 or b <= 0: + return None + return math.log(b / a) + + +def did_effect( + treated: Mapping[int, float], + controls: Mapping[str, Mapping[int, float]], + year: int, + *, + code: str, + kind: str, + outcome: str, + pre: int = 2, + post: int = 2, +) -> Effect: + years = set(treated) | {y for s in controls.values() for y in s} + pre_years, post_years = windows(year, available=years, pre=pre, post=post) + t_pre, t_post = _window_means(treated, pre_years, post_years) + t_change = log_change(treated, pre_years, post_years) + changes = [ + c + for c in (log_change(s, pre_years, post_years) for s in controls.values()) + if c is not None + ] + c_change = _mean(changes) + effect = ( + t_change - c_change if t_change is not None and c_change is not None else None + ) + lo = hi = None + if effect is not None and len(changes) >= 2: + rng = random.Random(_SEED) + n = len(changes) + draws = sorted( + t_change - sum(rng.choice(changes) for _ in range(n)) / n + for _ in range(_BOOT) + ) + lo, hi = draws[int(0.025 * _BOOT)], draws[int(0.975 * _BOOT) - 1] + return Effect( + code=code, + year=year, + kind=kind, + outcome=outcome, + pre_years=pre_years, + post_years=post_years, + treated_pre=t_pre, + treated_post=t_post, + treated_change=t_change, + control_change=c_change, + effect=effect, + ci_low=lo, + ci_high=hi, + n_controls=len(changes), + ) + + +def pre_post( + series: Mapping[int, float], year: int, *, pre: int = 2, post: int = 2 +) -> PrePost: + """A single series' window means around *year* — for outcomes that + have no code-level control set (occupation wages, employment).""" + pre_years, post_years = windows(year, available=series, pre=pre, post=post) + a, b = _window_means(series, pre_years, post_years) + return PrePost( + year, pre_years, post_years, a, b, log_change(series, pre_years, post_years) + ) + + +# ── panels over the replica ─────────────────────────────────────────── + + +def utilization_panel( + con: Any, + code: str, + year: int, + kind: str, + *, + outcome: str = "n_services", + pre: int = 2, + post: int = 2, + window: int = 2, + max_controls: int = 60, +) -> Panel: + """The exposure row, its ``pfs.utilization`` series, the never-treated + control series (``pfs.exposure.control_codes``) and the estimate.""" + from pfs.codetables import read_exposures + from pfs.exposure import control_codes + from pfs.utilization import read_series + + code = code.upper() + exposure = next( + ( + r + for r in read_exposures(con, code=code) + if r.year == year and r.kind == kind + ), + None, + ) + try: + ctrl = control_codes(con, code, year, window=window, limit=max_controls) + except Exception: # noqa: BLE001 — no exposure calendar yet + ctrl = [] + ctrl_codes = [c for c, *_ in ctrl] + rows = read_series(con, [code, *ctrl_codes]) + series: dict[str, dict[int, float]] = {} + for r in rows: + v = getattr(r, outcome) + if v is not None: + series.setdefault(r.hcpcs, {})[r.year] = float(v) + treated = series.get(code, {}) + controls = {c: series[c] for c in ctrl_codes if c in series} + scope = "code" + pre_years, _ = windows(year, available=range(1900, 2200), pre=pre, post=post) + if ( + exposure is not None + and exposure.family + and not any(y in treated for y in pre_years) + ): + fam_series = family_series(con, exposure.family, outcome=outcome) + if any(y in fam_series for y in pre_years): + treated, scope = fam_series, "family" + effect = did_effect( + treated, + controls, + year, + code=code, + kind=kind, + outcome=outcome, + pre=pre, + post=post, + ) + effect = Effect(**{**effect.__dict__, "scope": scope}) + return Panel(exposure=exposure, effect=effect, treated=treated, controls=controls) + + +def family_series( + con: Any, family: str, *, outcome: str = "n_services" +) -> dict[int, float]: + """The family's member codes (``pfs.code_family``) summed per year.""" + from pfs.utilization import read_series + + try: + codes = [ + r[0] + for r in con.execute( + "SELECT DISTINCT code FROM pfs.code_family WHERE key = ?", [family] + ).fetchall() + ] + except Exception: # noqa: BLE001 — no family table on this replica + return {} + out: dict[int, float] = {} + for r in read_series(con, codes): + v = getattr(r, outcome) + if v is not None: + out[r.year] = out.get(r.year, 0.0) + float(v) + return out diff --git a/tests/notebooks/test_event_study_nb.py b/tests/notebooks/test_event_study_nb.py new file mode 100644 index 0000000..29426f6 --- /dev/null +++ b/tests/notebooks/test_event_study_nb.py @@ -0,0 +1,72 @@ +"""notebooks/event_study.py — structure, headless degradation, live fixture panels (#696).""" + +from __future__ import annotations + +import ast +import importlib.util +import re +from pathlib import Path + +import pytest + +NB = Path(__file__).resolve().parents[2] / "notebooks" / "event_study.py" +_ROOT = Path(__file__).resolve().parents[2] +_HAS_DATA = (_ROOT / "data" / "replica" / "aco.ro.duckdb").exists() and ( + _ROOT / "data" / "bib.sqlite" +).exists() +_live = pytest.mark.skipif(not _HAS_DATA, reason="needs the live replica and bib") + + +def _load(): + spec = importlib.util.spec_from_file_location("event_study_nb", NB) + mod = importlib.util.module_from_spec(spec) + spec.loader.exec_module(mod) + return mod + + +def test_notebook_is_a_marimo_app(): + assert _load().app.__class__.__name__ == "App" + + +def test_cells_are_anonymous_and_banners_present(): + src = NB.read_text() + names = [n.name for n in ast.parse(src).body if isinstance(n, ast.FunctionDef)] + assert names and set(names) == {"_"} + assert re.findall(r"# ── (\S+)\. ", src) == ["0", "1", "2", "3", "4", "5"] + + +def test_no_causal_language(): + src = NB.read_text().lower() + assert ( + "caused" not in src and "causal estimate" in src + ) # the disclaimer, nothing else + + +def test_headless_run_degrades_without_data(monkeypatch): + monkeypatch.setenv("STACK_DUCKDB_REPLICA", "1") + mod = _load() + import conf.connect as cc + + monkeypatch.setattr( + cc, + "duckdb", + lambda *a, **k: (_ for _ in ()).throw(FileNotFoundError("no replica")), + ) + monkeypatch.setattr( + cc, "bib", lambda *a, **k: (_ for _ in ()).throw(FileNotFoundError("no bib")) + ) + outputs, _ = mod.app.run() + rendered = "\n".join(o._repr_html_() for o in outputs if hasattr(o, "_repr_html_")) + assert "Not built yet" in rendered and "Traceback" not in rendered + + +@_live +def test_ccm_2015_and_apcm_2025_panels_render_with_anchors(): + mod = _load() + outputs, _ = mod.app.run() + rendered = "\n".join(o._repr_html_() for o in outputs if hasattr(o, "_repr_html_")) + assert "Traceback" not in rendered + # section 4 lists every fixture with its anchor: CCM 2015 (DE2VH9PD) and APCM 2025 (JJ6AM5HJ) + assert "99490 becomes-payable 2015" in rendered and "DE2VH9PD" in rendered + assert "G0556 becomes-payable 2025" in rendered and "JJ6AM5HJ" in rendered + assert "no post-period on file yet" in rendered diff --git a/tests/pfs/test_event_study.py b/tests/pfs/test_event_study.py new file mode 100644 index 0000000..c62f4f6 --- /dev/null +++ b/tests/pfs/test_event_study.py @@ -0,0 +1,263 @@ +"""pfs.event_study — pre/post windows around exposures with never-treated controls (#696).""" + +from __future__ import annotations + +import datetime as dt +import math + +import duckdb +import pytest + +from pfs.codetables import ExposureRow, ensure_tables, write_exposures +from pfs.event_study import ( + Effect, + PrePost, + did_effect, + log_change, + pre_post, + utilization_panel, + windows, +) +from pfs.utilization import UtilizationRow +from pfs.utilization import write as write_util + + +class TestWindows: + def test_pre_before_post_from_exposure_year(self): + pre, post = windows(2017, available=range(2013, 2025)) + assert pre == (2015, 2016) and post == (2017, 2018) + + def test_truncated_at_the_data_edge(self): + pre, post = windows(2024, available=range(2013, 2025), pre=3, post=3) + assert pre == (2021, 2022, 2023) and post == (2024,) + assert windows(2013, available=range(2013, 2025)) == ((), (2013, 2014)) + + +class TestLogChange: + def test_ratio_of_means(self): + s = {2015: 100.0, 2016: 100.0, 2017: 200.0, 2018: 200.0} + assert log_change(s, (2015, 2016), (2017, 2018)) == pytest.approx(math.log(2)) + + def test_missing_side_is_none(self): + assert log_change({2017: 1.0}, (2015, 2016), (2017,)) is None + assert log_change({2015: 1.0}, (2015,), (2017,)) is None + assert ( + log_change({2015: 0.0, 2017: 1.0}, (2015,), (2017,)) is None + ) # zero pre → undefined + + +class TestDidEffect: + def _controls(self, n=20, growth=0.1): + return { + f"C{i:04d}": { + 2015: 100.0, + 2016: 100.0, + 2017: 100.0 * (1 + growth + 0.01 * (i % 5)), + 2018: 100.0 * (1 + growth + 0.01 * (i % 5)), + } + for i in range(n) + } + + def test_effect_is_treated_change_minus_control_change(self): + treated = {2015: 100.0, 2016: 100.0, 2017: 300.0, 2018: 300.0} + e = did_effect( + treated, + self._controls(), + 2017, + code="X", + kind="becomes-payable", + outcome="services", + ) + assert isinstance(e, Effect) + assert e.pre_years == (2015, 2016) and e.post_years == (2017, 2018) + assert e.treated_change == pytest.approx(math.log(3)) + assert 0.09 < e.control_change < 0.14 + assert e.effect == pytest.approx(e.treated_change - e.control_change) + assert e.ci_low is not None and e.ci_low < e.effect < e.ci_high + assert e.n_controls == 20 + assert e.effect_pct == pytest.approx(math.exp(e.effect) - 1) + + def test_ci_is_deterministic_and_shrinks_with_more_controls(self): + treated = {2015: 100.0, 2016: 100.0, 2017: 150.0, 2018: 150.0} + a = did_effect( + treated, self._controls(10), 2017, code="X", kind="k", outcome="o" + ) + b = did_effect( + treated, self._controls(10), 2017, code="X", kind="k", outcome="o" + ) + c = did_effect( + treated, self._controls(200), 2017, code="X", kind="k", outcome="o" + ) + assert (a.ci_low, a.ci_high) == (b.ci_low, b.ci_high) + assert (c.ci_high - c.ci_low) < (a.ci_high - a.ci_low) + + def test_too_few_controls_gives_no_ci(self): + treated = {2015: 1.0, 2016: 1.0, 2017: 2.0} + e = did_effect( + treated, + {"C1": {2015: 1.0, 2017: 1.1}}, + 2017, + code="X", + kind="k", + outcome="o", + ) + assert e.n_controls == 1 and e.ci_low is None and e.effect is not None + e0 = did_effect(treated, {}, 2017, code="X", kind="k", outcome="o") + assert e0.n_controls == 0 and e0.control_change is None and e0.effect is None + + def test_treated_without_post_data(self): + e = did_effect( + {2015: 1.0, 2016: 1.0}, + self._controls(), + 2017, + code="X", + kind="k", + outcome="o", + ) + assert e.treated_change is None and e.effect is None and e.n_controls == 20 + + +class TestPrePost: + def test_means_and_change(self): + p = pre_post({2015: 10.0, 2016: 12.0, 2017: 20.0, 2018: 24.0}, 2017) + assert isinstance(p, PrePost) + assert p.pre_mean == 11.0 and p.post_mean == 22.0 + assert p.change == pytest.approx(math.log(2)) + + def test_empty(self): + p = pre_post({}, 2017) + assert p.pre_mean is None and p.change is None + + +RVU_COLS = "hcpcs VARCHAR, mod VARCHAR, description VARCHAR, status_code VARCHAR, non_fac_total DOUBLE, year INTEGER" + + +def _util(code, year, services): + return UtilizationRow( + year, + code, + "d", + "N", + "O", + 10, + 100, + float(services), + float(services), + 1.0, + 2.0, + 3.0, + 4.0, + "ds", + "K", + ) + + +def _exp(code, year, kind, status="A", band=1, item_key="ITEM", p_id=5): + return ExposureRow( + code, + "FAM", + dt.date(year, 1, 1), + year, + kind, + status, + "", + band, + 1.0, + 0.0, + item_key, + p_id, + 9, + "created", + False, + "", + ) + + +@pytest.fixture +def con(): + c = duckdb.connect(":memory:") + ensure_tables(c) + c.execute(f"CREATE TABLE pfs.rvu ({RVU_COLS})") + rows = [] + # treated code T appears payable 2017; controls C1..C6 payable across the span in the same band + for y in range(2015, 2021): + for i in range(1, 7): + rows.append((f"C000{i}", None, "c", "A", 1.0, y)) + if y >= 2017: + rows.append(("T0001", None, "t", "A", 1.0, y)) + c.executemany("INSERT INTO pfs.rvu VALUES (?,?,?,?,?,?)", rows) + write_exposures(c, [_exp("T0001", 2017, "becomes-payable")]) + util = [] + for y in range(2015, 2021): + for i in range(1, 7): + util.append(_util(f"C000{i}", y, 100 * (1 + 0.05 * (y - 2015)))) + if y >= 2017: + util.append(_util("T0001", y, 500 + 100 * (y - 2017))) + for y in range(2015, 2021): + write_util(c, y, [r for r in util if r.year == y]) + yield c + c.close() + + +class TestUtilizationPanel: + def test_panel_uses_control_set_and_carries_anchor(self, con): + panel = utilization_panel(con, "T0001", 2017, "becomes-payable") + assert panel.effect.code == "T0001" and panel.effect.outcome == "n_services" + assert panel.effect.n_controls == 6 + assert ( + panel.exposure is not None + and panel.exposure.item_key == "ITEM" + and panel.exposure.p_id == 5 + ) + # treated has no pre-period utilization (it did not exist) → change undefined, stated not faked + assert panel.effect.treated_change is None and panel.effect.scope == "code" + assert set(panel.controls) == {f"C000{i}" for i in range(1, 7)} + assert panel.treated == {2017: 500.0, 2018: 600.0, 2019: 700.0, 2020: 800.0} + + def test_panel_with_pre_period(self, con): + write_exposures( + con, + [_exp("T0001", 2017, "becomes-payable"), _exp("C0001", 2018, "revalued")], + ) + panel = utilization_panel(con, "C0001", 2018, "revalued") + assert panel.effect.treated_change is not None + assert panel.effect.n_controls == 5 # the other controls, not itself + assert panel.effect.effect == pytest.approx( + 0.0, abs=1e-9 + ) # it grows like its controls + + def test_new_code_falls_back_to_the_family_series(self, con): + # T0001 has no pre-period of its own; its family (FAM) also holds P0001, + # billed before 2017 — the family total stands in as the treated series. + con.executemany( + "INSERT INTO pfs.code_family VALUES (?,?,?,?,?,?,?,?,?)", + [ + ("FAM", "Fam", "T0001", "successor", 2017, None, "", 0, ""), + ("FAM", "Fam", "P0001", "predecessor", 2015, 2016, "", 0, ""), + ], + ) + for y in (2015, 2016): + write_util( + con, + y, + [ + _util(f"C000{i}", y, 100 * (1 + 0.05 * (y - 2015))) + for i in range(1, 7) + ] + + [_util("P0001", y, 400.0)], + ) + panel = utilization_panel(con, "T0001", 2017, "becomes-payable") + assert panel.effect.scope == "family" + assert panel.treated[2015] == 400.0 and panel.treated[2017] == 500.0 + assert ( + panel.effect.treated_change is not None and panel.effect.effect is not None + ) + assert panel.effect.n_controls == 6 + + def test_unknown_exposure_still_returns_a_panel(self, con): + panel = utilization_panel(con, "ZZZZZ", 2018, "ends") + assert ( + panel.exposure is None + and panel.effect.n_controls == 0 + and panel.treated == {} + )