92ae58a2ec
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.
556 lines
21 KiB
Python
556 lines
21 KiB
Python
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
|