Districts with Similar Health Needs

# ---------------------------------------------------------------------------
# district_health.py
#
# Capstone project: grouping districts with similar health needs, so that
# each group gets the programme it needs instead of one plan for all.
#
# CBSE Class XII, Artificial Intelligence (843), Part C - Capstone Project.
# SDG 3 (Good Health and Well-being) and SDG 10 (Reduced Inequalities).
#
# Data Science Methodology, 10 steps. UNSUPERVISED LEARNING: the districts
# have no labels. K-means clustering finds the groups; the number of groups
# is chosen with the elbow method and the silhouette score, and the result
# is checked for stability. Each group is then described by its weakest
# indicators.
#
# Libraries: pandas, numpy, matplotlib, scikit-learn, joblib.
# ---------------------------------------------------------------------------

# %% Step 1 - Business Understanding
print("=" * 64)
print("DISTRICTS WITH SIMILAR HEALTH NEEDS")
print("=" * 64)
print("Problem : one state-wide plan treats every district the same,")
print("          though their health needs are very different")
print("Goal    : find groups of districts with similar needs, and what")
print("          each group needs most")
print()

# %% Step 2 - Analytic Approach
print("Analytic approach : descriptive - clustering (unsupervised learning)")
print("There are no labels: nobody has said which district belongs where")
print()

# %% Step 3 - Data Requirements
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, silhouette_samples, adjusted_rand_score
import joblib

INDICATORS = ["institutional_births_pct", "full_immunisation_pct", "stunted_children_pct", "anaemic_women_pct",
              "improved_sanitation_pct", "clean_fuel_pct", "women_literacy_pct", "health_insurance_pct"]
SHORT = ["Births in\nhospital", "Full\nimmunisation", "Stunted\nchildren", "Anaemic\nwomen",
         "Sanitation", "Clean\nfuel", "Women's\nliteracy", "Health\ninsurance"]
LOWER_IS_BETTER = ["stunted_children_pct", "anaemic_women_pct"]

# %% Step 4 - Data Collection
data = pd.read_csv("districts.csv")
print("Districts collected :", len(data))
print("Indicators          :", len(INDICATORS))
print()

# %% Step 5 - Data Understanding
print("--- Step 5: Data Understanding ---")
print("Blank cells in each indicator:")
print(data[INDICATORS].isnull().sum().to_string())
print("Values above 100 percent :", int((data[INDICATORS] > 100).sum().sum()))
print()

# Chart 1: how the indicators move together.
corr = data[INDICATORS].corr()
plt.figure(figsize=(8, 6.5))
plt.imshow(corr, cmap="RdBu", vmin=-1, vmax=1)
for a in range(len(INDICATORS)):
    for b in range(len(INDICATORS)):
        plt.text(b, a, round(corr.iloc[a, b], 2), ha="center", va="center", fontsize=8)
plt.xticks(range(len(INDICATORS)), SHORT, fontsize=8)
plt.yticks(range(len(INDICATORS)), SHORT, fontsize=8)
plt.colorbar(label="Correlation")
plt.title("How the indicators move together")
plt.tight_layout()
plt.savefig("chart1_correlation.png")
plt.close()

# %% Step 6 - Data Preparation
print("--- Step 6: Data Preparation ---")
before = len(data)
data = data.drop_duplicates(subset="district").reset_index(drop=True)
print("District entered twice, removed :", before - len(data))
typo = data[INDICATORS] > 100
print("Values typed without the decimal point, fixed :", int(typo.sum().sum()))
for col in INDICATORS:
    data.loc[data[col] > 100, col] = data[col] / 10        # 912 -> 91.2
medians = data[INDICATORS].median()
print("Blank cells filled with the indicator's median :", int(data[INDICATORS].isnull().sum().sum()))
data[INDICATORS] = data[INDICATORS].fillna(medians)
print()
print("Spread (standard deviation) of each indicator:")
print(data[INDICATORS].std().round(1).to_string())
# K-means measures distance. Without scaling, the indicators that vary most
# (clean fuel, insurance) would decide the groups on their own.
scaler = StandardScaler()
X = scaler.fit_transform(data[INDICATORS])
print("Every indicator rescaled to mean 0 and standard deviation 1")
print()

# Chart 2: the spread of every indicator.
plt.figure(figsize=(10, 4.5))
plt.boxplot([data[c] for c in INDICATORS])
plt.xticks(range(1, len(INDICATORS) + 1), SHORT, fontsize=8)
plt.ylabel("Percent")
plt.title("The eight indicators across all districts")
plt.tight_layout()
plt.savefig("chart2_spread.png")
plt.close()

# %% Step 7 - Modelling
print("--- Step 7: Modelling ---")
print("k   inertia  silhouette")
inertias = []
silhouettes = []
ks = range(2, 9)
for k in ks:
    km = KMeans(n_clusters=k, n_init=10, random_state=0).fit(X)
    inertias.append(km.inertia_)
    silhouettes.append(silhouette_score(X, km.labels_))
    print(k, " ", round(km.inertia_, 1), "  ", round(silhouettes[-1], 4))
# When two values of k have almost the same silhouette, the larger one is
# kept: more groups mean more targeted programmes, at no cost in fit.
best = max(silhouettes)
K = max(k for k, s in zip(ks, silhouettes) if s >= best - 0.01)
print("Highest silhouette :", round(best, 4))
print("Number of groups chosen :", K)
print("(the most groups with a silhouette within 0.01 of the highest)")
model = KMeans(n_clusters=K, n_init=10, random_state=0).fit(X)
print()

# Chart 3: the elbow and the silhouette.
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(list(ks), inertias, marker="o", color="#4e79a7")
plt.axvline(K, color="gray", linestyle="--")
plt.title("Elbow method: inertia")
plt.xlabel("Number of groups (k)")
plt.subplot(1, 2, 2)
plt.plot(list(ks), silhouettes, marker="o", color="#e15759")
plt.axvline(K, color="gray", linestyle="--")
plt.title("Silhouette score")
plt.xlabel("Number of groups (k)")
plt.tight_layout()
plt.savefig("chart3_elbow_silhouette.png")
plt.close()

# Name the groups in order of need: Group 1 has the greatest need.
sign = np.array([-1 if c in LOWER_IS_BETTER else 1 for c in INDICATORS])
wellbeing = (model.cluster_centers_ * sign).mean(axis=1)    # high = better off
order = np.argsort(wellbeing)
group_of = {old: new + 1 for new, old in enumerate(order)}
data["group"] = [group_of[g] for g in model.labels_]

# %% Step 8 - Evaluation
print("--- Step 8: Evaluation ---")
print("Silhouette score :", round(silhouette_score(X, model.labels_), 3))
same = []
for seed in range(1, 11):
    other = KMeans(n_clusters=K, n_init=10, random_state=seed).fit(X)
    same.append(adjusted_rand_score(model.labels_, other.labels_))
print("Lowest agreement with 10 runs from other starting points (adjusted Rand index) :", round(min(same), 3))
print("Districts in each group:")
print(data["group"].value_counts().sort_index().to_string())
print("Group sizes :", ", ".join(str(n) for n in data["group"].value_counts().sort_index()))
print()
profile = data.groupby("group")[INDICATORS].mean().round(1)
print("Average of each indicator in each group:")
print(profile.T.to_string())
print()
state = data[INDICATORS].mean()
need = {}
for g in profile.index:
    gap = (profile.loc[g] - state) / data[INDICATORS].std() * sign   # negative = worse than average
    weakest = gap.sort_values().index[:2]
    need[g] = [w for w in weakest if gap[w] < 0]
    names = [SHORT[INDICATORS.index(w)].replace("\n", " ").lower() for w in need[g]]
    print("Group", g, "- weakest compared with the state :", ", ".join(names) if names else "none - above average on all")
sil = silhouette_samples(X, model.labels_)
data["silhouette"] = sil
print("Districts that fit their group poorly (silhouette below 0) :", int((sil < 0).sum()))
print()

# Chart 4: each group's profile.
z = (profile - state) / data[INDICATORS].std() * sign
plt.figure(figsize=(10, 3.2 + 0.4 * K))
plt.imshow(z.values, cmap="RdYlGn", vmin=-2, vmax=2, aspect="auto")
for a in range(K):
    for b in range(len(INDICATORS)):
        plt.text(b, a, profile.iloc[a, b], ha="center", va="center", fontsize=9)
plt.xticks(range(len(INDICATORS)), SHORT, fontsize=8)
plt.yticks(range(K), ["Group " + str(g) for g in profile.index])
plt.colorbar(label="Better (green) or worse (red) than the state average")
plt.title("What each group looks like (numbers are group averages, %)")
plt.tight_layout()
plt.savefig("chart4_profiles.png")
plt.close()

# Chart 5: the groups on two of the indicators.
colours = ["#e15759", "#f28e2b", "#edc948", "#59a14f", "#4e79a7", "#b07aa1", "#9c755f", "#bab0ac"]
plt.figure(figsize=(7.5, 5))
for g in range(1, K + 1):
    part = data[data["group"] == g]
    plt.scatter(part["clean_fuel_pct"], part["stunted_children_pct"], s=18, color=colours[g - 1], label="Group " + str(g))
plt.xlabel("Households using clean cooking fuel, %")
plt.ylabel("Stunted children, %")
plt.title("The districts on two of the eight indicators")
plt.legend()
plt.tight_layout()
plt.savefig("chart5_scatter.png")
plt.close()

# Chart 6: how well each district fits its group.
plt.figure(figsize=(8, 4.5))
y0 = 0
for g in range(1, K + 1):
    s = np.sort(data.loc[data["group"] == g, "silhouette"].values)
    plt.barh(range(y0, y0 + len(s)), s, height=1.0, color=colours[g - 1], label="Group " + str(g))
    y0 = y0 + len(s) + 5
plt.axvline(0, color="black", linewidth=0.8)
plt.title("Silhouette of every district, by group")
plt.xlabel("Silhouette (near 1 = fits its group well, below 0 = closer to another group)")
plt.yticks([])
plt.legend(loc="lower right")
plt.tight_layout()
plt.savefig("chart6_silhouette.png")
plt.close()

# %% Step 9 - Deployment
print("--- Step 9: Deployment ---")
joblib.dump({"medians": medians, "scaler": scaler, "model": model, "group_of": group_of}, "district_groups.pkl")
saved = joblib.load("district_groups.pkl")


def group_for(values):
    row = pd.DataFrame([values], columns=INDICATORS).fillna(saved["medians"])
    label = saved["model"].predict(saved["scaler"].transform(row))[0]
    return saved["group_of"][label]


new_district = [78, 64, 36, 60, 64, 33, 66, None]
g = group_for(new_district)
print("A newly formed district, insurance figure not yet known : Group", g)
print("Programme focus for Group", g, ":", ", ".join(SHORT[INDICATORS.index(w)].replace("\n", " ").lower() for w in need[g]))
print()

# %% Step 10 - Feedback
print("--- Step 10: Feedback ---")
nxt = pd.read_csv("districts_next_survey.csv")
nxt["group"] = [group_for(list(r)) for r in nxt[INDICATORS].values]
both = data[["district", "group"]].merge(nxt[["district", "group"]], on="district", suffixes=("_before", "_after"))
moved = both[both["group_before"] != both["group_after"]]
print("Districts in the next survey        :", len(nxt))
print("Districts that changed group        :", len(moved))
print("  moved to a better-off group       :", int((moved["group_after"] > moved["group_before"]).sum()))
print("  moved to a group in greater need  :", int((moved["group_after"] < moved["group_before"]).sum()))
print("Districts in Group 1 before         :", int((both["group_before"] == 1).sum()))
print("Districts in Group 1 after          :", int((both["group_after"] == 1).sum()))
print("Districts that jumped two or more groups:")
jump = moved[moved["group_after"] - moved["group_before"] >= 2]
for _, r in jump.iterrows():
    print("  ", r["district"], ": Group", r["group_before"], "-> Group", r["group_after"])

# Chart 7: where each district went.
table = pd.crosstab(both["group_before"], both["group_after"]).reindex(index=range(1, K + 1), columns=range(1, K + 1), fill_value=0)
plt.figure(figsize=(6, 5))
plt.imshow(table.values, cmap="Blues")
for a in range(K):
    for b in range(K):
        plt.text(b, a, table.iloc[a, b], ha="center", va="center",
                 color="white" if table.iloc[a, b] > table.values.max() / 2 else "black")
plt.xticks(range(K), ["Group " + str(g) for g in range(1, K + 1)])
plt.yticks(range(K), ["Group " + str(g) for g in range(1, K + 1)])
plt.xlabel("Group in the next survey")
plt.ylabel("Group before")
plt.title("Feedback: which group each district is in now")
plt.tight_layout()
plt.savefig("chart7_next_survey.png")
plt.close()

================================================================
DISTRICTS WITH SIMILAR HEALTH NEEDS
================================================================
Problem : one state-wide plan treats every district the same,
          though their health needs are very different
Goal    : find groups of districts with similar needs, and what
          each group needs most

Analytic approach : descriptive - clustering (unsupervised learning)
There are no labels: nobody has said which district belongs where

Districts collected : 301
Indicators          : 8

--- Step 5: Data Understanding ---
Blank cells in each indicator:
institutional_births_pct    0
full_immunisation_pct       4
stunted_children_pct        0
anaemic_women_pct           0
improved_sanitation_pct     0
clean_fuel_pct              4
women_literacy_pct          0
health_insurance_pct        5
Values above 100 percent : 3

--- Step 6: Data Preparation ---
District entered twice, removed : 1
Values typed without the decimal point, fixed : 3
Blank cells filled with the indicator's median : 12

Spread (standard deviation) of each indicator:
institutional_births_pct    12.4
full_immunisation_pct       13.1
stunted_children_pct         8.6
anaemic_women_pct            8.7
improved_sanitation_pct     15.9
clean_fuel_pct              21.3
women_literacy_pct          13.3
health_insurance_pct        15.6
Every indicator rescaled to mean 0 and standard deviation 1

--- Step 7: Modelling ---
k   inertia  silhouette
2   1197.7    0.3996
3   743.4    0.4115
4   514.9    0.4109
5   491.2    0.3383
6   476.6    0.2407
7   453.1    0.1968
8   437.8    0.1291
Highest silhouette : 0.4115
Number of groups chosen : 4
(the most groups with a silhouette within 0.01 of the highest)

--- Step 8: Evaluation ---
Silhouette score : 0.411
Lowest agreement with 10 runs from other starting points (adjusted Rand index) : 1.0
Districts in each group:
group
1    70
2    84
3    86
4    60
Group sizes : 70, 84, 86, 60

Average of each indicator in each group:
group                        1     2     3     4
institutional_births_pct  64.3  89.0  84.2  94.6
full_immunisation_pct     55.8  79.6  71.1  89.9
stunted_children_pct      41.2  39.5  30.3  21.3
anaemic_women_pct         63.8  61.3  51.7  45.5
improved_sanitation_pct   51.2  68.1  82.4  91.1
clean_fuel_pct            24.2  37.5  55.6  80.2
women_literacy_pct        55.4  71.0  78.6  90.3
health_insurance_pct      20.4  33.2  46.2  59.5

Group 1 - weakest compared with the state : births in hospital, full immunisation
Group 2 - weakest compared with the state : stunted children, anaemic women
Group 3 - weakest compared with the state : full immunisation
Group 4 - weakest compared with the state : none - above average on all
Districts that fit their group poorly (silhouette below 0) : 0

--- Step 9: Deployment ---
A newly formed district, insurance figure not yet known : Group 2
Programme focus for Group 2 : stunted children, anaemic women

--- Step 10: Feedback ---
Districts in the next survey        : 300
Districts that changed group        : 32
  moved to a better-off group       : 32
  moved to a group in greater need  : 0
Districts in Group 1 before         : 70
Districts in Group 1 after          : 56
Districts that jumped two or more groups:
   District 118 : Group 1 -> Group 3
   District 392 : Group 1 -> Group 3