Tribal Health Analytics¶

# Databricks notebook source
# MAGIC %md
# MAGIC # Silver → Gold: Health Disparity Index Calculation
# MAGIC Calculates composite health disparity indices comparing AI/AN
# MAGIC age-adjusted rates to national benchmarks using CDC public data.
# MAGIC
# MAGIC **DATA NOTE:** This analysis uses ONLY publicly available aggregate
# MAGIC federal data. No tribal-specific or patient-level data is accessed.

# COMMAND ----------

from pyspark.sql import functions as F

# Read Silver CDC WONDER mortality data (age-adjusted rates)
df_mortality = spark.read.format("delta").load(
    "Tables/silver_tribal_cdc_mortality_rates"
)

# Read Silver BRFSS prevalence data
df_brfss = spark.read.format("delta").load(
    "Tables/silver_tribal_brfss_prevalence"
)

# Define health disparity domains and national benchmarks
# Source: CDC WONDER, age-adjusted rates per 100,000 (2020-2022)
national_benchmarks = {
    "all_cause_mortality": 835.4,
    "diabetes_mortality": 25.4,
    "chronic_liver_disease": 12.8,
    "heart_disease_mortality": 173.8,
    "suicide_mortality": 14.0,
    "unintentional_injury": 57.6,
}

# Calculate disparity ratios (AI/AN rate / national rate)
for condition, benchmark in national_benchmarks.items():
    df_mortality = df_mortality.withColumn(
        f"{condition}_disparity_ratio",
        F.when(
            F.col("race_ethnicity") == "American Indian or Alaska Native",
            F.col(f"{condition}_rate") / F.lit(benchmark)
        )
    )

# Aggregate to IHS Area level
disparity_cols = [f"{c}_disparity_ratio" for c in national_benchmarks.keys()]

df_hdi = (
    df_mortality
    .filter(F.col("race_ethnicity") == "American Indian or Alaska Native")
    .groupBy("ihs_area", "state", "reporting_year")
    .agg(
        *[F.avg(c).alias(c) for c in disparity_cols],
        F.count("*").alias("data_points"),
    )
    .withColumn(
        "composite_hdi",
        F.round(
            (F.col("all_cause_mortality_disparity_ratio") * 0.25 +
             F.col("diabetes_mortality_disparity_ratio") * 0.20 +
             F.col("heart_disease_mortality_disparity_ratio") * 0.15 +
             F.col("chronic_liver_disease_disparity_ratio") * 0.15 +
             F.col("suicide_mortality_disparity_ratio") * 0.15 +
             F.col("unintentional_injury_disparity_ratio") * 0.10),
            3
        )
    )
    .withColumn(
        "disparity_level",
        F.when(F.col("composite_hdi") >= 2.0, "SEVERE")
        .when(F.col("composite_hdi") >= 1.5, "SIGNIFICANT")
        .when(F.col("composite_hdi") >= 1.2, "MODERATE")
        .otherwise("NEAR_PARITY")
    )
)

df_hdi.write.format("delta").mode("overwrite").option(
    "overwriteSchema", "true"
).saveAsTable("lh_gold.gold_tribal_health_disparity_index")

severe = df_hdi.filter(F.col("disparity_level") == "SEVERE").count()
print(f"IHS Area-years with SEVERE health disparities: {severe}")

# Databricks notebook source
# MAGIC %md
# MAGIC # Gold: Chronic Disease Prevalence Trends by IHS Area
# MAGIC Models multi-year chronic disease prevalence using BRFSS data
# MAGIC to identify accelerating disease burden in AI/AN service areas.
# MAGIC
# MAGIC **DATA NOTE:** BRFSS provides aggregate prevalence estimates.
# MAGIC No individual-level or tribal-specific records are accessed.

# COMMAND ----------

from pyspark.sql import functions as F
from pyspark.sql.window import Window

# Read Silver BRFSS prevalence data (AI/AN filtered, multi-year pooled)
df_brfss = spark.read.format("delta").load(
    "Tables/silver_tribal_brfss_prevalence"
)

# Focus on key chronic conditions
conditions = [
    "diabetes_prevalence_pct",
    "heart_disease_prevalence_pct",
    "obesity_prevalence_pct",
    "smoking_prevalence_pct",
    "heavy_drinking_prevalence_pct",
    "no_health_insurance_pct",
]

# Calculate 3-year rolling averages to smooth small sample effects
area_window = Window.partitionBy("ihs_area", "condition").orderBy("survey_year")
rolling_window = Window.partitionBy("ihs_area", "condition").orderBy("survey_year").rowsBetween(-2, 0)

# Unpivot conditions for trend analysis
df_unpivoted = df_brfss.selectExpr(
    "ihs_area", "state", "survey_year",
    "stack(6, "
    "'diabetes', diabetes_prevalence_pct, "
    "'heart_disease', heart_disease_prevalence_pct, "
    "'obesity', obesity_prevalence_pct, "
    "'smoking', smoking_prevalence_pct, "
    "'heavy_drinking', heavy_drinking_prevalence_pct, "
    "'uninsured', no_health_insurance_pct"
    ") as (condition, prevalence_pct)"
)

df_trends = (
    df_unpivoted
    .withColumn("rolling_avg_3yr", F.avg("prevalence_pct").over(rolling_window))
    .withColumn("prev_year_avg", F.lag("rolling_avg_3yr").over(area_window))
    .withColumn(
        "trend_direction",
        F.when(F.col("rolling_avg_3yr") > F.col("prev_year_avg") * 1.02, "INCREASING")
        .when(F.col("rolling_avg_3yr") < F.col("prev_year_avg") * 0.98, "DECREASING")
        .otherwise("STABLE")
    )
    .withColumn(
        "acceleration_flag",
        F.when(
            (F.col("trend_direction") == "INCREASING") &
            (F.col("prevalence_pct") > F.col("rolling_avg_3yr")),
            "ACCELERATING"
        ).otherwise("NORMAL")
    )
)

df_trends.write.format("delta").mode("overwrite").option(
    "overwriteSchema", "true"
).saveAsTable("lh_gold.gold_tribal_chronic_disease_trends")

accelerating = df_trends.filter(F.col("acceleration_flag") == "ACCELERATING").count()
print(f"Condition-area-years with accelerating prevalence: {accelerating}")

# Databricks notebook source
# MAGIC %md
# MAGIC # Gold: IHS Facility Utilization Benchmarking
# MAGIC Calculates facility utilization rates and identifies capacity
# MAGIC gaps using publicly available IHS workload statistics.
# MAGIC
# MAGIC **DATA NOTE:** Uses only aggregate IHS published workload data.
# MAGIC No facility-specific patient data or tribal records accessed.

# COMMAND ----------

from pyspark.sql import functions as F

# Read Silver IHS workload data (from published IHS fact sheets)
df_ihs = spark.read.format("delta").load(
    "Tables/silver_tribal_ihs_workload"
)

# Read HRSA health workforce data for provider supply context
df_hrsa = spark.read.format("delta").load(
    "Tables/silver_tribal_hrsa_workforce"
)

# Calculate utilization metrics per IHS Area
df_utilization = (
    df_ihs
    .groupBy("ihs_area", "fiscal_year")
    .agg(
        F.sum("user_population").alias("total_users"),
        F.sum("outpatient_visits").alias("total_op_visits"),
        F.sum("inpatient_days").alias("total_ip_days"),
        F.sum("dental_visits").alias("total_dental_visits"),
        F.sum("community_health_contacts").alias("total_ch_contacts"),
        F.sum("er_visits").alias("total_er_visits"),
        F.countDistinct("facility_id").alias("facility_count"),
    )
    .withColumn(
        "op_visits_per_user",
        F.round(F.col("total_op_visits") / F.col("total_users"), 2)
    )
    .withColumn(
        "dental_visits_per_user",
        F.round(F.col("total_dental_visits") / F.col("total_users"), 3)
    )
    .withColumn(
        "er_rate_per_1000",
        F.round(F.col("total_er_visits") / F.col("total_users") * 1000, 1)
    )
)

# Benchmark against national averages (AHRQ MEPS data)
# National avg outpatient visits: ~4.0/person/year
# National avg dental visits: ~0.44/person/year
# National avg ER visits: ~429 per 1,000
df_benchmarked = (
    df_utilization
    .withColumn("op_benchmark_ratio",
        F.round(F.col("op_visits_per_user") / 4.0, 2))
    .withColumn("dental_benchmark_ratio",
        F.round(F.col("dental_visits_per_user") / 0.44, 2))
    .withColumn("er_benchmark_ratio",
        F.round(F.col("er_rate_per_1000") / 429, 2))
    .withColumn(
        "capacity_assessment",
        F.when(
            (F.col("op_benchmark_ratio") < 0.7) &
            (F.col("er_benchmark_ratio") > 1.3),
            "CAPACITY_GAP"  # Low outpatient + high ER = access problem
        )
        .when(F.col("op_benchmark_ratio") < 0.7, "UNDERSERVED")
        .when(F.col("er_benchmark_ratio") > 1.5, "ER_OVERRELIANCE")
        .otherwise("ADEQUATE")
    )
)

df_benchmarked.write.format("delta").mode("overwrite").option(
    "overwriteSchema", "true"
).saveAsTable("lh_gold.gold_tribal_facility_utilization")

gaps = df_benchmarked.filter(F.col("capacity_assessment") == "CAPACITY_GAP").count()
print(f"IHS Area-years with capacity gaps: {gaps}")

# Cross-domain: USDA food access vs. diabetes prevalence in AI/AN areas
df_cross_food = (
    df_usda_food_access
    .filter(F.col("tract_overlaps_tribal_area") == True)
    .join(df_tribal_chronic_disease, on=["state", "county_fips"])
    .filter(F.col("condition") == "diabetes")
    .groupBy("ihs_area", "state")
    .agg(
        F.avg("low_access_pct").alias("avg_food_desert_pct"),
        F.avg("prevalence_pct").alias("avg_diabetes_prevalence"),
        F.count("*").alias("tracts_analyzed"),
    )
    .withColumn(
        "food_diabetes_correlation",
        F.corr("avg_food_desert_pct", "avg_diabetes_prevalence")
            .over(Window.partitionBy(F.lit(1)))
    )
)

# Cross-domain: EPA environmental burden vs. tribal health disparities
df_cross_ej = (
    df_epa_ejscreen
    .filter(F.col("pct_aian") > 0.10)  # Tracts with >10% AI/AN population
    .join(df_tribal_health_disparity, on=["state_fips", "county_fips"])
    .select(
        "state", "county", "ihs_area",
        "pm25_score", "ozone_score", "lead_paint_score",
        "superfund_proximity_score", "rmp_proximity_score",
        "composite_hdi", "disparity_level",
    )
    .withColumn(
        "environmental_health_risk",
        F.when(
            (F.col("superfund_proximity_score") > 80) &
            (F.col("composite_hdi") >= 1.5),
            "HIGH_RISK_OVERLAP"
        ).otherwise("STANDARD")
    )
)