"""Phase 3 - Insight Discovery (extending Agrawal & Choudhary KDD 2013).
Implements the validated 3-tier coefficient-of-variation framework:
- CV_b coefficient of variation of provider billing (avg_submitted_charge)
- CV_p coefficient of variation of provider payment
- CV_nb coefficient of variation of normalized billing
Adds descriptive cross-dataset provider comparisons and urban-rural
stratification. Writes JSONs to site/data/ so the insights.html page can render
choropleths, scatter plots, and bar charts.
"""
import json
from collections import OrderedDict
from pathlib import Path
import numpy as np
import pandas as pd
from .config import (
DATA_PROC, SITE_DATA, MIN_PROVIDERS_PER_SERVICE,
)
# ---------------------------------------------------------------------------
# Analysis metadata and local stability rules
# ---------------------------------------------------------------------------
INPATIENT_PERIOD = "FY2023"
OUTPATIENT_PERIOD = "CY2023"
# The shared config has a service-level rule but no within-state rule. Five
# providers is used here so state means are not driven by one or two hospitals.
MIN_PROVIDERS_PER_STATE_SERVICE = 5
MIN_STATES_FOR_CORRELATION = 4
PRIMARY_GEOGRAPHIC_EXCLUDED_STATES = ("MD",)
_PAYMENT_MEASURE_LABELS = {
"avg_total_payment": "average total payment (beneficiary and third-party amounts included)",
"avg_medicare_payment": "average Medicare payment",
"avg_allowed_amount": "average Medicare allowed amount (including beneficiary share)",
}
def _period_for(dataset_name: str) -> str:
"""Return the CMS period represented by a dataset/analysis track name."""
if "outpatient" in str(dataset_name).lower():
return OUTPATIENT_PERIOD
return INPATIENT_PERIOD
def _payment_measure(payment_col: str) -> str:
return _PAYMENT_MEASURE_LABELS.get(payment_col, payment_col)
def _desc_column(df: pd.DataFrame, service_col: str, desc_col: str = None) -> str:
"""Resolve a service description column without relying on one code name."""
if desc_col and desc_col in df.columns:
return desc_col
candidate = service_col.replace("_code", "_desc")
return candidate if candidate in df.columns else None
def _provider_count(df: pd.DataFrame) -> int:
if "provider_ccn" in df.columns:
return int(df["provider_ccn"].nunique())
return int(len(df))
def _eligible_service_count(df: pd.DataFrame, service_col: str,
billing_col: str, payment_col: str) -> int:
required = [service_col, billing_col, payment_col]
if any(col not in df.columns for col in required):
return 0
valid = df.dropna(subset=required)
if valid.empty:
return 0
if "provider_ccn" in valid.columns:
counts = valid.groupby(service_col)["provider_ccn"].nunique()
else:
counts = valid.groupby(service_col).size()
return int((counts >= MIN_PROVIDERS_PER_SERVICE).sum())
def _analysis_metadata(df: pd.DataFrame, service_col: str, billing_col: str,
payment_col: str, dataset_name: str,
period: str = None, service_universe: str = None,
comparison_eligibility: str = None,
comparison_eligible: bool = False) -> OrderedDict:
"""Build the repeated provenance fields used by every Phase 3 artifact."""
universe_count = int(df[service_col].nunique()) if service_col in df.columns else 0
return OrderedDict(
dataset=dataset_name,
period=period or _period_for(dataset_name),
service_column=service_col,
service_universe=service_universe or (
f"all observed {service_col} codes in the supplied analysis frame"
),
service_universe_count=universe_count,
eligible_service_count=_eligible_service_count(
df, service_col, billing_col, payment_col),
billing_column=billing_col,
payment_column=payment_col,
payment_measure=_payment_measure(payment_col),
minimum_provider_rule=(
f"at least {MIN_PROVIDERS_PER_SERVICE} valid providers per service"
),
minimum_providers_per_service=int(MIN_PROVIDERS_PER_SERVICE),
comparison_eligible=bool(comparison_eligible),
comparison_eligibility=(comparison_eligibility or (
"Exploratory result; not a like-for-like historical comparison"
)),
)
# ---------------------------------------------------------------------------
# Helpers
# ---------------------------------------------------------------------------
def _cv(series: pd.Series) -> float:
"""Coefficient of variation (std / |mean|). NaN-safe. Returns None for degenerate."""
s = series.dropna()
if len(s) < 2:
return None
mu = s.mean()
if mu == 0 or pd.isna(mu):
return None
return float(s.std(ddof=1) / abs(mu))
def _log2_ratio(numerator: float, denominator: float) -> float:
"""log2(x / y) with safety for zero/negatives. NaN inputs → NaN."""
if numerator is None or denominator is None or numerator <= 0 or denominator <= 0:
return float("nan")
return float(np.log2(numerator / denominator))
# ---------------------------------------------------------------------------
# 3a. 3-tier CV per service
# ---------------------------------------------------------------------------
def compute_cv_table(df: pd.DataFrame, service_col: str, desc_col: str,
billing_col: str, payment_col: str, dataset_name: str,
period: str = None, service_universe: str = None,
comparison_eligibility: str = None,
comparison_eligible: bool = False
) -> pd.DataFrame:
"""Per-service CV_b, CV_p, CV_nb with summary stats.
Filters to services with >= MIN_PROVIDERS_PER_SERVICE valid providers.
The broad service universe is intentional for exploration. Metadata on each
row records the eligibility rule and payment concept so rankings are not
mistaken for a controlled historical comparison.
"""
valid = df.dropna(subset=[service_col, billing_col, payment_col]).copy()
metadata = _analysis_metadata(
df, service_col, billing_col, payment_col, dataset_name,
period=period, service_universe=service_universe,
comparison_eligibility=comparison_eligibility,
comparison_eligible=comparison_eligible,
)
rows = []
grouped = valid.groupby(service_col, sort=False)
for code, sub in grouped:
if _provider_count(sub) < MIN_PROVIDERS_PER_SERVICE:
continue
resolved_desc_col = _desc_column(sub, service_col, desc_col)
desc = (sub[resolved_desc_col].dropna().iloc[0]
if resolved_desc_col and sub[resolved_desc_col].notna().any()
else None)
us_b_mu = sub[billing_col].mean()
us_p_mu = sub[payment_col].mean()
cv_b = _cv(sub[billing_col])
cv_p = _cv(sub[payment_col])
# Per-provider normalized billing: billing * (us_p_mu / provider_payment)
# (Agrawal's NF_h = μ_p,h / μ_p,us — we apply the inverse to billing)
sub = sub.copy()
nf = sub[payment_col] / us_p_mu # hospital's payment rel to US
nf = nf.replace([0, np.inf, -np.inf], np.nan)
sub["normalized_billing"] = (sub[billing_col] / nf).replace([np.inf, -np.inf], np.nan)
cv_nb = _cv(sub["normalized_billing"])
ratio = (sub[billing_col] / sub[payment_col]).replace(
[np.inf, -np.inf], np.nan)
volume_col = next(
(col for col in ("total_discharges", "apc_services") if col in sub),
None,
)
total_volume = int(sub[volume_col].sum()) if volume_col else 0
row = OrderedDict(
service_code=code,
service_desc=desc,
dataset=dataset_name,
n_providers=int(_provider_count(sub)),
total_volume=total_volume,
mean_billing=float(round(us_b_mu, 2)),
mean_payment=float(round(us_p_mu, 2)),
median_charge_to_payment=(
float(round(ratio.median(), 3)) if ratio.notna().any() else None
),
cv_b=float(round(cv_b, 4)) if cv_b is not None else None,
cv_p=float(round(cv_p, 4)) if cv_p is not None else None,
cv_nb=float(round(cv_nb, 4)) if cv_nb is not None else None,
)
row.update(metadata)
rows.append(row)
out = pd.DataFrame(rows)
if not out.empty:
out = out.sort_values("cv_b", ascending=False).reset_index(drop=True)
return out
def five_rankings(cv_table: pd.DataFrame) -> dict:
"""Agrawal's 5 ranking lists — one for each 'interesting' axis.
Returns dict: {rank_label: [{service, value, rank}, ...]}.
"""
if cv_table.empty:
return {}
rankings = OrderedDict()
rankings["most_common_by_volume"] = _rank(cv_table, "total_volume", desc=True)
rankings["most_expensive_by_billing"] = _rank(cv_table, "mean_billing", desc=True)
rankings["highest_cv_b"] = _rank(cv_table, "cv_b", desc=True)
rankings["highest_cv_p"] = _rank(cv_table, "cv_p", desc=True)
rankings["highest_cv_nb"] = _rank(cv_table, "cv_nb", desc=True)
return rankings
def _rank(df: pd.DataFrame, col: str, desc: bool, top_n: int = 10) -> list:
sub = df.dropna(subset=[col]).sort_values(col, ascending=not desc).head(top_n)
rows = []
metadata_fields = (
"period", "service_universe", "payment_column", "payment_measure",
"minimum_provider_rule", "comparison_eligible",
"comparison_eligibility",
)
for i, (_, row) in enumerate(sub.iterrows()):
item = OrderedDict(
service_code=row["service_code"],
service_desc=row["service_desc"],
dataset=row["dataset"],
value=float(round(row[col], 4)),
rank=i + 1,
)
for field in metadata_fields:
if field in row.index:
item[field] = row[field]
rows.append(item)
return rows
# ---------------------------------------------------------------------------
# 3b. State normalization factors + B/P/NB maps
# ---------------------------------------------------------------------------
def state_maps(df: pd.DataFrame, service_col: str, service_code: str,
billing_col: str, payment_col: str,
dataset_name: str = "unknown", desc_col: str = None,
period: str = None, service_universe: str = None,
comparison_eligibility: str = None,
comparison_eligible: bool = False,
excluded_states: tuple = ()) -> dict:
"""Per-state B/P/NB maps for one service on the log2 scale.
State cells with fewer than MIN_PROVIDERS_PER_STATE_SERVICE providers are
omitted. The national means remain the means of all valid providers for the
service, while the cell-level eligibility rule is recorded in metadata.
"""
valid = df.dropna(subset=[service_col, billing_col, payment_col]).copy()
comparison_excluded_states = tuple(excluded_states or ())
if comparison_excluded_states and "state" in valid.columns:
valid = valid[~valid["state"].isin(comparison_excluded_states)].copy()
sub = valid[valid[service_col] == service_code]
if _provider_count(sub) < MIN_PROVIDERS_PER_SERVICE:
return {}
us_b_mu = sub[billing_col].mean()
us_p_mu = sub[payment_col].mean()
out = {}
excluded_states = {}
n_states_observed = 0
for state, g in sub.groupby("state"):
n_states_observed += 1
if _provider_count(g) < MIN_PROVIDERS_PER_STATE_SERVICE:
excluded_states[state] = {
"n_providers": int(_provider_count(g)),
"included": False,
"exclusion_reason": "below_min_state_service_n",
"minimum_providers_per_state_service": int(
MIN_PROVIDERS_PER_STATE_SERVICE
),
}
continue
b_mu = g[billing_col].mean()
p_mu = g[payment_col].mean()
b_map = _log2_ratio(b_mu, us_b_mu)
p_map = _log2_ratio(p_mu, us_p_mu)
# NB_s = log2(B_s / P_s) — direct from Agrawal
nb_map = _log2_ratio(b_mu / us_b_mu, p_mu / us_p_mu) if (
b_mu > 0 and us_b_mu > 0 and p_mu > 0 and us_p_mu > 0) else float("nan")
out[state] = {
"n_providers": int(_provider_count(g)),
"mean_billing": float(round(b_mu, 2)),
"mean_payment": float(round(p_mu, 2)),
"b_map": round(b_map, 4) if not pd.isna(b_map) else None,
"p_map": round(p_map, 4) if not pd.isna(p_map) else None,
"nb_map": round(nb_map, 4) if not pd.isna(nb_map) else None,
"minimum_providers_per_state_service": int(
MIN_PROVIDERS_PER_STATE_SERVICE),
"primary_geographic_excluded_states": list(
PRIMARY_GEOGRAPHIC_EXCLUDED_STATES
),
"state_cell_eligible": True,
}
metadata = _analysis_metadata(
df, service_col, billing_col, payment_col, dataset_name,
period=period, service_universe=service_universe,
comparison_eligibility=comparison_eligibility,
comparison_eligible=comparison_eligible,
)
metadata.update({
"analysis": "state-level B/P/NB maps",
"aggregation_unit": "provider means aggregated within state-service cells",
"minimum_providers_per_state_service": int(
MIN_PROVIDERS_PER_STATE_SERVICE),
"state_service_minimum_rule": (
f"at least {MIN_PROVIDERS_PER_STATE_SERVICE} providers per state-service cell"
),
"n_states_observed": int(n_states_observed),
"n_states_eligible": int(len(out)),
"comparison_excluded_states": list(comparison_excluded_states),
})
return {
"service_code": service_code,
"us_mean_billing": float(round(us_b_mu, 2)),
"us_mean_payment": float(round(us_p_mu, 2)),
"n_providers_total": int(_provider_count(sub)),
"states": out,
"excluded_states": excluded_states,
"comparison_excluded_states": list(comparison_excluded_states),
"metadata": metadata,
**metadata,
}
def build_choropleths(df: pd.DataFrame, service_col: str,
billing_col: str, payment_col: str,
service_codes: list, dataset_name: str,
desc_col: str = None, period: str = None,
service_universe: str = None,
comparison_eligibility: str = None,
comparison_eligible: bool = False,
excluded_states: tuple = ()) -> list:
"""Produce choropleth-ready JSON per service."""
out = []
for code in service_codes:
m = state_maps(
df, service_col, code, billing_col, payment_col,
dataset_name=dataset_name, desc_col=desc_col, period=period,
service_universe=service_universe,
comparison_eligibility=comparison_eligibility,
comparison_eligible=comparison_eligible,
excluded_states=excluded_states,
)
if m:
if desc_col and desc_col in df.columns:
descriptions = df.loc[
df[service_col] == code, desc_col
].dropna()
if not descriptions.empty:
m["service_desc"] = str(descriptions.iloc[0])
m["dataset"] = dataset_name
out.append(m)
return out
# ---------------------------------------------------------------------------
# 3c. Billing-payment correlation
# ---------------------------------------------------------------------------
def _correlation_summary(frame: pd.DataFrame, value_col: str,
metadata: dict) -> dict:
"""Summarize one explicitly identified correlation level."""
summary = OrderedDict(
n_services=int(len(frame)),
min_corr=None,
max_corr=None,
median_corr=None,
negative_corr_count=0,
)
if not frame.empty:
values = frame[value_col].dropna()
if not values.empty:
summary.update({
"min_corr": float(round(values.min(), 4)),
"max_corr": float(round(values.max(), 4)),
"median_corr": float(round(values.median(), 4)),
"negative_corr_count": int((values < 0).sum()),
})
if "n_states" in frame.columns and frame["n_states"].notna().any():
summary.update({
"minimum_states_used": int(frame["n_states"].min()),
"services_at_minimum_states": int(
(frame["n_states"] == frame["n_states"].min()).sum()
),
})
summary.update(metadata)
return summary
def _fisher_correlation_interval(correlation: float, n: int) -> tuple[float | None, float | None]:
"""Return a 95% Fisher-z interval for a Pearson correlation."""
if correlation is None or n <= 3 or abs(correlation) >= 1:
return None, None
standard_error = 1 / np.sqrt(n - 3)
center = np.arctanh(np.clip(correlation, -0.999999, 0.999999))
margin = 1.959964 * standard_error
return tuple(
round(float(np.tanh(value)), 4)
for value in (center - margin, center + margin)
)
def correlation_analysis(df: pd.DataFrame, service_col: str,
billing_col: str, payment_col: str,
desc_col: str = None, dataset_name: str = "unknown",
period: str = None, service_universe: str = None,
comparison_eligibility: str = None,
comparison_eligible: bool = False,
excluded_states: tuple = ()) -> dict:
"""Return separate state-ecological and provider-level Pearson results.
The state result is a correlation of provider means after aggregation to
state-service cells. It is therefore ecological, not a provider-level
statistic. State cells below MIN_PROVIDERS_PER_STATE_SERVICE are excluded.
The provider result aggregates duplicate rows to one provider-service mean
before computing Pearson's r.
"""
valid = df.dropna(subset=[service_col, billing_col, payment_col]).copy()
excluded_states = tuple(excluded_states or ())
if excluded_states and "state" in valid.columns:
valid = valid[~valid["state"].isin(excluded_states)].copy()
base_metadata = _analysis_metadata(
df, service_col, billing_col, payment_col, dataset_name,
period=period, service_universe=service_universe,
comparison_eligibility=comparison_eligibility,
comparison_eligible=comparison_eligible,
)
resolved_desc_col = _desc_column(valid, service_col, desc_col)
state_metadata = OrderedDict(base_metadata)
state_metadata.update({
"analysis": "state-level ecological Pearson correlation",
"correlation_type": "state-level ecological Pearson correlation",
"correlation_level": "state-level ecological",
"correlation_method": "Pearson r",
"aggregation_unit": "state-service cells from provider means",
"minimum_providers_per_state_service": int(
MIN_PROVIDERS_PER_STATE_SERVICE),
"state_service_minimum_rule": (
f"at least {MIN_PROVIDERS_PER_STATE_SERVICE} providers per state-service cell"
),
"minimum_state_cells_per_service": int(MIN_STATES_FOR_CORRELATION),
"comparison_excluded_states": list(excluded_states),
"comparison_eligible": False,
"comparison_eligibility": (
"Not directly comparable with a provider-level statistic: this is "
"an ecological state-service aggregation with a broad 2023 service "
"universe and explicit within-state cell threshold"
),
})
provider_metadata = OrderedDict(base_metadata)
provider_metadata.update({
"analysis": "provider-level Pearson correlation",
"correlation_type": "provider-level Pearson correlation",
"correlation_level": "provider-level",
"correlation_method": "Pearson r",
"aggregation_unit": "provider-service observations, provider means",
"comparison_excluded_states": list(excluded_states),
"comparison_eligible": False,
"comparison_eligibility": (
"Provider-level unit is explicit, but period and service universe "
"remain 2023 exploratory rather than a controlled historical match"
),
})
state_records = []
provider_records = []
grouped = valid.groupby(service_col, sort=False)
for code, sub in grouped:
if "provider_ccn" in sub.columns:
provider_tab = sub.groupby("provider_ccn", as_index=False).agg(
b_mu=(billing_col, "mean"),
p_mu=(payment_col, "mean"),
**({"state": ("state", "first")} if "state" in sub.columns else {}),
)
else:
provider_columns = [billing_col, payment_col]
if "state" in sub.columns:
provider_columns.append("state")
provider_tab = sub[provider_columns].rename(
columns={billing_col: "b_mu", payment_col: "p_mu"})
provider_tab["b_mu"] = provider_tab["b_mu"].astype(float)
provider_tab["p_mu"] = provider_tab["p_mu"].astype(float)
n_providers = int(len(provider_tab))
if n_providers < MIN_PROVIDERS_PER_SERVICE:
continue
desc = (sub[resolved_desc_col].dropna().iloc[0]
if resolved_desc_col and sub[resolved_desc_col].notna().any()
else None)
us_b_mu = provider_tab["b_mu"].mean()
us_p_mu = provider_tab["p_mu"].mean()
# Provider-level correlation is kept separate from the ecological result.
provider_corr = provider_tab["b_mu"].corr(
provider_tab["p_mu"], method="pearson")
if pd.notna(provider_corr):
provider_row = OrderedDict(
service_code=code,
service_desc=desc,
n_providers=n_providers,
provider_level_pearson_r=round(float(provider_corr), 4),
)
provider_row.update(provider_metadata)
provider_records.append(provider_row)
# The current map/correlation is a state-level ecological statistic.
if "state" not in provider_tab.columns:
continue
state_all = provider_tab.groupby("state", as_index=False).agg(
b_mu=("b_mu", "mean"),
p_mu=("p_mu", "mean"),
n=("b_mu", "size"),
)
state_tab = state_all[
state_all["n"] >= MIN_PROVIDERS_PER_STATE_SERVICE
].copy()
if len(state_tab) < MIN_STATES_FOR_CORRELATION:
continue
state_corr = state_tab["b_mu"].corr(
state_tab["p_mu"], method="pearson")
if pd.isna(state_corr):
continue
nb = ((state_tab["b_mu"] / us_b_mu) /
(state_tab["p_mu"] / us_p_mu)).replace(
[0, np.inf, -np.inf], np.nan)
nb = nb[(nb > 0) & nb.notna()]
correlation_ci_low, correlation_ci_high = _fisher_correlation_interval(
float(state_corr), len(state_tab)
)
state_row = OrderedDict(
service_code=code,
service_desc=desc,
n_providers=n_providers,
n_states=int(len(state_tab)),
n_states_observed=int(len(state_all)),
n_state_cells_excluded=int(len(state_all) - len(state_tab)),
correlation_b_vs_p=round(float(state_corr), 4),
state_ecological_pearson_r=round(float(state_corr), 4),
correlation_ci_low=correlation_ci_low,
correlation_ci_high=correlation_ci_high,
correlation_interval="95% Fisher-z interval",
median_state_nb=(
round(float(np.median(np.log2(nb))), 4) if not nb.empty else None
),
)
state_row.update(state_metadata)
state_records.append(state_row)
state_df = pd.DataFrame(state_records)
provider_df = pd.DataFrame(provider_records)
if not state_df.empty:
state_df = state_df.sort_values(
"state_ecological_pearson_r").reset_index(drop=True)
if not provider_df.empty:
provider_df = provider_df.sort_values(
"provider_level_pearson_r").reset_index(drop=True)
state_summary = _correlation_summary(
state_df, "state_ecological_pearson_r", state_metadata)
provider_summary = _correlation_summary(
provider_df, "provider_level_pearson_r", provider_metadata)
return {
"metadata": {
"state_ecological": state_metadata,
"provider_level": provider_metadata,
},
# Keep per_service as the existing state-level site entry point.
"per_service": state_df.to_dict(orient="records"),
"state_ecological_per_service": state_df.to_dict(orient="records"),
"provider_level_per_service": provider_df.to_dict(orient="records"),
"summary": state_summary,
"provider_level_summary": provider_summary,
"lowest_5": state_df.head(5).to_dict(orient="records"),
"highest_5": state_df.tail(5).to_dict(orient="records"),
"provider_level_lowest_5": provider_df.head(5).to_dict(orient="records"),
"provider_level_highest_5": provider_df.tail(5).to_dict(orient="records"),
"historical_comparison": (
"No direct historical benchmark is reported: aggregation level, "
"service universe, payment field, and period must be aligned before "
"a numerical comparison is valid."
),
}
# ---------------------------------------------------------------------------
# 3d. Cross-dataset provider top-quartile comparison
# ---------------------------------------------------------------------------
def _provider_proxy_table(df: pd.DataFrame, service_col: str,
billing_col: str, payment_col: str,
prefix: str) -> tuple[pd.DataFrame, dict]:
"""Aggregate an unweighted provider mean of valid charge/payment ratios."""
required = ["provider_ccn", service_col, billing_col, payment_col]
missing = [col for col in required if col not in df.columns]
if missing:
raise KeyError(f"missing columns for provider proxy: {missing}")
universe_count = int(df[service_col].nunique())
valid = df.dropna(subset=required).copy()
valid["_charge_payment_proxy"] = (
valid[billing_col] / valid[payment_col]
).replace([np.inf, -np.inf], np.nan)
valid = valid[
(valid[payment_col] > 0)
& (valid[billing_col] > 0)
& valid["_charge_payment_proxy"].notna()
].copy()
metadata = {
"service_column": service_col,
"service_universe": f"all observed {service_col} codes in the supplied frame",
"service_universe_count": universe_count,
"eligible_service_universe_count": _eligible_service_count(
df, service_col, billing_col, payment_col),
"billing_column": billing_col,
"payment_column": payment_col,
"payment_measure": _payment_measure(payment_col),
"ratio_definition": (
f"row-level {billing_col} / {payment_col}, then unweighted mean "
"across each provider's valid service rows"
),
"observed_service_rows": int(len(df)),
"valid_payment_rows": int(len(valid)),
"providers_observed": int(df["provider_ccn"].nunique()),
"providers_with_valid_proxy": int(valid["provider_ccn"].nunique()),
}
if valid.empty:
return pd.DataFrame(columns=["provider_ccn"]), metadata
aggregations = {
f"{prefix}_mean_charge_to_payment_proxy": (
"_charge_payment_proxy", "mean"),
f"{prefix}_median_charge_to_payment_proxy": (
"_charge_payment_proxy", "median"),
f"{prefix}_service_breadth": (service_col, "nunique"),
f"{prefix}_service_opportunities": (
"_charge_payment_proxy", "size"),
}
for column in ("provider_name", "state", "census_region", "urban_rural"):
if column in valid.columns:
aggregations[column] = (column, "first")
provider = valid.groupby("provider_ccn", as_index=False).agg(**aggregations)
breadth_col = f"{prefix}_service_breadth"
opportunity_col = f"{prefix}_service_opportunities"
provider[f"{prefix}_service_breadth_pct"] = (
100 * provider[breadth_col] / universe_count if universe_count else 0.0
).round(2)
provider[f"{prefix}_payment_column"] = payment_col
provider[f"{prefix}_payment_measure"] = _payment_measure(payment_col)
provider[f"{prefix}_payment_denominator"] = payment_col
provider[f"{prefix}_billing_column"] = billing_col
provider[f"{prefix}_opportunity_definition"] = (
"valid observed provider-service rows with a positive payment denominator"
)
provider[f"{prefix}_service_universe_count"] = universe_count
provider[f"{prefix}_eligible_service_universe_count"] = metadata[
"eligible_service_universe_count"]
provider[f"{prefix}_valid_payment_rows"] = provider[opportunity_col]
return provider, metadata
def cross_dataset_top_quartile_providers(
inp_df: pd.DataFrame, out_df: pd.DataFrame,
inp_billing_col: str = "avg_submitted_charge",
inp_payment_col: str = "avg_medicare_payment",
out_billing_col: str = "avg_submitted_charge",
out_payment_col: str = "avg_allowed_amount",
inp_service_col: str = "drg_code",
out_service_col: str = "apc_code",
inp_period: str = INPATIENT_PERIOD,
out_period: str = OUTPATIENT_PERIOD) -> dict:
"""Compare providers in both datasets using explicitly defined proxies.
A selected provider is in the top quartile of the provider-level mean
charge-to-payment proxy in both datasets. This is a descriptive cohort, not
a measure of service-level CV or an explanation of causes.
"""
inp_p, inp_meta = _provider_proxy_table(
inp_df, inp_service_col, inp_billing_col, inp_payment_col, "inp")
out_p, out_meta = _provider_proxy_table(
out_df, out_service_col, out_billing_col, out_payment_col, "out")
observed_inp = set(inp_df["provider_ccn"].dropna())
observed_out = set(out_df["provider_ccn"].dropna())
observed_both = observed_inp & observed_out
merged = inp_p.merge(out_p, on="provider_ccn", how="inner",
suffixes=("", "_out"))
for column in ("provider_name", "state", "census_region", "urban_rural"):
out_column = f"{column}_out"
if column in merged.columns and out_column in merged.columns:
merged[column] = merged[column].fillna(merged[out_column])
merged = merged.drop(columns=[out_column])
elif column not in merged.columns and out_column in merged.columns:
merged = merged.rename(columns={out_column: column})
inp_proxy = "inp_mean_charge_to_payment_proxy"
out_proxy = "out_mean_charge_to_payment_proxy"
if merged.empty:
inp_q3 = out_q3 = None
merged["inpatient_top_quartile"] = pd.Series(dtype=bool)
merged["outpatient_top_quartile"] = pd.Series(dtype=bool)
merged["top_quartile_in_both"] = pd.Series(dtype=bool)
else:
inp_q3 = merged[inp_proxy].quantile(0.75)
out_q3 = merged[out_proxy].quantile(0.75)
merged["inpatient_top_quartile"] = merged[inp_proxy] >= inp_q3
merged["outpatient_top_quartile"] = merged[out_proxy] >= out_q3
merged["top_quartile_in_both"] = (
merged["inpatient_top_quartile"]
& merged["outpatient_top_quartile"]
)
comparison_denominator = int(len(merged))
merged["comparison_denominator"] = comparison_denominator
merged["comparison_eligibility"] = (
"Providers with valid provider-level mean proxies in both datasets; "
"descriptive cross-period comparison only"
)
merged["inpatient_period"] = inp_period
merged["outpatient_period"] = out_period
selected = merged[merged["top_quartile_in_both"]].sort_values(
by=[inp_proxy, out_proxy], ascending=False)
pct = lambda count: (
float(round(100 * count / comparison_denominator, 2))
if comparison_denominator else None
)
summary = OrderedDict(
providers_in_both=int(len(observed_both)),
providers_with_valid_proxy_in_both=comparison_denominator,
comparison_denominator=comparison_denominator,
provider_comparison_denominator=comparison_denominator,
inpatient_provider_denominator=int(inp_meta["providers_with_valid_proxy"]),
outpatient_provider_denominator=int(out_meta["providers_with_valid_proxy"]),
opportunity_denominator=comparison_denominator,
inpatient_opportunity_rows=int(inp_meta["valid_payment_rows"]),
outpatient_opportunity_rows=int(out_meta["valid_payment_rows"]),
inpatient_q3_threshold=(
float(round(inp_q3, 3)) if inp_q3 is not None else None
),
outpatient_q3_threshold=(
float(round(out_q3, 3)) if out_q3 is not None else None
),
inpatient_top_quartile_count=int(
merged["inpatient_top_quartile"].sum()),
outpatient_top_quartile_count=int(
merged["outpatient_top_quartile"].sum()),
top_quartile_in_both_count=int(
merged["top_quartile_in_both"].sum()),
inpatient_top_quartile_pct=pct(
int(merged["inpatient_top_quartile"].sum())),
outpatient_top_quartile_pct=pct(
int(merged["outpatient_top_quartile"].sum())),
top_quartile_in_both_pct=pct(
int(merged["top_quartile_in_both"].sum())),
)
metadata = {
"analysis": "providers in the top quartile of provider-level mean charge-to-payment proxy in both datasets",
"period": {"inpatient": inp_period, "outpatient": out_period},
"periods": {"inpatient": inp_period, "outpatient": out_period},
"service_universe": {
"inpatient": inp_meta["service_universe"],
"outpatient": out_meta["service_universe"],
},
"service_universe_counts": {
"inpatient": inp_meta["service_universe_count"],
"outpatient": out_meta["service_universe_count"],
},
"eligible_service_universe_counts": {
"inpatient": inp_meta["eligible_service_universe_count"],
"outpatient": out_meta["eligible_service_universe_count"],
},
"payment_denominator": {
"inpatient": {
"column": inp_payment_col,
"measure": _payment_measure(inp_payment_col),
"ratio": f"{inp_billing_col} / {inp_payment_col}",
},
"outpatient": {
"column": out_payment_col,
"measure": _payment_measure(out_payment_col),
"ratio": f"{out_billing_col} / {out_payment_col}",
},
},
"provider_aggregation": (
"unweighted mean of valid row-level charge/payment ratios; no service "
"standardization or causal interpretation is implied"
),
"opportunity_definition": (
"each observed provider-service row with a positive payment denominator"
),
"comparison_denominator": comparison_denominator,
"inpatient_opportunity_rows": inp_meta["valid_payment_rows"],
"outpatient_opportunity_rows": out_meta["valid_payment_rows"],
"comparison_eligibility": (
"Descriptive cross-period comparison among providers with valid proxies "
"in both datasets; not a causal or service-standardized comparison"
),
"comparison_eligible": False,
"minimum_provider_rule": (
f"service-level eligibility elsewhere uses at least {MIN_PROVIDERS_PER_SERVICE} "
"providers; provider proxy means retain valid observed opportunities"
),
"service_breadth_fields": {
"inpatient": "inp_service_breadth",
"outpatient": "out_service_breadth",
},
"opportunity_fields": {
"inpatient": "inp_service_opportunities",
"outpatient": "out_service_opportunities",
},
}
scatter_columns = [
"provider_ccn", "provider_name", "state", "census_region",
"urban_rural", inp_proxy, "inp_median_charge_to_payment_proxy",
"inp_service_breadth", "inp_service_breadth_pct",
"inp_service_opportunities", "inp_service_universe_count",
"inp_eligible_service_universe_count", "inp_payment_column",
"inp_payment_measure", "inp_payment_denominator", out_proxy,
"out_median_charge_to_payment_proxy", "out_service_breadth",
"out_service_breadth_pct", "out_service_opportunities",
"out_service_universe_count", "out_eligible_service_universe_count",
"out_payment_column", "out_payment_measure", "out_payment_denominator",
"inpatient_top_quartile", "outpatient_top_quartile",
"top_quartile_in_both", "comparison_denominator",
"comparison_eligibility", "inpatient_period", "outpatient_period",
]
scatter_columns = [col for col in scatter_columns if col in merged.columns]
top_columns = scatter_columns
return {
"metadata": metadata,
"summary": summary,
"scatter_data": merged[scatter_columns].to_dict(orient="records"),
"top_quartile_providers": selected[top_columns].head(50).to_dict(
orient="records"),
}
# ---------------------------------------------------------------------------
# Orchestration
# ---------------------------------------------------------------------------
def save_json(obj, name: str):
SITE_DATA.mkdir(parents=True, exist_ok=True)
path = SITE_DATA / name
with open(path, "w") as f:
json.dump(obj, f, indent=2, default=_json_default)
print(f" wrote {path.name} ({path.stat().st_size:,} bytes)")
def _json_default(o):
if isinstance(o, (np.integer,)):
return int(o)
if isinstance(o, (np.floating,)):
v = float(o)
return None if pd.isna(v) else v
if isinstance(o, (np.ndarray,)):
return o.tolist()
if pd.isna(o):
return None
raise TypeError(f"not serializable: {type(o)}")
def run() -> dict:
"""Main entrypoint — reads parquet snapshots produced by run_pipeline.py."""
print("[Phase 3] Loading feature parquets ...")
inp = pd.read_parquet(DATA_PROC / "inpatient_features.parquet")
out = pd.read_parquet(DATA_PROC / "outpatient_features.parquet")
print(f" inpatient: {len(inp):,} rows | outpatient: {len(out):,} rows")
artifacts = {}
# 3a — 3-tier CV tables
print("[Phase 3] Tier-3 CV tables ...")
inp_payment_col = "avg_medicare_payment"
inp_total_payment_col = "avg_total_payment"
out_payment_col = "avg_allowed_amount"
inp_service_universe = (
"all observed inpatient DRG codes in the cleaned FY2023 feature frame"
)
out_service_universe = (
"all observed outpatient APC codes with cost-observed rows in the cleaned CY2023 feature frame"
)
inp_comparable_note = (
"Primary inpatient track uses the average Medicare payment field. The "
"FY2023 period, broad service universe, and >=30-provider rule are "
"recorded explicitly rather than treated as a controlled historical match"
)
inp_total_note = (
"Sensitivity track using average total payment; not the primary "
"Medicare-payment comparable track"
)
out_exploratory_note = (
"Outpatient extension using average Medicare allowed amount; CY2023 and "
"cost-observed suppression universe are explicit"
)
inp_cv = compute_cv_table(inp, "drg_code", "drg_desc",
"avg_submitted_charge", inp_payment_col, "inpatient",
period=INPATIENT_PERIOD,
service_universe=inp_service_universe,
comparison_eligibility=inp_comparable_note,
comparison_eligible=False)
inp_total_cv = compute_cv_table(
inp, "drg_code", "drg_desc", "avg_submitted_charge",
inp_total_payment_col, "inpatient_total_payment",
period=INPATIENT_PERIOD, service_universe=inp_service_universe,
comparison_eligibility=inp_total_note, comparison_eligible=False)
out_cv = compute_cv_table(out, "apc_code", "apc_desc",
"avg_submitted_charge", out_payment_col, "outpatient",
period=OUTPATIENT_PERIOD,
service_universe=out_service_universe,
comparison_eligibility=out_exploratory_note,
comparison_eligible=False)
print(f" inpatient: {len(inp_cv)} DRGs ≥30 providers")
print(f" inpatient total-payment sensitivity: {len(inp_total_cv)} DRGs")
print(f" outpatient: {len(out_cv)} APCs ≥30 providers")
# Five rankings
inp_ranks = five_rankings(inp_cv)
inp_total_ranks = five_rankings(inp_total_cv)
out_ranks = five_rankings(out_cv)
artifacts["cv_tables"] = {
"inpatient": inp_cv.to_dict(orient="records"),
"inpatient_total_payment": inp_total_cv.to_dict(orient="records"),
"outpatient": out_cv.to_dict(orient="records"),
"metadata": {
"inpatient": _analysis_metadata(
inp, "drg_code", "avg_submitted_charge", inp_payment_col,
"inpatient", period=INPATIENT_PERIOD,
service_universe=inp_service_universe,
comparison_eligibility=inp_comparable_note),
"inpatient_total_payment": _analysis_metadata(
inp, "drg_code", "avg_submitted_charge", inp_total_payment_col,
"inpatient_total_payment", period=INPATIENT_PERIOD,
service_universe=inp_service_universe,
comparison_eligibility=inp_total_note),
"outpatient": _analysis_metadata(
out, "apc_code", "avg_submitted_charge", out_payment_col,
"outpatient", period=OUTPATIENT_PERIOD,
service_universe=out_service_universe,
comparison_eligibility=out_exploratory_note),
},
}
artifacts["five_rankings"] = {
"inpatient": inp_ranks,
"inpatient_total_payment": inp_total_ranks,
"outpatient": out_ranks,
"metadata": {
"inpatient": artifacts["cv_tables"]["metadata"]["inpatient"],
"inpatient_total_payment": artifacts["cv_tables"]["metadata"][
"inpatient_total_payment"],
"outpatient": artifacts["cv_tables"]["metadata"]["outpatient"],
},
}
# 3b — choropleths for top services from each ranking
print("[Phase 3] State choropleths (B/P/NB maps) ...")
# Take top-N from each ranking (deduplicated)
top_codes_inp = set()
for r in inp_ranks.values():
for s in r:
top_codes_inp.add(s["service_code"])
top_codes_out = set()
for r in out_ranks.values():
for s in r:
top_codes_out.add(s["service_code"])
inp_choropleths = build_choropleths(
inp, "drg_code", "avg_submitted_charge", inp_payment_col,
sorted(top_codes_inp), "inpatient", "drg_desc",
period=INPATIENT_PERIOD, service_universe=inp_service_universe,
comparison_eligibility=inp_comparable_note,
excluded_states=PRIMARY_GEOGRAPHIC_EXCLUDED_STATES)
inp_total_choropleths = build_choropleths(
inp, "drg_code", "avg_submitted_charge", inp_total_payment_col,
sorted(top_codes_inp), "inpatient_total_payment", "drg_desc",
period=INPATIENT_PERIOD, service_universe=inp_service_universe,
comparison_eligibility=inp_total_note,
excluded_states=PRIMARY_GEOGRAPHIC_EXCLUDED_STATES)
out_choropleths = build_choropleths(
out, "apc_code", "avg_submitted_charge", out_payment_col,
sorted(top_codes_out), "outpatient", "apc_desc",
period=OUTPATIENT_PERIOD, service_universe=out_service_universe,
comparison_eligibility=out_exploratory_note)
print(f" inpatient choropleths: {len(inp_choropleths)} top DRGs")
print(f" inpatient total-payment sensitivity choropleths: "
f"{len(inp_total_choropleths)} top DRGs")
print(f" outpatient choropleths: {len(out_choropleths)} top APCs")
artifacts["choropleths"] = {
"inpatient": inp_choropleths,
"inpatient_total_payment": inp_total_choropleths,
"outpatient": out_choropleths,
"metadata": {
"inpatient": artifacts["cv_tables"]["metadata"]["inpatient"],
"inpatient_total_payment": artifacts["cv_tables"]["metadata"][
"inpatient_total_payment"],
"outpatient": artifacts["cv_tables"]["metadata"]["outpatient"],
"minimum_providers_per_state_service": int(
MIN_PROVIDERS_PER_STATE_SERVICE),
},
}
# 3c — billing-payment correlation
print("[Phase 3] Billing-payment correlation analysis ...")
inp_corr = correlation_analysis(inp, "drg_code",
"avg_submitted_charge", inp_payment_col,
desc_col="drg_desc", dataset_name="inpatient",
period=INPATIENT_PERIOD,
service_universe=inp_service_universe,
comparison_eligibility=inp_comparable_note,
excluded_states=PRIMARY_GEOGRAPHIC_EXCLUDED_STATES)
inp_total_corr = correlation_analysis(
inp, "drg_code", "avg_submitted_charge", inp_total_payment_col,
desc_col="drg_desc", dataset_name="inpatient_total_payment",
period=INPATIENT_PERIOD, service_universe=inp_service_universe,
comparison_eligibility=inp_total_note,
excluded_states=PRIMARY_GEOGRAPHIC_EXCLUDED_STATES)
out_corr = correlation_analysis(out, "apc_code",
"avg_submitted_charge", out_payment_col,
desc_col="apc_desc", dataset_name="outpatient",
period=OUTPATIENT_PERIOD,
service_universe=out_service_universe,
comparison_eligibility=out_exploratory_note,
excluded_states=PRIMARY_GEOGRAPHIC_EXCLUDED_STATES)
inp_corr_including_md = correlation_analysis(
inp, "drg_code", "avg_submitted_charge", inp_payment_col,
desc_col="drg_desc", dataset_name="inpatient_including_md",
period=INPATIENT_PERIOD, service_universe=inp_service_universe,
comparison_eligibility="Sensitivity including Maryland; not the primary geographic summary",
)
inp_total_corr_including_md = correlation_analysis(
inp, "drg_code", "avg_submitted_charge", inp_total_payment_col,
desc_col="drg_desc", dataset_name="inpatient_total_payment_including_md",
period=INPATIENT_PERIOD, service_universe=inp_service_universe,
comparison_eligibility="Sensitivity including Maryland; not the primary geographic summary",
)
print(f" inpatient: {inp_corr.get('summary', {}).get('n_services', 0)} services, "
f"median corr = {inp_corr.get('summary', {}).get('median_corr', 'n/a')}")
print(f" inpatient total-payment sensitivity: "
f"{inp_total_corr.get('summary', {}).get('n_services', 0)} services")
print(f" outpatient: {out_corr.get('summary', {}).get('n_services', 0)} services, "
f"median corr = {out_corr.get('summary', {}).get('median_corr', 'n/a')}")
artifacts["correlations"] = {
"inpatient": inp_corr,
"inpatient_total_payment": inp_total_corr,
"outpatient": out_corr,
"metadata": {
"inpatient": inp_corr["metadata"],
"inpatient_total_payment": inp_total_corr["metadata"],
"outpatient": out_corr["metadata"],
"primary_geographic_excluded_states": list(
PRIMARY_GEOGRAPHIC_EXCLUDED_STATES
),
},
}
artifacts["maryland_sensitivity"] = {
"primary_excluded_states": list(PRIMARY_GEOGRAPHIC_EXCLUDED_STATES),
"note": (
"Primary inpatient maps and correlations exclude Maryland because its "
"all-payer rate-setting system is not comparable to the standard Medicare context."
),
"inpatient": {
"excluded": inp_corr["summary"],
"included": inp_corr_including_md["summary"],
"service_comparison": [
{
"service_code": row.get("service_code"),
"correlation_excluded_md": row.get("state_ecological_pearson_r"),
}
for row in inp_corr.get("state_ecological_per_service", [])
],
},
"inpatient_total_payment": {
"excluded": inp_total_corr["summary"],
"included": inp_total_corr_including_md["summary"],
},
}
# 3d — cross-dataset provider top-quartile comparison
print("[Phase 3] Cross-dataset provider top-quartile comparison ...")
top_quartile_providers = cross_dataset_top_quartile_providers(
inp, out,
inp_billing_col="avg_submitted_charge",
inp_payment_col=inp_payment_col,
out_billing_col="avg_submitted_charge",
out_payment_col=out_payment_col,
inp_service_col="drg_code", out_service_col="apc_code",
inp_period=INPATIENT_PERIOD, out_period=OUTPATIENT_PERIOD)
print(
f" valid provider denominator: "
f"{top_quartile_providers['summary']['comparison_denominator']:,}, "
f"top quartile in both: "
f"{top_quartile_providers['summary']['top_quartile_in_both_count']:,} "
f"({top_quartile_providers['summary']['top_quartile_in_both_pct']}%)"
)
artifacts["cross_dataset_top_quartile_providers"] = top_quartile_providers
# Persist JSON for site consumer
print("[Phase 3] Serializing to site/data/ ...")
save_json({
"inpatient": artifacts["cv_tables"]["inpatient"],
"inpatient_total_payment": artifacts["cv_tables"][
"inpatient_total_payment"],
"outpatient": artifacts["cv_tables"]["outpatient"],
"metadata": artifacts["cv_tables"]["metadata"],
}, "cv_tables.json")
save_json(artifacts["five_rankings"], "five_rankings.json")
save_json({
"inpatient": inp_choropleths,
"inpatient_total_payment": inp_total_choropleths,
"outpatient": out_choropleths,
"metadata": artifacts["choropleths"]["metadata"],
}, "regional_choropleths.json")
save_json(artifacts["correlations"], "billing_payment_correlation.json")
save_json(artifacts["maryland_sensitivity"], "maryland_sensitivity.json")
save_json(top_quartile_providers, "cross_dataset_variators.json")
print("[Phase 3] DONE.")
return artifacts
if __name__ == "__main__": # pragma: no cover
run()