"""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()