analysis_outliers.py

1270 lines · 56,485 bytes · view raw

"""Phase 4 — Project-defined high-charge flag analysis.

The primary cohort is defined within each service as providers at or above the
75th percentile of submitted charge.  IQR and z-score screens provide separate
statistical comparisons.  The outputs distinguish provider-service flag counts
from payment-dollar measures and include explicit eligibility denominators for
the cross-dataset provider overlap.
"""
from collections import OrderedDict

import numpy as np
import pandas as pd

from .config import (
    DATA_PROC, MIN_PROVIDERS_PER_SERVICE,
    IQR_MULTIPLIER, ZSCORE_MODERATE, ZSCORE_EXTREME,
)

from .analysis_insights import save_json


HIGH_CHARGE_METHOD = "project_high_charge_top_quartile"
HIGH_CHARGE_QUANTILE = 0.75
CONCENTRATION_TARGET_SHARE = 0.40
FLAG_PREVIEW_LIMIT = 2000
COHORT_PREVIEW_LIMIT = 50
BREADTH_BIN_LABELS = ("1", "2-4", "5-9", "10+")
DATASET_PERIODS = {"inpatient": "FY2023", "outpatient": "CY2023"}


def _period_from_frames(*frames) -> str:
    """Read a period label when frames carry one; use the project period otherwise."""
    period_columns = ("period", "data_period", "year", "data_year", "reporting_year")
    for frame in frames:
        if not isinstance(frame, pd.DataFrame):
            continue
        for column in period_columns:
            if column not in frame.columns:
                continue
            values = frame[column].dropna().astype(str).unique()
            if len(values) == 1:
                return values[0]
    return "2023"


def _analysis_metadata(period: str = "dataset-specific") -> OrderedDict:
    """Return shared method metadata for every Phase 4 output."""
    return OrderedDict(
        method=HIGH_CHARGE_METHOD,
        method_label="Project-defined high-charge cohort",
        inspired_by=(
            "OIG top-quartile screening approach; no external source data is "
            "required for this project method"
        ),
        flag_label="high-charge flag",
        flag_plural="high-charge flags",
        flag_definition=(
            "A high-charge flag is a provider-service row at or above the "
            "within-service 75th percentile of avg_submitted_charge."
        ),
        thresholds=OrderedDict(
            min_providers_per_service=int(MIN_PROVIDERS_PER_SERVICE),
            high_charge_quantile=HIGH_CHARGE_QUANTILE,
            iqr_multiplier=float(IQR_MULTIPLIER),
            zscore_moderate=float(ZSCORE_MODERATE),
            zscore_extreme=float(ZSCORE_EXTREME),
            concentration_target_share=CONCENTRATION_TARGET_SHARE,
        ),
        selection_metric="avg_submitted_charge",
        proxy_measure=OrderedDict(
            inpatient="avg_submitted_charge / avg_total_payment",
            outpatient="avg_submitted_charge / avg_allowed_amount",
        ),
        proxy_selection_note=(
            "The proxy reuses submitted charge as its numerator after the cohort "
            "is selected on submitted charge; proxy premiums are post-selection descriptions."
        ),
        payment_measure=OrderedDict(
            inpatient="avg_total_payment",
            outpatient="avg_allowed_amount",
        ),
        period=str(period),
        periods=OrderedDict(DATASET_PERIODS),
        period_note=(
            "Inpatient is FY2023 and outpatient is CY2023; cross-dataset "
            "results are descriptive rather than same-period or causal comparisons."
        ),
    )


def _stable_sort(df: pd.DataFrame, columns: list) -> pd.DataFrame:
    """Sort records deterministically without changing their values."""
    if df.empty:
        return df.reset_index(drop=True)
    sort_columns = [column for column in columns if column in df.columns]
    if not sort_columns:
        return df.reset_index(drop=True)
    return df.sort_values(sort_columns, kind="mergesort", na_position="last").reset_index(drop=True)


def _rounded(value, digits: int):
    if value is None or pd.isna(value):
        return None
    return float(round(float(value), digits))


def _mean(series: pd.Series, digits: int):
    if series is None or series.dropna().empty:
        return None
    return _rounded(series.mean(), digits)


def _premium_pct(high_value, low_value):
    if high_value is None or low_value is None or pd.isna(high_value) or pd.isna(low_value):
        return None
    if low_value == 0:
        return None
    return _rounded(100 * (high_value / low_value - 1), 1)


def _rate(numerator: int, denominator: int):
    if not denominator:
        return 0.0
    return float(numerator / denominator)


def _breadth_rate_profile(
    provider_ids: set,
    flagged_providers: set,
    opportunity_counts: dict,
    n_bins: int = 5,
) -> tuple[dict, dict]:
    """Estimate flag rates by provider opportunity-breadth bin."""
    values = pd.Series(
        {
            provider: float(opportunity_counts.get(provider, 0))
            for provider in sorted(provider_ids)
        },
        dtype=float,
    )
    if values.empty:
        return {}, {}
    rank = values.rank(method="first")
    bin_count = min(n_bins, len(values))
    bins = pd.qcut(rank, q=bin_count, labels=False, duplicates="drop")
    assignments = {provider: int(bin_id) for provider, bin_id in bins.items()}
    rates = {}
    details = {}
    for bin_id in sorted(set(assignments.values())):
        members = [provider for provider, value in assignments.items() if value == bin_id]
        flagged_count = sum(provider in flagged_providers for provider in members)
        rates[bin_id] = _rate(flagged_count, len(members))
        details[str(bin_id)] = OrderedDict(
            provider_count=int(len(members)),
            flagged_provider_count=int(flagged_count),
            flag_rate=float(rates[bin_id]),
            mean_eligible_service_opportunities=float(
                round(values.loc[members].mean(), 2)
            ),
        )
    return assignments, details


def _breadth_matched_expected_overlap(
    shared_providers: set,
    inp_flagged: set,
    out_flagged: set,
    inp_opportunities: dict,
    out_opportunities: dict,
) -> tuple[float, dict]:
    """Compute an independence null while matching opportunity breadth."""
    inp_bins, inp_details = _breadth_rate_profile(
        shared_providers, inp_flagged, inp_opportunities
    )
    out_bins, out_details = _breadth_rate_profile(
        shared_providers, out_flagged, out_opportunities
    )
    expected = sum(
        inp_details[str(inp_bins[provider])]["flag_rate"]
        * out_details[str(out_bins[provider])]["flag_rate"]
        for provider in shared_providers
    )
    return float(expected), OrderedDict(
        benchmark="independence conditional on shared provider opportunity-breadth quintiles",
        n_bins=max(len(inp_details), len(out_details)),
        inpatient_bins=inp_details,
        outpatient_bins=out_details,
        expected_flagged_provider_overlap=float(round(expected, 4)),
    )


def _flagged_rows(flags_df: pd.DataFrame, flag_column: str = "is_high_charge_flag") -> pd.DataFrame:
    """Return only rows carrying a true high-charge flag."""
    if not isinstance(flags_df, pd.DataFrame) or flags_df.empty:
        return flags_df.copy() if isinstance(flags_df, pd.DataFrame) else pd.DataFrame()
    if flag_column in flags_df.columns:
        return flags_df[flags_df[flag_column].fillna(False)].copy()
    # This fallback keeps callers that already pass a flagged-only table useful.
    return flags_df.copy()


def _records(df: pd.DataFrame, limit: int = None) -> list:
    """Convert a stable, optionally truncated frame to JSON records."""
    if not isinstance(df, pd.DataFrame) or df.empty:
        return []
    if limit is not None:
        df = df.head(limit)
    return df.to_dict(orient="records")

# ---------------------------------------------------------------------------
# Helpers
# ---------------------------------------------------------------------------

def _per_service(df: pd.DataFrame, service_col: str, billing_col: str,
                 payment_col: str) -> list:
    """Return list of (service_code, sub_df) for services with ≥30 providers."""
    out = []
    for code, sub in df.groupby(service_col, as_index=False, sort=True):
        if len(sub) >= MIN_PROVIDERS_PER_SERVICE:
            out.append((code, sub))
    return out


# ---------------------------------------------------------------------------
# 4a. Project-defined high-charge cohort — top quartile per service
# ---------------------------------------------------------------------------

def high_charge_cohort_method(df: pd.DataFrame, service_col: str, desc_col: str,
                              billing_col: str, payment_col: str,
                              volume_col: str, dataset_name: str) -> dict:
    """Flag high-charge provider-service rows and compare the two cohorts."""
    flags = []
    cohort_rows = []
    for code, sub in _per_service(df, service_col, billing_col, payment_col):
        q3 = sub[billing_col].quantile(HIGH_CHARGE_QUANTILE)
        if pd.isna(q3):
            continue
        desc = sub[desc_col].iloc[0]
        sub = sub.copy()
        sub["is_high_charge_flag"] = sub[billing_col] >= q3
        sub["outlier_method"] = HIGH_CHARGE_METHOD
        sub["service_code"] = code
        sub["service_desc"] = desc
        sub["eligible_service"] = True
        sub["eligible_provider_count"] = int(sub["provider_ccn"].nunique())
        sub["high_charge_threshold"] = q3
        proxy = (
            sub["charge_to_payment_ratio"]
            if "charge_to_payment_ratio" in sub.columns
            else sub[billing_col].div(sub[payment_col]).replace([np.inf, -np.inf], np.nan)
        )
        sub["charge_to_payment_proxy"] = proxy
        flag_columns = [
            "provider_ccn", "provider_name", "state", "census_region",
            "urban_rural", "service_code", "service_desc",
            "is_high_charge_flag", "outlier_method", "eligible_service",
            "eligible_provider_count", "high_charge_threshold",
            billing_col, payment_col, volume_col, "charge_to_payment_proxy",
        ]
        flags.append(sub[flag_columns].rename(columns={
            billing_col: "submitted_charge",
            payment_col: "payment",
            volume_col: "volume",
        }))
        # Cohort aggregate row
        hi = sub[sub["is_high_charge_flag"]]
        lo = sub[~sub["is_high_charge_flag"]]
        if len(hi) == 0 or len(lo) == 0:
            continue
        high_proxy = _mean(hi["charge_to_payment_proxy"], 3)
        low_proxy = _mean(lo["charge_to_payment_proxy"], 3)
        cohort_rows.append(OrderedDict(
            dataset=dataset_name,
            period=DATASET_PERIODS.get(dataset_name, _period_from_frames(df)),
            service_code=code,
            service_desc=desc,
            n_providers=int(len(sub)),
            n_high_charge_flags=int(len(hi)),
            n_high_outliers=int(len(hi)),
            high_mean_charge=_mean(hi[billing_col], 2),
            low_mean_charge=_mean(lo[billing_col], 2),
            charge_premium_pct=_premium_pct(hi[billing_col].mean(), lo[billing_col].mean()),
            high_mean_charge_to_payment_proxy=high_proxy,
            low_mean_charge_to_payment_proxy=low_proxy,
            charge_to_payment_proxy_premium_pct=_premium_pct(high_proxy, low_proxy),
            proxy_selection_note=(
                "Post-selection descriptive premium: cohort selected on submitted charge "
                "and proxy numerator is submitted charge."
            ),
            high_mean_volume=_mean(hi[volume_col], 1),
            low_mean_volume=_mean(lo[volume_col], 1),
            volume_ratio=_rounded(hi[volume_col].mean() / lo[volume_col].mean(), 3)
                if lo[volume_col].mean() > 0 else None,
            selection_metric=billing_col,
            payment_measure=payment_col,
            high_charge_quantile=HIGH_CHARGE_QUANTILE,
            high_charge_threshold=_rounded(q3, 2),
        ))
    flags_df = pd.concat(flags, ignore_index=True) if flags else pd.DataFrame()
    flags_df = _stable_sort(flags_df, ["service_code", "provider_ccn"])
    cohort_df = _stable_sort(pd.DataFrame(cohort_rows), ["service_code", "service_desc"])
    return {
        "metadata": _analysis_metadata(
            DATASET_PERIODS.get(dataset_name, _period_from_frames(df))
        ),
        "flags": flags_df,
        "cohort_comparison": cohort_df,
    }


# ---------------------------------------------------------------------------
# 4b. Statistical methods — IQR + z-score (per service)
# ---------------------------------------------------------------------------

def iqr_method(df: pd.DataFrame, service_col: str, billing_col: str) -> pd.DataFrame:
    """Per service, flag providers outside [Q1 - 1.5*IQR, Q3 + 1.5*IQR] on billing."""
    flags = []
    for code, sub in _per_service(df, service_col, billing_col, billing_col):
        q1 = sub[billing_col].quantile(0.25)
        q3 = sub[billing_col].quantile(0.75)
        iqr = q3 - q1
        lo_thr = q1 - IQR_MULTIPLIER * iqr
        hi_thr = q3 + IQR_MULTIPLIER * iqr
        sub = sub.copy()
        sub["is_iqr_outlier"] = (sub[billing_col] < lo_thr) | (sub[billing_col] > hi_thr)
        sub["iqr_direction"] = np.where(sub[billing_col] > hi_thr, "high",
                                        np.where(sub[billing_col] < lo_thr, "low", "none"))
        flags.append(sub[sub["is_iqr_outlier"]][[
            "provider_ccn", "provider_name", "state", "census_region", "urban_rural",
            service_col, billing_col, "is_iqr_outlier", "iqr_direction"]].assign(method="iqr"))
    if not flags:
        return pd.DataFrame()
    out = pd.concat(flags, ignore_index=True)
    out = out.rename(columns={service_col: "service_code", billing_col: "submitted_charge"})
    return _stable_sort(out, ["service_code", "provider_ccn"])


def zscore_method(df: pd.DataFrame, service_col: str, billing_col: str,
                  threshold: float = None) -> pd.DataFrame:
    """Per service, flag providers |z| > threshold on billing."""
    threshold = threshold or ZSCORE_MODERATE
    flags = []
    for code, sub in _per_service(df, service_col, billing_col, billing_col):
        mu = sub[billing_col].mean()
        sigma = sub[billing_col].std(ddof=1)
        if sigma == 0 or pd.isna(sigma):
            continue
        sub = sub.copy()
        sub["z"] = (sub[billing_col] - mu) / sigma
        sub["is_z_outlier"] = sub["z"].abs() > threshold
        sub["z_direction"] = np.where(sub["z"] > threshold, "high",
                                       np.where(sub["z"] < -threshold, "low", "none"))
        flags.append(sub[sub["is_z_outlier"]][[
            "provider_ccn", "provider_name", "state", "census_region", "urban_rural",
            service_col, billing_col, "z", "is_z_outlier", "z_direction"]].assign(method=f"zscore_{threshold}"))
    if not flags:
        return pd.DataFrame()
    out = pd.concat(flags, ignore_index=True)
    out = out.rename(columns={service_col: "service_code", billing_col: "submitted_charge"})
    return _stable_sort(out, ["service_code", "provider_ccn"])


def method_agreement_matrix(inp_high_charge_flags, inp_iqr_flags, inp_z_flags,
                            out_high_charge_flags, out_iqr_flags, out_z_flags) -> dict:
    """Return provider-level agreement counts for the three screening methods."""
    def _providers(*flag_dfs):
        sets = [set(df["provider_ccn"].unique()) if isinstance(df, pd.DataFrame) and not df.empty else set()
                for df in flag_dfs]
        return sets

    results = {}
    for ds, high_charge, iqr, z in [
        ("inpatient", inp_high_charge_flags, inp_iqr_flags, inp_z_flags),
        ("outpatient", out_high_charge_flags, out_iqr_flags, out_z_flags),
    ]:
        high_charge = _flagged_rows(high_charge)
        high_s, i_s, z_s = _providers(high_charge, iqr, z)
        three = len(high_s & i_s & z_s)
        any_one = len(high_s | i_s | z_s)
        results[ds] = {
            "high_charge_only": len(high_s - i_s - z_s),
            "iqr_only": len(i_s - high_s - z_s),
            "z_only": len(z_s - high_s - i_s),
            "high_charge_and_iqr": len((high_s & i_s) - z_s),
            "high_charge_and_z": len((high_s & z_s) - i_s),
            "iqr_and_z": len((i_s & z_s) - high_s),
            "all_three": three,
            "any_method": any_one,
            "high_charge_method": HIGH_CHARGE_METHOD,
        }
    return {"metadata": _analysis_metadata(), **results}


# ---------------------------------------------------------------------------
# 4c. Concentration analysis
# ---------------------------------------------------------------------------

def concentration_analysis(flags_df: pd.DataFrame, dataset_name: str,
                            top_n: int = 16) -> dict:
    """Rank services by provider-service high-charge flag counts."""
    flags = _flagged_rows(flags_df)
    if flags.empty:
        counts = pd.DataFrame(columns=[
            "service_code", "service_desc", "provider_service_flag_count",
            "share_of_total_provider_service_flags", "cumulative_share",
        ])
    else:
        pair_columns = [column for column in ("provider_ccn", "service_code") if column in flags.columns]
        flags = flags.drop_duplicates(pair_columns) if len(pair_columns) == 2 else flags
        counts = (
            flags.groupby(["service_code", "service_desc"], dropna=False, as_index=False)
            .size()
            .rename(columns={"size": "provider_service_flag_count"})
        )
        counts = _stable_sort(
            counts,
            ["provider_service_flag_count", "service_code", "service_desc"],
        )
        counts = counts.sort_values(
            ["provider_service_flag_count", "service_code", "service_desc"],
            ascending=[False, True, True],
            kind="mergesort",
            na_position="last",
        ).reset_index(drop=True)
        total = int(counts["provider_service_flag_count"].sum())
        counts["share_of_total_provider_service_flags"] = (
            counts["provider_service_flag_count"] / total if total else 0.0
        )
        counts["cumulative_share"] = counts["share_of_total_provider_service_flags"].cumsum()

    total = int(counts["provider_service_flag_count"].sum()) if not counts.empty else 0
    top16_share = (
        float(round(counts["share_of_total_provider_service_flags"].head(16).sum(), 4))
        if total else 0.0
    )
    crossing_target = (
        int((counts["cumulative_share"] < CONCENTRATION_TARGET_SHARE).sum()) + 1
        if total else 0
    )
    crossing_target = min(crossing_target, len(counts)) if counts.size else 0
    cumulative_top16 = (
        float(round(counts["cumulative_share"].iloc[15] if len(counts) >= 16
                   else counts["cumulative_share"].iloc[-1], 4))
        if total else 0.0
    )

    top_services = []
    if not counts.empty:
        preview = counts.head(top_n).copy()
        preview["share_of_total"] = preview["share_of_total_provider_service_flags"]
        preview["outlier_count"] = preview["provider_service_flag_count"]
        top_services = preview.to_dict(orient="records")

    return {
        "dataset": dataset_name,
        "metadata": _analysis_metadata(),
        "flag_unit": "provider-service pair",
        "flag_label": "high-charge flag",
        "total_flagged_provider_service_pairs": total,
        "total_outlier_events": total,
        "n_services_with_high_charge_flags": int(len(counts)),
        "n_services_with_outliers": int(len(counts)),
        "top16_provider_service_flag_share": top16_share,
        "our_top16_share": top16_share,
        "smallest_n_accounting_for_target_share": int(crossing_target),
        "smallest_n_accounting_for_40_pct": int(crossing_target),
        "target_share": CONCENTRATION_TARGET_SHARE,
        "cumulative_share_at_top16": cumulative_top16,
        "payment_dollar_benchmark_comparable": False,
        "oig_payment_dollar_benchmark": {
            "top16_share": 0.41,
            "unit": "payment dollars",
            "comparable": False,
            "comparison_status": "not_comparable",
            "used_for_calculation": False,
            "reason": (
                "The historical benchmark is a payment-dollar share, while "
                "this output counts provider-service high-charge flags."
            ),
        },
        "top_services": top_services,
    }


# ---------------------------------------------------------------------------
# 4d. Cross-dataset high-charge overlap
# ---------------------------------------------------------------------------

def _eligible_opportunity_info(frame: pd.DataFrame, service_col: str,
                               billing_col: str, desc_col: str) -> dict:
    """Build eligible-service and provider-service opportunity denominators."""
    empty_pairs = pd.DataFrame(columns=["provider_ccn", "service_code"])
    if not isinstance(frame, pd.DataFrame) or frame.empty:
        return {
            "available": False,
            "eligible_codes": set(),
            "eligible_provider_set": set(),
            "provider_eligible_counts": pd.Series(dtype=int),
            "eligible_pairs": empty_pairs,
            "by_service": [],
        }
    required = ["provider_ccn", service_col, billing_col]
    if any(column not in frame.columns for column in required):
        return {
            "available": False,
            "eligible_codes": set(),
            "eligible_provider_set": set(),
            "provider_eligible_counts": pd.Series(dtype=int),
            "eligible_pairs": empty_pairs,
            "by_service": [],
        }

    working = frame.dropna(subset=["provider_ccn", service_col, billing_col]).copy()
    working = working.drop_duplicates(["provider_ccn", service_col])
    service_counts = (
        working.groupby(service_col, sort=True)["provider_ccn"]
        .nunique()
        .rename("eligible_provider_count")
        .reset_index()
    )
    eligible_codes = set(
        service_counts.loc[
            service_counts["eligible_provider_count"] >= MIN_PROVIDERS_PER_SERVICE,
            service_col,
        ].tolist()
    )
    eligible = working[working[service_col].isin(eligible_codes)].copy()
    eligible = eligible.rename(columns={service_col: "service_code"})
    provider_counts = (
        eligible.groupby("provider_ccn")["service_code"].nunique()
        if not eligible.empty else pd.Series(dtype=int)
    )

    by_service = []
    if not eligible.empty:
        desc_by_service = (
            working.groupby(service_col, sort=True)[desc_col].first()
            if desc_col in working.columns else pd.Series(dtype=object)
        )
        by_service_df = (
            eligible.groupby("service_code", sort=True)
            .agg(
                eligible_provider_count=("provider_ccn", "nunique"),
                eligible_service_opportunities=("provider_ccn", "size"),
            )
            .reset_index()
        )
        for _, row in by_service_df.iterrows():
            by_service.append(OrderedDict(
                service_code=row["service_code"],
                service_desc=desc_by_service.get(row["service_code"]),
                eligible_provider_count=int(row["eligible_provider_count"]),
                eligible_service_opportunities=int(row["eligible_service_opportunities"]),
            ))

    return {
        "available": True,
        "eligible_codes": eligible_codes,
        "eligible_provider_set": set(provider_counts.index),
        "provider_eligible_counts": provider_counts,
        "eligible_pairs": eligible[["provider_ccn", "service_code"]].drop_duplicates(),
        "by_service": by_service,
    }


def _flag_opportunity_info(flags_df: pd.DataFrame) -> dict:
    """Use a flagged table as a clearly labeled denominator fallback."""
    empty = {
        "available": False,
        "eligible_codes": set(),
        "eligible_provider_set": set(),
        "provider_eligible_counts": pd.Series(dtype=int),
        "eligible_pairs": pd.DataFrame(columns=["provider_ccn", "service_code"]),
        "by_service": [],
    }
    if not isinstance(flags_df, pd.DataFrame) or flags_df.empty:
        return empty
    required = {"provider_ccn", "service_code"}
    if not required.issubset(flags_df.columns):
        return empty
    pairs = flags_df.dropna(subset=["provider_ccn", "service_code"])[
        ["provider_ccn", "service_code"]
    ].drop_duplicates()
    if pairs.empty:
        return empty
    desc_by_service = (
        flags_df.groupby("service_code", sort=True)["service_desc"].first()
        if "service_desc" in flags_df.columns else pd.Series(dtype=object)
    )
    by_service_df = (
        pairs.groupby("service_code", sort=True)
        .agg(
            eligible_provider_count=("provider_ccn", "nunique"),
            eligible_service_opportunities=("provider_ccn", "size"),
        )
        .reset_index()
    )
    by_service = [OrderedDict(
        service_code=row["service_code"],
        service_desc=desc_by_service.get(row["service_code"]),
        eligible_provider_count=int(row["eligible_provider_count"]),
        eligible_service_opportunities=int(row["eligible_service_opportunities"]),
    ) for _, row in by_service_df.iterrows()]
    provider_counts = pairs.groupby("provider_ccn")["service_code"].nunique()
    return {
        "available": False,
        "eligible_codes": set(pairs["service_code"]),
        "eligible_provider_set": set(provider_counts.index),
        "provider_eligible_counts": provider_counts,
        "eligible_pairs": pairs,
        "by_service": by_service,
    }


def _universe_stats(frame: pd.DataFrame, service_col: str, payment_col: str,
                    eligible_codes: set, available: bool) -> dict:
    """Summarize full and cost-observed provider-service universes."""
    if not available or not isinstance(frame, pd.DataFrame):
        return {"available": False}
    if frame.empty or "provider_ccn" not in frame.columns or service_col not in frame.columns:
        return {"available": True, "rows": 0, "providers": 0, "services": 0,
                "provider_service_pairs": 0, "eligible_service_opportunities": 0,
                "eligible_providers": 0, "eligible_services": 0}

    pairs = frame.dropna(subset=["provider_ccn", service_col]).drop_duplicates(
        ["provider_ccn", service_col]
    )
    eligible = pairs[pairs[service_col].isin(eligible_codes)]
    cost_observed = (
        frame[frame[payment_col].notna()]
        if payment_col in frame.columns else None
    )
    cost_pairs = (
        cost_observed.dropna(subset=["provider_ccn", service_col]).drop_duplicates(
            ["provider_ccn", service_col]
        ) if cost_observed is not None else None
    )
    return OrderedDict(
        available=True,
        rows=int(len(frame)),
        providers=int(frame["provider_ccn"].nunique()),
        services=int(frame[service_col].nunique()),
        provider_service_pairs=int(len(pairs)),
        eligible_services=int(
            pairs.loc[pairs[service_col].isin(eligible_codes), service_col].nunique()
        ),
        eligible_service_opportunities=int(len(eligible)),
        eligible_providers=int(eligible["provider_ccn"].nunique()),
        cost_observed_rows=(int(len(cost_observed)) if cost_observed is not None else None),
        cost_observed_provider_service_pairs=(
            int(len(cost_pairs)) if cost_pairs is not None else None
        ),
    )


def _breadth_bin(count: int) -> str:
    if count <= 1:
        return "1"
    if count <= 4:
        return "2-4"
    if count <= 9:
        return "5-9"
    return "10+"


def _breadth_bins(breadth: pd.Series) -> list:
    counts = breadth.value_counts() if isinstance(breadth, pd.Series) else pd.Series(dtype=int)
    total = int(len(breadth)) if isinstance(breadth, pd.Series) else 0
    rows = []
    for label in BREADTH_BIN_LABELS:
        provider_count = int(sum(
            int(frequency) for breadth_value, frequency in counts.items()
            if _breadth_bin(int(breadth_value)) == label
        ))
        rows.append(OrderedDict(
            breadth_bin=label,
            provider_count=provider_count,
            share_of_flagged_providers=_rate(provider_count, total),
        ))
    return rows


def _method_flags(flags_df: pd.DataFrame, method: str) -> pd.DataFrame:
    if not isinstance(flags_df, pd.DataFrame) or flags_df.empty:
        return flags_df.copy() if isinstance(flags_df, pd.DataFrame) else pd.DataFrame()
    if "outlier_method" in flags_df.columns:
        return flags_df[flags_df["outlier_method"] == method].copy()
    return flags_df.copy()


def cross_dataset_outliers(inp_flags: pd.DataFrame, out_flags: pd.DataFrame,
                          method: str = HIGH_CHARGE_METHOD,
                          inp_frame: pd.DataFrame = None,
                          out_frame: pd.DataFrame = None,
                          inp_full_frame: pd.DataFrame = None,
                          out_full_frame: pd.DataFrame = None,
                          metadata: dict = None,
                          ) -> dict:
    """Measure provider overlap using high-charge flags and full denominators."""
    inp_method_flags = _method_flags(inp_flags, method)
    out_method_flags = _method_flags(out_flags, method)
    inp_sub = _flagged_rows(inp_method_flags)
    out_sub = _flagged_rows(out_method_flags)

    inp_info = _eligible_opportunity_info(
        inp_frame, "drg_code", "avg_submitted_charge", "drg_desc"
    )
    out_info = _eligible_opportunity_info(
        out_frame, "apc_code", "avg_submitted_charge", "apc_desc"
    )
    # When raw frames are unavailable, retain a transparent flag-only fallback.
    inp_frame_available = inp_info["available"]
    out_frame_available = out_info["available"]
    if not inp_info["available"]:
        inp_info = _eligible_opportunity_info(
            inp_method_flags, "service_code", "submitted_charge", "service_desc"
        )
        if not inp_info["available"]:
            inp_info = _flag_opportunity_info(inp_sub)
    if not out_info["available"]:
        out_info = _eligible_opportunity_info(
            out_method_flags, "service_code", "submitted_charge", "service_desc"
        )
        if not out_info["available"]:
            out_info = _flag_opportunity_info(out_sub)
    inp_denominator_basis = (
        "eligible provider-service rows in the input frame"
        if inp_frame_available
        else (
            "complete cohort flag table"
            if inp_info["available"] else "flag table fallback; input frame unavailable"
        )
    )
    out_denominator_basis = (
        "eligible provider-service rows in the input frame"
        if out_frame_available
        else (
            "complete cohort flag table"
            if out_info["available"] else "flag table fallback; input frame unavailable"
        )
    )

    inp_p = set(inp_sub["provider_ccn"].dropna().unique()) if not inp_sub.empty else set()
    out_p = set(out_sub["provider_ccn"].dropna().unique()) if not out_sub.empty else set()
    common = inp_p & out_p

    inp_pairs = (inp_sub[["provider_ccn", "service_code"]].drop_duplicates()
                 if not inp_sub.empty else pd.DataFrame(columns=["provider_ccn", "service_code"]))
    out_pairs = (out_sub[["provider_ccn", "service_code"]].drop_duplicates()
                 if not out_sub.empty else pd.DataFrame(columns=["provider_ccn", "service_code"]))
    inp_breadth = (inp_pairs.groupby("provider_ccn")["service_code"].nunique()
                   if not inp_pairs.empty else pd.Series(dtype=int))
    out_breadth = (out_pairs.groupby("provider_ccn")["service_code"].nunique()
                   if not out_pairs.empty else pd.Series(dtype=int))

    inp_eligible_provider_count = len(inp_info["eligible_provider_set"])
    out_eligible_provider_count = len(out_info["eligible_provider_set"])
    inp_flagged_provider_rate = _rate(len(inp_p), inp_eligible_provider_count)
    out_flagged_provider_rate = _rate(len(out_p), out_eligible_provider_count)
    eligible_both = (
        inp_info["eligible_provider_set"] & out_info["eligible_provider_set"]
    )
    marginal_expected_overlap = (
        len(eligible_both) * inp_flagged_provider_rate * out_flagged_provider_rate
    )
    shared_inp_flagged = inp_p & eligible_both
    shared_out_flagged = out_p & eligible_both
    shared_inp_rate = _rate(len(shared_inp_flagged), len(eligible_both))
    shared_out_rate = _rate(len(shared_out_flagged), len(eligible_both))
    expected_overlap = len(eligible_both) * shared_inp_rate * shared_out_rate
    breadth_expected_overlap, breadth_detail = _breadth_matched_expected_overlap(
        eligible_both,
        inp_p,
        out_p,
        inp_info["provider_eligible_counts"],
        out_info["provider_eligible_counts"],
    )

    inp_universe = OrderedDict(
        full=_universe_stats(
            inp_full_frame, "drg_code", "avg_total_payment",
            inp_info["eligible_codes"], inp_full_frame is not None,
        ),
        cost_observed=_universe_stats(
            inp_frame, "drg_code", "avg_total_payment",
            inp_info["eligible_codes"], inp_frame is not None,
        ),
    )
    out_universe = OrderedDict(
        full=_universe_stats(
            out_full_frame, "apc_code", "avg_allowed_amount",
            out_info["eligible_codes"], out_full_frame is not None,
        ),
        cost_observed=_universe_stats(
            out_frame, "apc_code", "avg_allowed_amount",
            out_info["eligible_codes"], out_frame is not None,
        ),
    )

    # Build provider metadata keyed by provider_ccn.
    meta_frames = []
    for frame in (inp_sub, out_sub):
        if frame.empty or "provider_ccn" not in frame.columns:
            continue
        columns = [
            column for column in (
                "provider_name", "state", "census_region", "urban_rural",
            ) if column in frame.columns
        ]
        if columns:
            meta_frames.append(frame.drop_duplicates("provider_ccn").set_index("provider_ccn")[columns])
    meta = (
        pd.concat(meta_frames).groupby(level=0).first()
        if meta_frames else pd.DataFrame()
    )
    inp_provider_eligible = inp_info["provider_eligible_counts"]
    out_provider_eligible = out_info["provider_eligible_counts"]
    rows = []
    for ccn in sorted(common, key=str):
        inp_flag_count = int(inp_breadth.get(ccn, 0))
        out_flag_count = int(out_breadth.get(ccn, 0))
        inp_eligible_count = int(inp_provider_eligible.get(ccn, inp_flag_count))
        out_eligible_count = int(out_provider_eligible.get(ccn, out_flag_count))
        row = OrderedDict(
            provider_ccn=ccn,
            provider_name=meta.loc[ccn, "provider_name"]
                if ccn in meta.index and "provider_name" in meta.columns else None,
            state=meta.loc[ccn, "state"]
                if ccn in meta.index and "state" in meta.columns else None,
            census_region=meta.loc[ccn, "census_region"]
                if ccn in meta.index and "census_region" in meta.columns else None,
            urban_rural=meta.loc[ccn, "urban_rural"]
                if ccn in meta.index and "urban_rural" in meta.columns else None,
            inpatient_high_charge_services_count=inp_flag_count,
            outpatient_high_charge_services_count=out_flag_count,
            inpatient_eligible_services_count=inp_eligible_count,
            outpatient_eligible_services_count=out_eligible_count,
            inpatient_provider_flag_rate=_rate(inp_flag_count, inp_eligible_count),
            outpatient_provider_flag_rate=_rate(out_flag_count, out_eligible_count),
            inpatient_breadth_bin=_breadth_bin(inp_flag_count),
            outpatient_breadth_bin=_breadth_bin(out_flag_count),
            total_high_charge_services_count=inp_flag_count + out_flag_count,
            # Existing generic count fields remain useful to current consumers.
            inp_outlier_services_count=inp_flag_count,
            out_outlier_services_count=out_flag_count,
        )
        rows.append(row)
    rows = sorted(
        rows,
        key=lambda row: (-row["total_high_charge_services_count"], str(row["provider_ccn"])),
    )

    inp_summary = OrderedDict(
        flagged_provider_service_pairs=int(len(inp_pairs)),
        eligible_service_opportunities=int(len(inp_info["eligible_pairs"])),
        provider_service_flag_rate=_rate(len(inp_pairs), len(inp_info["eligible_pairs"])),
        provider_service_flag_rate_pct=float(round(
            100 * _rate(len(inp_pairs), len(inp_info["eligible_pairs"])), 2
        )),
        denominator_basis=inp_denominator_basis,
        eligible_service_count=int(len(inp_info["eligible_codes"])),
        flagged_provider_count=int(len(inp_p)),
        eligible_provider_count=int(inp_eligible_provider_count),
        provider_flag_rate=inp_flagged_provider_rate,
        provider_flag_rate_pct=float(round(100 * inp_flagged_provider_rate, 2)),
    )
    out_summary = OrderedDict(
        flagged_provider_service_pairs=int(len(out_pairs)),
        eligible_service_opportunities=int(len(out_info["eligible_pairs"])),
        provider_service_flag_rate=_rate(len(out_pairs), len(out_info["eligible_pairs"])),
        provider_service_flag_rate_pct=float(round(
            100 * _rate(len(out_pairs), len(out_info["eligible_pairs"])), 2
        )),
        denominator_basis=out_denominator_basis,
        eligible_service_count=int(len(out_info["eligible_codes"])),
        flagged_provider_count=int(len(out_p)),
        eligible_provider_count=int(out_eligible_provider_count),
        provider_flag_rate=out_flagged_provider_rate,
        provider_flag_rate_pct=float(round(100 * out_flagged_provider_rate, 2)),
    )
    expected_detail = OrderedDict(
        benchmark="conditional independence within the shared eligible-provider universe",
        eligible_provider_overlap=int(len(eligible_both)),
        inpatient_provider_flag_rate=shared_inp_rate,
        outpatient_provider_flag_rate=shared_out_rate,
        expected_flagged_provider_overlap=float(round(expected_overlap, 4)),
        marginal_expected_flagged_provider_overlap=float(round(marginal_expected_overlap, 4)),
        breadth_matched_expected_flagged_provider_overlap=float(
            round(breadth_expected_overlap, 4)
        ),
        observed_flagged_provider_overlap=int(len(common)),
        observed_to_expected_ratio=(
            float(round(len(common) / expected_overlap, 4)) if expected_overlap else None
        ),
        observed_to_breadth_matched_expected_ratio=(
            float(round(len(common) / breadth_expected_overlap, 4))
            if breadth_expected_overlap else None
        ),
        assumption=(
            "Within the shared eligible-provider universe, the two flag indicators "
            "are treated as independent; service-opportunity breadth is not matched."
        ),
        breadth_matched_benchmark=breadth_detail,
    )
    output_metadata = metadata or _analysis_metadata(
        _period_from_frames(inp_frame, out_frame, inp_full_frame, out_full_frame)
    )
    return {
        "metadata": output_metadata,
        "summary": {
            "method": method,
            "high_charge_method": HIGH_CHARGE_METHOD,
            "inpatient_flagged_providers": int(len(inp_p)),
            "outpatient_flagged_providers": int(len(out_p)),
            "inpatient_eligible_providers": int(inp_eligible_provider_count),
            "outpatient_eligible_providers": int(out_eligible_provider_count),
            "inpatient_provider_flag_rate": inp_flagged_provider_rate,
            "outpatient_provider_flag_rate": out_flagged_provider_rate,
            "shared_eligible_provider_count": int(len(eligible_both)),
            "shared_inpatient_provider_flag_rate": shared_inp_rate,
            "shared_outpatient_provider_flag_rate": shared_out_rate,
            "inpatient_flagged_provider_service_pairs": int(len(inp_pairs)),
            "outpatient_flagged_provider_service_pairs": int(len(out_pairs)),
            "inpatient_eligible_service_opportunities": int(len(inp_info["eligible_pairs"])),
            "outpatient_eligible_service_opportunities": int(len(out_info["eligible_pairs"])),
            "inpatient_outlier_providers": int(len(inp_p)),
            "outpatient_outlier_providers": int(len(out_p)),
            "outliers_in_both": int(len(common)),
            "pct_of_inpatient": float(round(100 * len(common) / max(1, len(inp_p)), 2)),
            "pct_of_outpatient": float(round(100 * len(common) / max(1, len(out_p)), 2)),
            "expected_overlap_under_independence": float(round(expected_overlap, 4)),
            "marginal_expected_overlap_under_independence": float(
                round(marginal_expected_overlap, 4)
            ),
            "expected_overlap_under_breadth_matched_independence": float(
                round(breadth_expected_overlap, 4)
            ),
            "expected_overlap_under_independence_detail": expected_detail,
        },
        "eligible_service_opportunity_denominators": {
            "inpatient": inp_info["by_service"],
            "outpatient": out_info["by_service"],
        },
        "dataset_denominators": {
            "inpatient": inp_summary,
            "outpatient": out_summary,
        },
        "breadth_bins": {
            "inpatient": _breadth_bins(inp_breadth),
            "outpatient": _breadth_bins(out_breadth),
        },
        "universe_metadata": {
            "inpatient": inp_universe,
            "outpatient": out_universe,
        },
        "top_cross_dataset_outliers": rows[:50],
    }


# ---------------------------------------------------------------------------
# 4e. Outlier profiling
# ---------------------------------------------------------------------------

def profile_outliers(inp_flags: pd.DataFrame, out_flags: pd.DataFrame) -> dict:
    """Geographic + breadth profile of outlier providers per dataset.
    For each dataset we report, by census_region and urban_rural, the count of
    distinct outlier providers + their mean number of services they're outliers in.
    """
    def _prof(flags_df, ds_name):
        flags_df = _flagged_rows(flags_df)
        if flags_df.empty:
            return {"dataset": ds_name, "regions": [], "urban_rural": [], "breadth": {}}
        # Distinct providers per (geography)
        per_prov_breadth = flags_df.groupby(["provider_ccn", "census_region", "urban_rural"])["service_code"].nunique().reset_index(name="service_count")
        region_agg = (per_prov_breadth.groupby("census_region")
                      .agg(n_outlier_providers=("provider_ccn", "nunique"),
                           median_services_outlier_in=("service_count", "median"))
                      .reset_index())
        urban_agg = (per_prov_breadth.groupby("urban_rural")
                     .agg(n_outlier_providers=("provider_ccn", "nunique"),
                          median_services_outlier_in=("service_count", "median"))
                     .reset_index())
        breadth = per_prov_breadth["service_count"]
        return {
            "dataset": ds_name,
            "regions": region_agg.to_dict(orient="records"),
            "urban_rural": urban_agg.to_dict(orient="records"),
            "breadth": {
                "median_services_outlier_in": float(round(breadth.median(), 2)),
                "p90_services_outlier_in": float(round(breadth.quantile(0.9), 2)),
                "max_services_outlier_in": int(breadth.max()),
                "n_providers_outlier_in_one_service_only": int((breadth == 1).sum()),
                "n_providers_outlier_in_5plus_services": int((breadth >= 5).sum()),
            },
        }
    return {
        "metadata": _analysis_metadata(),
        "inpatient": _prof(inp_flags, "inpatient"),
        "outpatient": _prof(out_flags, "outpatient"),
    }


def _audit_frame(frame: pd.DataFrame, dataset: str) -> pd.DataFrame:
    """Add dataset identity and deterministic ordering to an audit table."""
    if not isinstance(frame, pd.DataFrame):
        frame = pd.DataFrame()
    frame = frame.copy()
    if "dataset" in frame.columns:
        frame["dataset"] = dataset
    else:
        frame.insert(0, "dataset", dataset)
    return _stable_sort(frame, ["dataset", "service_code", "provider_ccn"])


def _write_audit_table(frame: pd.DataFrame, stem: str) -> str:
    """Write a complete audit table to processed data, with CSV fallback."""
    parquet_path = DATA_PROC / f"{stem}.parquet"
    try:
        frame.to_parquet(parquet_path, index=False)
        return str(parquet_path)
    except Exception:
        csv_path = DATA_PROC / f"{stem}.csv"
        frame.to_csv(csv_path, index=False)
        return str(csv_path)


def _preview_metadata(metadata: dict, full_counts: dict, preview_counts: dict,
                      limit: int, sort_order: list) -> OrderedDict:
    """Attach explicit size and ordering information to a site preview."""
    out = OrderedDict(metadata)
    out.update(
        preview=True,
        truncated=any(
            preview_counts.get(dataset, 0) < full_counts.get(dataset, 0)
            for dataset in full_counts
        ),
        preview_limit_per_dataset=int(limit),
        full_row_counts={key: int(value) for key, value in full_counts.items()},
        preview_row_counts={key: int(value) for key, value in preview_counts.items()},
        sort_order=sort_order,
    )
    return out


# ---------------------------------------------------------------------------
# Orchestration
# ---------------------------------------------------------------------------

def run() -> dict:
    print("[Phase 4] Loading feature parquets ...")
    inp = pd.read_parquet(DATA_PROC / "inpatient_features.parquet")
    out = pd.read_parquet(DATA_PROC / "outpatient_features.parquet")
    try:
        out_full = pd.read_parquet(DATA_PROC / "outpatient_full.parquet")
    except FileNotFoundError:
        out_full = None
    print(f"      inpatient: {len(inp):,} rows | outpatient: {len(out):,} rows")

    output_metadata = _analysis_metadata()

    # 4a — project-defined high-charge cohort (top quartile per service)
    print("[Phase 4] 4a — Project-defined high-charge cohort (top quartile per service) ...")
    inp_high_charge = high_charge_cohort_method(
        inp, "drg_code", "drg_desc", "avg_submitted_charge", "avg_total_payment",
        "total_discharges", "inpatient"
    )
    out_high_charge = high_charge_cohort_method(
        out, "apc_code", "apc_desc", "avg_submitted_charge", "avg_allowed_amount",
        "apc_services", "outpatient"
    )
    inp_high_charge_flags = _flagged_rows(inp_high_charge["flags"])
    out_high_charge_flags = _flagged_rows(out_high_charge["flags"])
    n_inp_flags = len(inp_high_charge_flags)
    n_out_flags = len(out_high_charge_flags)
    print(f"      inpatient cohort: {n_inp_flags:,} provider-service high-charge flags")
    print(f"      outpatient cohort: {n_out_flags:,} provider-service high-charge flags")

    # Cohort comparison headline
    inp_cohort = inp_high_charge["cohort_comparison"]
    out_cohort = out_high_charge["cohort_comparison"]
    median_inp_premium = inp_cohort["charge_premium_pct"].median() if not inp_cohort.empty else None
    median_out_premium = out_cohort["charge_premium_pct"].median() if not out_cohort.empty else None
    median_inp_proxy_prem = (
        inp_cohort["charge_to_payment_proxy_premium_pct"].median()
        if not inp_cohort.empty else None
    )
    median_out_proxy_prem = (
        out_cohort["charge_to_payment_proxy_premium_pct"].median()
        if not out_cohort.empty else None
    )
    print(f"      inpatient high-charge charge premium (median): {median_inp_premium}%")
    print(f"      outpatient high-charge charge premium (median): {median_out_premium}%")
    print(f"      inpatient charge-to-payment proxy premium (median): {median_inp_proxy_prem}%")
    print(f"      outpatient charge-to-payment proxy premium (median): {median_out_proxy_prem}%")

    # 4b — IQR + z-score methods
    print("[Phase 4] 4b — IQR + z-score (per service) ...")
    inp_iqr_flags = iqr_method(inp, "drg_code", "avg_submitted_charge")
    out_iqr_flags = iqr_method(out, "apc_code", "avg_submitted_charge")
    inp_z2_flags = zscore_method(inp, "drg_code", "avg_submitted_charge", ZSCORE_MODERATE)
    out_z2_flags = zscore_method(out, "apc_code", "avg_submitted_charge", ZSCORE_MODERATE)
    print(f"      inpatient: IQR {len(inp_iqr_flags):,} | z>2 {len(inp_z2_flags):,}")
    print(f"      outpatient: IQR {len(out_iqr_flags):,} | z>2 {len(out_z2_flags):,}")

    agreement = method_agreement_matrix(
        inp_high_charge_flags, inp_iqr_flags, inp_z2_flags,
        out_high_charge_flags, out_iqr_flags, out_z2_flags,
    )
    agreement["metadata"] = output_metadata
    print(f"      3-way agreement (all three methods): "
          f"inpatient={agreement['inpatient']['all_three']} | "
          f"outpatient={agreement['outpatient']['all_three']}")
    print(f"      total inpatient providers flagged by any method: {agreement['inpatient']['any_method']}")

    # 4c — Concentration analysis
    print("[Phase 4] 4c — Concentration analysis (provider-service high-charge flags) ...")
    inp_conc = concentration_analysis(inp_high_charge_flags, "inpatient")
    out_conc = concentration_analysis(out_high_charge_flags, "outpatient")
    inp_conc["metadata"] = output_metadata
    out_conc["metadata"] = output_metadata
    print(f"      inpatient: top 16 provider-service flag share = "
          f"{inp_conc['top16_provider_service_flag_share'] * 100:.1f}% | "
          f"smallest N for 40% = {inp_conc['smallest_n_accounting_for_target_share']}")
    print(f"      outpatient: top 16 provider-service flag share = "
          f"{out_conc['top16_provider_service_flag_share'] * 100:.1f}% | "
          f"smallest N for 40% = {out_conc['smallest_n_accounting_for_target_share']}")
    print("      historical payment-dollar benchmark is not comparable to these flag counts")

    # 4d — Cross-dataset high-charge overlap
    print("[Phase 4] 4d — Cross-dataset high-charge flag overlap ...")
    cross = cross_dataset_outliers(
        inp_high_charge_flags,
        out_high_charge_flags,
        method=HIGH_CHARGE_METHOD,
        inp_frame=inp,
        out_frame=out,
        inp_full_frame=inp,
        out_full_frame=out_full,
        metadata=output_metadata,
    )
    print(f"      inpatient flagged providers: {cross['summary']['inpatient_flagged_providers']:,}")
    print(f"      outpatient flagged providers: {cross['summary']['outpatient_flagged_providers']:,}")
    print(f"      providers flagged in BOTH: {cross['summary']['outliers_in_both']:,} "
          f"({cross['summary']['pct_of_inpatient']}% of inpatient / "
          f"{cross['summary']['pct_of_outpatient']}% of outpatient)")
    print(f"      expected overlap under independence: "
          f"{cross['summary']['expected_overlap_under_independence']:.2f}")

    # 4e — Profiling
    print("[Phase 4] 4e — High-charge provider profiling (geography, breadth) ...")
    profile = profile_outliers(inp_high_charge_flags, out_high_charge_flags)
    print(f"      inpatient high-charge breadth: median "
          f"{profile['inpatient']['breadth']['median_services_outlier_in']} services "
          f"(max {profile['inpatient']['breadth']['max_services_outlier_in']})")

    # Serialize
    print("[Phase 4] Serializing to site/data/ ...")
    cohort_sort = ["service_code", "service_desc"]
    cohort_preview_inp = _stable_sort(inp_cohort, cohort_sort)
    cohort_preview_out = _stable_sort(out_cohort, cohort_sort)
    cohort_full_counts = {
        "inpatient": len(cohort_preview_inp),
        "outpatient": len(cohort_preview_out),
    }
    cohort_preview_counts = {
        "inpatient": min(COHORT_PREVIEW_LIMIT, len(cohort_preview_inp)),
        "outpatient": min(COHORT_PREVIEW_LIMIT, len(cohort_preview_out)),
    }
    save_json({
        "metadata": _preview_metadata(
            output_metadata, cohort_full_counts, cohort_preview_counts,
            COHORT_PREVIEW_LIMIT, cohort_sort,
        ),
        "inpatient": _records(cohort_preview_inp, COHORT_PREVIEW_LIMIT),
        "outpatient": _records(cohort_preview_out, COHORT_PREVIEW_LIMIT),
        "summary": {
            "inpatient_median_charge_premium_pct": _rounded(median_inp_premium, 2),
            "outpatient_median_charge_premium_pct": _rounded(median_out_premium, 2),
            "inpatient_median_charge_to_payment_proxy_premium_pct": _rounded(
                median_inp_proxy_prem, 2
            ),
            "outpatient_median_charge_to_payment_proxy_premium_pct": _rounded(
                median_out_proxy_prem, 2
            ),
        },
    }, "outlier_cohort_comparison.json")

    save_json({**agreement, "metadata": output_metadata}, "outlier_method_agreement.json")

    save_json({
        "metadata": output_metadata,
        "inpatient": inp_conc,
        "outpatient": out_conc,
    }, "outlier_concentration.json")

    save_json(cross, "cross_dataset_outliers.json")

    profile["metadata"] = output_metadata
    save_json(profile, "outlier_profile.json")

    flag_columns = [
        "provider_ccn", "provider_name", "state", "census_region", "urban_rural",
        "service_code", "service_desc", "is_high_charge_flag", "outlier_method",
        "eligible_service", "eligible_provider_count", "high_charge_threshold",
        "submitted_charge", "payment", "volume", "charge_to_payment_proxy",
    ]
    inp_site_flags = _stable_sort(inp_high_charge_flags, ["service_code", "provider_ccn"])
    out_site_flags = _stable_sort(out_high_charge_flags, ["service_code", "provider_ccn"])
    inp_site_flags = inp_site_flags[[column for column in flag_columns if column in inp_site_flags.columns]]
    out_site_flags = out_site_flags[[column for column in flag_columns if column in out_site_flags.columns]]
    flag_full_counts = {
        "inpatient": len(inp_site_flags),
        "outpatient": len(out_site_flags),
    }
    flag_preview_counts = {
        "inpatient": min(FLAG_PREVIEW_LIMIT, len(inp_site_flags)),
        "outpatient": min(FLAG_PREVIEW_LIMIT, len(out_site_flags)),
    }
    save_json({
        "metadata": _preview_metadata(
            output_metadata, flag_full_counts, flag_preview_counts,
            FLAG_PREVIEW_LIMIT, ["service_code", "provider_ccn"],
        ),
        "inpatient": _records(inp_site_flags, FLAG_PREVIEW_LIMIT),
        "outpatient": _records(out_site_flags, FLAG_PREVIEW_LIMIT),
    }, "outlier_flags.json")

    # Complete tables stay outside the website preview for auditability.
    high_charge_audit = pd.concat([
        _audit_frame(inp_high_charge["flags"], "inpatient"),
        _audit_frame(out_high_charge["flags"], "outpatient"),
    ], ignore_index=True, sort=False)
    cohort_audit = pd.concat([
        _audit_frame(inp_cohort, "inpatient"),
        _audit_frame(out_cohort, "outpatient"),
    ], ignore_index=True, sort=False)
    iqr_audit = pd.concat([
        _audit_frame(inp_iqr_flags, "inpatient"),
        _audit_frame(out_iqr_flags, "outpatient"),
    ], ignore_index=True, sort=False)
    zscore_audit = pd.concat([
        _audit_frame(inp_z2_flags, "inpatient"),
        _audit_frame(out_z2_flags, "outpatient"),
    ], ignore_index=True, sort=False)
    audit_artifacts = OrderedDict(
        outlier_flags_audit=_write_audit_table(high_charge_audit, "outlier_flags_audit"),
        outlier_cohort_comparison_audit=_write_audit_table(
            cohort_audit, "outlier_cohort_comparison_audit"
        ),
        outlier_iqr_flags_audit=_write_audit_table(iqr_audit, "outlier_iqr_flags_audit"),
        outlier_zscore_flags_audit=_write_audit_table(
            zscore_audit, "outlier_zscore_flags_audit"
        ),
    )

    print("[Phase 4] DONE.")
    return {
        "cohort_comparison": {"inpatient": inp_cohort, "outpatient": out_cohort},
        "agreement": agreement,
        "concentration": {"inpatient": inp_conc, "outpatient": out_conc},
        "cross_dataset": cross,
        "profile": profile,
        "metadata": output_metadata,
        "audit_artifacts": audit_artifacts,
    }


if __name__ == "__main__":  # pragma: no cover
    run()