stack/dev/scripts/build_product_susceptibility.py

"""Systematic product susceptibility scoring for kickback/fraud risk.

Scores each skin substitute product on structural factors that make
it more or less susceptible to kickback arrangements and fraudulent
billing schemes.

Addresses #260.

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

from __future__ import annotations

from pathlib import Path

import duckdb

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

FLAT_RATE = 127.28


def main() -> None:
    print("Building product susceptibility scores ...")
    con = duckdb.connect(str(DUCKDB_PATH))

    con.execute("DROP TABLE IF EXISTS skin_subs.product_susceptibility")
    con.execute("""
        CREATE TABLE skin_subs.product_susceptibility AS
        WITH product_base AS (
            SELECT
                p.hcpcs_code,
                p.product_name,
                p.manufacturer,
                p.category,
                p.fda_pathway,
                p.status,
                p.latest_asp,
                p.latest_payment_limit,
                p.quarters_on_market,
                p.first_quarter,
                p.avg_asp,
                p.max_asp,
                p.asp_range,
                p.flat_rate_delta,
                p.claim_lines,
                p.total_paid,
                p.unique_providers,
                p.unique_benes
            FROM skin_subs.products_enriched p
        ),
        -- Evidence quality: count of studies mentioning this product
        evidence AS (
            SELECT
                products AS product_key,
                count(*) AS study_count,
                count(*) FILTER (WHERE pub_type = 'rct') AS rct_count,
                count(*) FILTER (WHERE pub_type = 'meta-analysis') AS meta_count,
                count(*) FILTER (WHERE is_industry_linked) AS industry_funded,
                count(*) FILTER (WHERE stance = 'skeptical') AS skeptical_studies
            FROM skin_subs.study_characteristics
            WHERE products != ''
            GROUP BY products
        ),
        -- Anomalous price trajectory
        price_anomaly AS (
            SELECT
                hcpcs_code,
                count(*) FILTER (WHERE anomalous_spike) AS spike_count,
                max(abs(qoq_pct_change)) AS max_qoq_change
            FROM skin_subs.asp_trajectories
            GROUP BY hcpcs_code
        ),
        -- Provider concentration: how many providers use this product
        provider_profile AS (
            SELECT
                hcpcs_code,
                count(DISTINCT rendering_npi) AS provider_count,
                -- Is it used mostly by high-risk providers?
                count(DISTINCT rendering_npi) FILTER (
                    WHERE rendering_npi IN (
                        SELECT rendering_npi FROM skin_subs.anomaly_scores
                        WHERE risk_tier IN ('critical', 'high')
                    )
                ) AS high_risk_provider_count,
                -- Setting mix
                count(*) FILTER (WHERE place_of_service = '11') * 100.0
                    / nullif(count(*), 0) AS pct_office,
                count(*) FILTER (WHERE place_of_service IN ('31','32')) * 100.0
                    / nullif(count(*), 0) AS pct_snf
            FROM skin_subs.claims_synthetic
            WHERE claim_type = 'product'
            GROUP BY hcpcs_code
        )
        SELECT
            pb.*,

            -- === PRICE-BASED SUSCEPTIBILITY ===
            -- Higher ASP = larger kickback headroom (ASP+6% markup)
            CASE WHEN pb.latest_asp IS NOT NULL
                 THEN round(pb.latest_asp * 0.06, 2)
                 ELSE 0 END AS kickback_headroom_per_cm2,

            -- Flat rate delta: products losing most had most to protect
            COALESCE(pb.flat_rate_delta, 0) AS reclassification_loss,

            -- Price volatility score (0-1): normalized spike count + range
            round(COALESCE(pa.spike_count, 0) * 0.3
                  + LEAST(COALESCE(pa.max_qoq_change, 0) / 100.0, 1.0) * 0.7,
                  3) AS price_volatility_score,

            -- === EVIDENCE-BASED SUSCEPTIBILITY ===
            -- Low evidence = higher susceptibility
            COALESCE(ev.study_count, 0) AS study_count,
            COALESCE(ev.rct_count, 0) AS rct_count,
            COALESCE(ev.industry_funded, 0) AS industry_funded_studies,
            -- Price-to-evidence ratio: high price with weak evidence
            CASE WHEN COALESCE(ev.rct_count, 0) > 0
                 THEN round(COALESCE(pb.latest_asp, 0) / ev.rct_count, 2)
                 ELSE COALESCE(pb.latest_asp, 0)
                 END AS price_per_rct,

            -- === DISTRIBUTION-BASED SUSCEPTIBILITY ===
            COALESCE(pp.provider_count, 0) AS active_providers,
            COALESCE(pp.high_risk_provider_count, 0) AS high_risk_providers,
            round(COALESCE(pp.pct_office, 0), 1) AS pct_office_setting,
            round(COALESCE(pp.pct_snf, 0), 1) AS pct_snf_setting,

            -- === REGULATORY-BASED SUSCEPTIBILITY ===
            -- FDA pathway risk: 510(k) < PMA < HCT/P (less evidence required)
            CASE pb.fda_pathway
                WHEN 'PMA' THEN 0.3
                WHEN '510(k)' THEN 0.5
                WHEN 'HCT/P' THEN 0.8
                ELSE 0.7 END AS fda_pathway_risk,

            -- === COMPOSITE SUSCEPTIBILITY SCORE ===
            round(
                -- Price component (30%): normalized ASP
                LEAST(COALESCE(pb.latest_asp, 0) / 1000.0, 1.0) * 0.15
                + LEAST(COALESCE(pb.flat_rate_delta, 0) / 3000.0, 1.0) * 0.15

                -- Evidence component (25%): inverse evidence quality
                + CASE WHEN COALESCE(ev.rct_count, 0) = 0 THEN 0.25
                       WHEN ev.rct_count <= 2 THEN 0.15
                       ELSE 0.05 END

                -- Distribution component (25%): office/SNF + high-risk providers
                + LEAST(COALESCE(pp.pct_office, 0) / 100.0, 1.0) * 0.10
                + LEAST(COALESCE(pp.pct_snf, 0) / 100.0, 1.0) * 0.05
                + CASE WHEN COALESCE(pp.high_risk_provider_count, 0) > 0
                       THEN 0.10 ELSE 0.0 END

                -- Price volatility (10%)
                + COALESCE(pa.spike_count, 0) * 0.02

                -- FDA pathway (10%)
                + CASE pb.fda_pathway
                    WHEN 'PMA' THEN 0.03
                    WHEN '510(k)' THEN 0.05
                    WHEN 'HCT/P' THEN 0.08
                    ELSE 0.07 END
            , 3) AS susceptibility_score,

            -- Risk tier
            'placeholder' AS susceptibility_tier

        FROM product_base pb
        LEFT JOIN evidence ev ON pb.product_name = ev.product_key
        LEFT JOIN price_anomaly pa ON pb.hcpcs_code = pa.hcpcs_code
        LEFT JOIN provider_profile pp ON pb.hcpcs_code = pp.hcpcs_code
        ORDER BY susceptibility_score DESC NULLS LAST
    """)

    # Update tiers based on score
    con.execute("""
        UPDATE skin_subs.product_susceptibility
        SET susceptibility_tier = CASE
            WHEN susceptibility_score >= 0.60 THEN 'critical'
            WHEN susceptibility_score >= 0.45 THEN 'high'
            WHEN susceptibility_score >= 0.30 THEN 'moderate'
            ELSE 'low' END
    """)

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

    # Tier distribution
    print("\n  Susceptibility tier distribution:")
    for r in con.execute("""
        SELECT susceptibility_tier, count(*) as n,
               round(avg(susceptibility_score), 3) as avg_score,
               round(avg(latest_asp), 2) as avg_asp
        FROM skin_subs.product_susceptibility
        GROUP BY susceptibility_tier ORDER BY avg_score DESC
    """).fetchall():
        print(
            f"    {r[0]:10s}  n={r[1]:3d}  avg_score={r[2]:6.3f}  avg_asp=${r[3] or 0:>8}"
        )

    # Top 20 most susceptible
    print("\n  Top 20 most susceptible products:")
    for r in con.execute("""
        SELECT hcpcs_code, product_name, manufacturer, category,
               susceptibility_score, susceptibility_tier,
               latest_asp, rct_count, active_providers, pct_office_setting
        FROM skin_subs.product_susceptibility
        WHERE susceptibility_score IS NOT NULL
        ORDER BY susceptibility_score DESC LIMIT 20
    """).fetchall():
        print(
            f"    {r[0]}  {(r[1] or ''):25s}  {(r[2] or ''):20s}  "
            f"score={r[4]:5.3f} ({r[5]})  asp=${r[6] or 0:>8}  "
            f"rcts={r[7]:2d}  provs={r[8]:3d}  office={r[9] or 0}%"
        )

    # Cross-reference: do susceptible products appear in enforcement?
    print("\n  Susceptibility by category:")
    for r in con.execute("""
        SELECT category, count(*) as n,
               round(avg(susceptibility_score), 3) as avg_score,
               round(avg(latest_asp), 2) as avg_asp,
               sum(rct_count) as total_rcts
        FROM skin_subs.product_susceptibility
        WHERE category IS NOT NULL
        GROUP BY category
        HAVING count(*) >= 2
        ORDER BY avg_score DESC
    """).fetchall():
        print(
            f"    {(r[0] or ''):40s}  n={r[1]:3d}  score={r[2]:6.3f}  "
            f"asp=${r[3] or 0:>8}  rcts={r[4]:3d}"
        )

    # Final table count
    print("\n  All skin_subs tables:")
    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()