stack/dev/scripts/build_skin_subs_anomaly.py

"""Build anomaly scoring and provider risk profiles for skin substitutes.

Addresses #242 (predatory pattern detection) and #243 (enforcement
action correlation).

Usage:
    uv run python dev/scripts/build_skin_subs_anomaly.py
"""

from __future__ import annotations

from pathlib import Path

import duckdb

ROOT = Path(__file__).resolve().parents[2]
DUCKDB_PATH = ROOT / "data" / "aco.duckdb"

# Known enforcement cases — providers we can match against
KNOWN_FRAUD_PATTERNS = {
    "jenson": {
        "specialty": "podiatry",
        "state": "TX",
        "setting": "office",
        "description": "USA v. Jenson — $90M, S.D. Tex., podiatry clinic",
    },
    "gehrke": {
        "specialty": "nurse_practitioner",
        "state": "AZ",
        "setting": "snf",
        "description": "USA v. Gehrke/King — $1.2B, D. Ariz., mobile wound care",
    },
    "vohra": {
        "specialty": "other",
        "state": "FL",
        "setting": "snf",
        "description": "Vohra Wound Physicians — $45M, S.D. Fla., EMR upcoding",
    },
}

# Setting risk weights (higher = more fraud-prone based on OIG findings)
SETTING_RISK = {
    "office": 1.5,  # Lowest oversight, wound care mills
    "snf": 1.3,  # Kickback-prone, captive patients
    "asc": 0.8,  # Moderate oversight
    "hopd": 0.5,  # Institutional controls, auditable
}


def build_anomaly_scores(con: duckdb.DuckDBPyConnection) -> None:
    """#242 — Composite anomaly scoring per provider."""
    print("\n=== Anomaly Scoring (#242) ===")

    con.execute("DROP TABLE IF EXISTS skin_subs.anomaly_scores")
    con.execute(f"""
        CREATE TABLE skin_subs.anomaly_scores AS
        WITH provider_stats AS (
            SELECT
                p.rendering_npi,
                p.provider_specialty,
                p.state,
                p.total_paid,
                p.unique_patients,
                p.paid_per_patient,
                p.total_claim_lines,
                p.unique_products,
                -- Dominant setting for this provider
                (SELECT place_of_service_description
                 FROM skin_subs.claims_synthetic c
                 WHERE c.rendering_npi = p.rendering_npi
                 GROUP BY place_of_service_description
                 ORDER BY sum(paid_amount) DESC LIMIT 1
                ) AS primary_setting,
                -- Avg units per product claim
                (SELECT round(avg(units), 1)
                 FROM skin_subs.claims_synthetic c
                 WHERE c.rendering_npi = p.rendering_npi
                   AND c.claim_type = 'product'
                ) AS avg_units
            FROM skin_subs.providers p
        ),
        peer_stats AS (
            SELECT
                provider_specialty,
                avg(total_paid) AS peer_avg_paid,
                stddev(total_paid) AS peer_std_paid,
                avg(paid_per_patient) AS peer_avg_per_patient,
                stddev(paid_per_patient) AS peer_std_per_patient,
                avg(unique_patients) AS peer_avg_patients,
                stddev(unique_patients) AS peer_std_patients,
                avg(avg_units) AS peer_avg_units,
                stddev(avg_units) AS peer_std_units
            FROM provider_stats
            GROUP BY provider_specialty
        )
        SELECT
            ps.rendering_npi,
            ps.provider_specialty,
            ps.state,
            ps.primary_setting,
            ps.total_paid,
            ps.unique_patients,
            ps.paid_per_patient,
            ps.total_claim_lines,
            ps.unique_products,
            ps.avg_units,
            -- Z-scores vs specialty peers
            CASE WHEN peer.peer_std_paid > 0
                 THEN round((ps.total_paid - peer.peer_avg_paid)
                            / peer.peer_std_paid, 2)
                 ELSE 0 END AS volume_zscore,
            CASE WHEN peer.peer_std_per_patient > 0
                 THEN round((ps.paid_per_patient - peer.peer_avg_per_patient)
                            / peer.peer_std_per_patient, 2)
                 ELSE 0 END AS intensity_zscore,
            CASE WHEN peer.peer_std_units > 0
                 THEN round((ps.avg_units - peer.peer_avg_units)
                            / peer.peer_std_units, 2)
                 ELSE 0 END AS units_zscore,
            -- Setting risk weight
            CASE ps.primary_setting
                WHEN 'office' THEN {SETTING_RISK["office"]}
                WHEN 'snf' THEN {SETTING_RISK["snf"]}
                WHEN 'asc' THEN {SETTING_RISK["asc"]}
                WHEN 'hopd' THEN {SETTING_RISK["hopd"]}
                ELSE 1.0 END AS setting_risk,
            -- Geographic risk (state per-bene spend z-score)
            (SELECT round((g.paid_per_bene - agg.m) / nullif(agg.s, 0), 2)
             FROM skin_subs.geographic_summary g,
                  (SELECT avg(paid_per_bene) AS m, stddev(paid_per_bene) AS s
                   FROM skin_subs.geographic_summary) agg
             WHERE g.state = ps.state
             LIMIT 1
            ) AS geo_risk_zscore
        FROM provider_stats ps
        LEFT JOIN peer_stats peer ON ps.provider_specialty = peer.provider_specialty
        ORDER BY volume_zscore DESC
    """)

    # Add composite risk score
    con.execute("""
        ALTER TABLE skin_subs.anomaly_scores
        ADD COLUMN composite_risk DOUBLE;

        UPDATE skin_subs.anomaly_scores
        SET composite_risk = round(
            (GREATEST(volume_zscore, 0) * 0.3
             + GREATEST(intensity_zscore, 0) * 0.3
             + GREATEST(units_zscore, 0) * 0.1
             + setting_risk * 0.2
             + GREATEST(COALESCE(geo_risk_zscore, 0), 0) * 0.1
            ), 2);
    """)

    # Add risk tier
    con.execute("""
        ALTER TABLE skin_subs.anomaly_scores
        ADD COLUMN risk_tier VARCHAR;

        UPDATE skin_subs.anomaly_scores
        SET risk_tier = CASE
            WHEN composite_risk >= 2.0 THEN 'critical'
            WHEN composite_risk >= 1.0 THEN 'high'
            WHEN composite_risk >= 0.5 THEN 'moderate'
            ELSE 'low' END;
    """)

    count = con.execute("SELECT count(*) FROM skin_subs.anomaly_scores").fetchone()[0]
    print(f"  Providers scored: {count}")

    # Tier distribution
    print("\n  Risk tier distribution:")
    for r in con.execute("""
        SELECT risk_tier, count(*) as n,
               round(avg(total_paid), 2) as avg_paid,
               round(avg(composite_risk), 2) as avg_score
        FROM skin_subs.anomaly_scores
        GROUP BY risk_tier ORDER BY avg_score DESC
    """).fetchall():
        print(
            f"    {r[0]:10s}  n={r[1]:3d}  avg_paid=${r[2]:>10,.2f}  avg_score={r[3]}"
        )

    # Top risk providers
    print("\n  Top 10 highest-risk providers:")
    for r in con.execute("""
        SELECT rendering_npi, provider_specialty, state, primary_setting,
               total_paid, paid_per_patient, composite_risk, risk_tier,
               volume_zscore, intensity_zscore
        FROM skin_subs.anomaly_scores
        ORDER BY composite_risk DESC LIMIT 10
    """).fetchall():
        print(
            f"    NPI {r[0]}  {r[1]:15s}  {r[2]}  {r[3]:8s}  "
            f"paid=${r[4]:>10,.2f}  per_pt=${r[5]:>8,.2f}  "
            f"risk={r[6]:5.2f} ({r[7]})  vol_z={r[8]:+5.2f}  int_z={r[9]:+5.2f}"
        )


def build_enforcement_correlation(con: duckdb.DuckDBPyConnection) -> None:
    """#243 — Correlate anomaly scores with known fraud patterns."""
    print("\n=== Enforcement Correlation (#243) ===")

    con.execute("DROP TABLE IF EXISTS skin_subs.provider_risk_profile")
    con.execute("""
        CREATE TABLE skin_subs.provider_risk_profile AS
        SELECT
            a.*,
            -- Match against known fraud patterns
            CASE
                WHEN a.provider_specialty = 'podiatry'
                     AND a.state = 'TX'
                     AND a.primary_setting = 'office'
                THEN 'jenson-like'
                WHEN a.primary_setting = 'snf'
                     AND a.provider_specialty IN ('nurse_practitioner', 'other')
                THEN 'gehrke-like'
                WHEN a.primary_setting = 'snf'
                     AND a.state = 'FL'
                THEN 'vohra-like'
                ELSE NULL
            END AS fraud_pattern_match,
            -- Benford analysis: first-digit distribution of paid amounts
            -- (precomputed per provider from claims)
            (SELECT round(
                -- Chi-squared statistic for Benford's law
                sum(power(observed - expected, 2) / expected), 4)
             FROM (
                SELECT
                    digit,
                    count(*) * 1.0 / nullif(total, 0) AS observed,
                    log10(1.0 + 1.0/digit) AS expected
                FROM (
                    SELECT
                        CAST(substr(CAST(CAST(round(abs(paid_amount)) AS INTEGER) AS VARCHAR), 1, 1) AS INTEGER) AS digit,
                        count(*) OVER () AS total
                    FROM skin_subs.claims_synthetic
                    WHERE rendering_npi = a.rendering_npi
                      AND paid_amount > 0
                ) digits
                WHERE digit BETWEEN 1 AND 9
                GROUP BY digit, total
             ) benford
            ) AS benford_chi2
        FROM skin_subs.anomaly_scores a
        ORDER BY composite_risk DESC
    """)

    count = con.execute(
        "SELECT count(*) FROM skin_subs.provider_risk_profile"
    ).fetchone()[0]
    print(f"  Risk profiles: {count}")

    # Fraud pattern matches
    print("\n  Fraud pattern matches:")
    for r in con.execute("""
        SELECT fraud_pattern_match, count(*) as n,
               round(avg(composite_risk), 2) as avg_risk,
               round(avg(total_paid), 2) as avg_paid
        FROM skin_subs.provider_risk_profile
        WHERE fraud_pattern_match IS NOT NULL
        GROUP BY fraud_pattern_match ORDER BY avg_risk DESC
    """).fetchall():
        print(
            f"    {r[0]:15s}  n={r[1]:3d}  avg_risk={r[2]:5.2f}  avg_paid=${r[3]:>10,.2f}"
        )

    # Do fraud-pattern providers score higher than peers?
    print("\n  Risk scores: fraud-pattern vs non-match:")
    for r in con.execute("""
        SELECT
            CASE WHEN fraud_pattern_match IS NOT NULL THEN 'matched'
                 ELSE 'unmatched' END AS group_name,
            count(*) as n,
            round(avg(composite_risk), 3) as avg_risk,
            round(avg(volume_zscore), 3) as avg_vol_z,
            round(avg(intensity_zscore), 3) as avg_int_z
        FROM skin_subs.provider_risk_profile
        GROUP BY group_name
    """).fetchall():
        print(
            f"    {r[0]:10s}  n={r[1]:3d}  risk={r[2]:6.3f}  vol_z={r[3]:+6.3f}  int_z={r[4]:+6.3f}"
        )

    # Benford analysis summary
    print(
        "\n  Benford's law chi² (lower = more conformant, >15.5 = suspicious at p<0.05):"
    )
    for r in con.execute("""
        SELECT risk_tier, count(*) as n,
               round(avg(benford_chi2), 4) as avg_chi2,
               count(*) FILTER (WHERE benford_chi2 > 15.5) as suspicious
        FROM skin_subs.provider_risk_profile
        GROUP BY risk_tier ORDER BY avg_chi2 DESC
    """).fetchall():
        print(f"    {r[0]:10s}  n={r[1]:3d}  avg_chi²={r[2]:8.4f}  suspicious={r[3]}")


def main() -> None:
    print("Building anomaly scoring and risk profiles ...")
    con = duckdb.connect(str(DUCKDB_PATH))

    build_anomaly_scores(con)
    build_enforcement_correlation(con)

    print("\n=== Final table count ===")
    for r in con.execute("""
        SELECT table_name FROM information_schema.tables
        WHERE table_schema = 'skin_subs' ORDER BY table_name
    """).fetchall():
        cnt = con.execute(f"SELECT count(*) FROM skin_subs.{r[0]}").fetchone()[0]
        print(f"  skin_subs.{r[0]:30s}: {cnt:>6} rows")

    con.close()
    print("\nDone.")


if __name__ == "__main__":
    main()