homelab/stack
1
0
Fork 0
Code Issues 57 Pull Requests 1 Actions Packages Projects 1 Releases Wiki Activity

feat(pfs): event study around code-family exposures — pfs.event_study + notebooks/event_study.py (refs #696)

Browse Source 
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.
This commit is contained in:
kert
2026-09-22 16:40:06 -04:00
parent 2d5d12b544
commit 92ae58a2ec
4 changed files with 1155 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

555
notebooks/event_study.py Normal file
View File

@@ -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

265
src/pfs/event_study.py Normal file
View File

@@ -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

72
tests/notebooks/test_event_study_nb.py Normal file
View File

@@ -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

263
tests/pfs/test_event_study.py Normal file
View File

@@ -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 == {}
)