03B — Model Evolution: V1 to V2

Mengapa enam cluster disederhanakan menjadi tiga tipologi

Notebook ini adalah POC portfolio yang dapat dijalankan ulang dari artefak lokal Pasuruan Lens. Kode produksi tetap berada di python/scripts/ dan python/src/pasuruan365/; notebook berfungsi sebagai narasi analitik yang ringkas, transparan, dan mudah dipresentasikan.

from pathlib import Path
import json
import sys

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns


def find_project_root(start: Path | None = None) -> Path:
    start = (start or Path.cwd()).resolve()
    for candidate in (start, *start.parents):
        if (candidate / "package.json").exists() and (candidate / "data").exists():
            return candidate
    raise RuntimeError("Root Pasuruan Lens tidak ditemukan. Jalankan notebook dari repository ini.")


ROOT = find_project_root()
sys.path.insert(0, str(ROOT / "python" / "src"))

pd.set_option("display.max_columns", 100)
pd.set_option("display.max_colwidth", 100)
sns.set_theme(style="whitegrid", context="notebook")
COLORS = {"A": "#2563eb", "B": "#f59e0b", "C": "#10b981"}

print(f"Project root: {ROOT}")
print(f"Python: {sys.version.split()[0]}")

Project root: PASURUAN365
Python: 3.12.6

Executive takeaway

V1 (k=6) adalah baseline yang valid dan cukup stabil ketika pipeline yang sama diulang. Namun, robustness review menemukan dua masalah: assignment sangat berubah ketika scaler diganti dan fitur poliklinik berbentuk laju log1p sangat zero-inflated. V2 tidak dipilih karena menang di semua metrik; V2 dipilih karena menjadi konfigurasi paling kuat di antara kandidat yang lolos gate semantik, cross-scaler, dan pipeline-consensus.

V1: 8 fitur → StandardScaler → PCA 6 → K-Means k=6
                  ↓ robustness review
5 feature sets × 2 scalers × PCA/original × 2 algorithms × k=3…8
                  ↓ semantic + stability gates
V2: 8 fitur → StandardScaler → PCA 7 → K-Means k=3

from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.metrics import adjusted_rand_score, normalized_mutual_info_score, silhouette_score
from sklearn.preprocessing import StandardScaler

eda = pd.read_csv(ROOT / "data/analytics/village_eda_features.csv", dtype={"village_id": str})
v1 = pd.read_csv(ROOT / "data/modeling/cluster_assignments.csv", dtype={"village_id": str})
v2 = pd.read_csv(ROOT / "data/modeling/cluster_assignments_v2.csv", dtype={"village_id": str})
v1_config = json.loads((ROOT / "data/modeling/selected_model_config.json").read_text())
v2_config = json.loads((ROOT / "data/modeling/selected_model_config_v2.json").read_text())
experiments = pd.read_csv(ROOT / "reports/task05b/candidate_experiments.csv")
finalists = pd.read_csv(ROOT / "reports/task05b/final_candidate_models.csv")
poly_review = pd.read_csv(ROOT / "reports/task05b/polyclinic_representation_comparison.csv")
feature_ablation = pd.read_csv(ROOT / "reports/task05b/feature_ablation.csv")

pd.Series({
    "unit analisis": eda.village_id.nunique(),
    "seluruh eksperimen": len(experiments),
    "finalis dengan evaluasi mendalam": len(finalists),
    "feature sets": experiments.feature_set.nunique(),
    "rentang k": f"{experiments.k.min()}–{experiments.k.max()}",
    "algoritma": ", ".join(sorted(experiments.algorithm.unique())),
}).to_frame("nilai")

	nilai
unit analisis	365
seluruh eksperimen	504
finalis dengan evaluasi mendalam	25
feature sets	5
rentang k	3–8
algoritma	AGGLOMERATIVE_WARD, KMEANS

1 — Reproduksi V1 dan V2

def reproduce(config, frame):
    values = frame.copy()
    if "has_polyclinic" in config["features"]:
        values["has_polyclinic"] = values["polyclinic_count"].gt(0).astype(int)
    X = values[config["features"]].astype(float)
    scaled = StandardScaler().fit_transform(X)
    n_components = config.get("number_of_pca_components", config.get("pca_components"))
    projected = PCA(n_components=n_components, random_state=config["random_state"]).fit_transform(scaled)
    labels = KMeans(
        n_clusters=config["k"], random_state=config["random_state"], n_init=20
    ).fit_predict(projected)
    return projected, labels


v1_space, v1_reproduced = reproduce(v1_config, eda)
v2_space, v2_reproduced = reproduce(v2_config, eda)

reproduction = pd.DataFrame({
    "model": ["V1", "V2"],
    "k": [v1_config["k"], v2_config["k"]],
    "ARI_vs_frozen_assignment": [
        adjusted_rand_score(v1.cluster_id, v1_reproduced),
        adjusted_rand_score(v2.cluster_id, v2_reproduced),
    ],
    "reproduced_silhouette": [
        silhouette_score(v1_space, v1_reproduced),
        silhouette_score(v2_space, v2_reproduced),
    ],
})
reproduction.round(6)

	model	k	ARI_vs_frozen_assignment	reproduced_silhouette
0	V1	6	1.0	0.222292
1	V2	3	1.0	0.175881

assert reproduction.ARI_vs_frozen_assignment.min() > .999
print("✓ V1 dan V2 dapat direproduksi dari feature layer dengan ARI permutation-invariant ≈ 1.")

✓ V1 dan V2 dapat direproduksi dari feature layer dengan ARI permutation-invariant ≈ 1.

ARI dipakai karena nomor cluster bersifat nominal: cluster 1 pada dua hasil fit tidak harus memiliki nomor yang sama untuk mewakili partisi yang sama.

2 — Apa yang memicu review V1?

v1_candidate_id = "T05B_A_CURRENT_8_KM_STANDARD_PCA6_K6"
v2_candidate_id = "T05B_C_BINARY_POLYCLINIC_8_KM_STANDARD_PCA7_K3"
v1_row = finalists.loc[finalists.model_id.eq(v1_candidate_id)].iloc[0]
v2_row = finalists.loc[finalists.model_id.eq(v2_candidate_id)].iloc[0]

trigger_summary = pd.Series({
    "V1 silhouette": v1_config["metrics"]["silhouette"],
    "V1 bootstrap stability ARI": v1_config["metrics"]["bootstrap_mean_ari"],
    "V1 Standard-vs-Robust ARI": v1_row["scaler_robustness"],
    "Polyclinic zero percentage": eda["polyclinic_count"].eq(0).mean() * 100,
    "Polyclinic IQR": eda["polyclinic_count"].quantile(.75) - eda["polyclinic_count"].quantile(.25),
}).round(3)
trigger_summary.to_frame("nilai")

	nilai
V1 silhouette	0.222
V1 bootstrap stability ARI	0.939
V1 Standard-vs-Robust ARI	0.270
Polyclinic zero percentage	86.849
Polyclinic IQR	0.000

ablation_plot = feature_ablation.sort_values("fragility_score_1_minus_ari")
fig, ax = plt.subplots(figsize=(9, 5))
bars = ax.barh(
    ablation_plot.removed_feature.str.replace("log1p_", "", regex=False),
    ablation_plot.fragility_score_1_minus_ari,
    color="#f59e0b",
)
ax.set(
    title="V1 berubah material ketika satu fitur dihapus",
    xlabel="Fragility = 1 − ARI (lebih tinggi = lebih sensitif)",
    ylabel="Fitur yang dihapus",
    xlim=(0, .55),
)
ax.bar_label(bars, fmt="%.2f", padding=3, fontsize=9)
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()

ablation_plot[["removed_feature", "domain", "ari_vs_v1", "fragility_score_1_minus_ari"]].round(3)

	removed_feature	domain	ari_vs_v1	fragility_score_1_minus_ari
7	cellular_operator_count	CONNECTIVITY	0.805	0.195
6	log1p_bts_per_10000_population	CONNECTIVITY	0.795	0.205
5	log1p_population_density	GEOGRAPHY	0.791	0.209
4	population_share_of_district	DEMOGRAPHY	0.781	0.219
3	log1p_polyclinic_per_10000_population	HEALTHCARE	0.583	0.417
2	sex_ratio	DEMOGRAPHY	0.582	0.418
1	has_banking_facility	ECONOMY	0.575	0.425
0	any_documented_disaster_event_2024	DISASTER	0.504	0.496

V1 memiliki separation dan repeatability yang baik, tetapi repeatability di dalam satu spesifikasi tidak sama dengan robustness terhadap pilihan preprocessing atau feature representation yang sama-sama masuk akal.

3 — Audit representasi poliklinik

poly_table = (
    poly_review.groupby("representation", as_index=False)
    .agg(
        cross_scaler_ARI=("scaler_agreement_ari", "first"),
        mean_silhouette=("silhouette", "mean"),
        min_cluster_size=("min_cluster_size", "min"),
    )
    .sort_values("cross_scaler_ARI", ascending=False)
)
poly_table.round(3)

	representation	cross_scaler_ARI	mean_silhouette	min_cluster_size
4	RAW_COUNT	0.346	0.200	23
2	EXCLUDED	0.309	0.192	29
1	CURRENT_LOG1P_RATE	0.270	0.194	28
0	BINARY_PRESENCE	0.265	0.196	27
3	LOG1P_COUNT	0.250	0.194	25

fig, ax = plt.subplots(figsize=(8, 4.5))
display_poly = poly_table.sort_values("cross_scaler_ARI")
bars = ax.barh(
    display_poly.representation.str.replace("_", " "),
    display_poly.cross_scaler_ARI,
    color=["#2563eb" if value == "BINARY_PRESENCE" else "#94a3b8" for value in display_poly.representation],
)
ax.axvline(.50, color="#dc2626", linestyle="--", linewidth=1.2, label="gate final = 0.50")
ax.set(
    title="Mengganti representasi saja belum menyelesaikan sensitivitas k=6",
    xlabel="ARI antara StandardScaler dan RobustScaler",
    ylabel="",
    xlim=(0, .55),
)
ax.bar_label(bars, fmt="%.3f", padding=3, fontsize=9)
ax.legend(loc="lower right")
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()

Keputusan semantik adalah mengganti log1p_polyclinic_per_10000_population dengan has_polyclinic. Presence tetap hanya proxy ketersediaan—bukan kapasitas, kualitas layanan, atau akses fisik. Grafik juga menunjukkan bahwa perubahan representasi saja belum cukup; k dan keseluruhan pipeline harus dievaluasi ulang.

4 — Ruang eksperimen dan hard gate

experiment_design = (
    experiments.groupby(["feature_set", "scaler", "algorithm"])
    .agg(configurations=("model_id", "count"), k_min=("k", "min"), k_max=("k", "max"))
    .reset_index()
)
experiment_design

	feature_set	scaler	algorithm	configurations	k_min	k_max
0	A_CURRENT_8	ROBUST	AGGLOMERATIVE_WARD	30	3	8
1	A_CURRENT_8	ROBUST	KMEANS	30	3	8
2	A_CURRENT_8	STANDARD	AGGLOMERATIVE_WARD	30	3	8
3	A_CURRENT_8	STANDARD	KMEANS	30	3	8
4	B_NO_POLYCLINIC_7	ROBUST	AGGLOMERATIVE_WARD	24	3	8
5	B_NO_POLYCLINIC_7	ROBUST	KMEANS	24	3	8
6	B_NO_POLYCLINIC_7	STANDARD	AGGLOMERATIVE_WARD	24	3	8
7	B_NO_POLYCLINIC_7	STANDARD	KMEANS	24	3	8
8	C_BINARY_POLYCLINIC_8	ROBUST	AGGLOMERATIVE_WARD	30	3	8
9	C_BINARY_POLYCLINIC_8	ROBUST	KMEANS	30	3	8
10	C_BINARY_POLYCLINIC_8	STANDARD	AGGLOMERATIVE_WARD	30	3	8
11	C_BINARY_POLYCLINIC_8	STANDARD	KMEANS	30	3	8
12	D_REDUCED_NO_POLY_NO_SEX_6	ROBUST	AGGLOMERATIVE_WARD	18	3	8
13	D_REDUCED_NO_POLY_NO_SEX_6	ROBUST	KMEANS	18	3	8
14	D_REDUCED_NO_POLY_NO_SEX_6	STANDARD	AGGLOMERATIVE_WARD	18	3	8
15	D_REDUCED_NO_POLY_NO_SEX_6	STANDARD	KMEANS	18	3	8
16	E_BINARY_POLY_NO_SEX_7	ROBUST	AGGLOMERATIVE_WARD	24	3	8
17	E_BINARY_POLY_NO_SEX_7	ROBUST	KMEANS	24	3	8
18	E_BINARY_POLY_NO_SEX_7	STANDARD	AGGLOMERATIVE_WARD	24	3	8
19	E_BINARY_POLY_NO_SEX_7	STANDARD	KMEANS	24	3	8

finalists["passes_gate"] = (
    finalists.algorithm.eq("KMEANS")
    & finalists.min_cluster_size.ge(15)
    & finalists.interpretability_rating.isin(["HIGH", "MEDIUM"])
    & finalists.feature_set.ne("A_CURRENT_8")
    & finalists.scaler_robustness.ge(.50)
    & finalists.pipeline_consensus.ge(.60)
)

top_score = finalists.nlargest(1, "robustness_score").iloc[0]
decision_rows = finalists.loc[
    finalists.model_id.isin([v1_candidate_id, v2_candidate_id, top_score.model_id]),
    [
        "model_id", "feature_set", "k", "silhouette", "bootstrap_ARI",
        "scaler_robustness", "feature_ablation_robustness", "pipeline_consensus",
        "parsimony_score", "robustness_score", "passes_gate", "selected_final",
    ],
].copy()
decision_rows["role"] = decision_rows.model_id.map({
    v1_candidate_id: "V1 baseline",
    v2_candidate_id: "V2 selected",
    top_score.model_id: "Highest raw score",
})
decision_rows.set_index("role").round(3)

	model_id	feature_set	k	silhouette	bootstrap_ARI	scaler_robustness	feature_ablation_robustness	pipeline_consensus	parsimony_score	robustness_score	passes_gate	selected_final
role
Highest raw score	T05B_D_REDUCED_NO_POLY_NO_SEX_6_KM_STANDARD_PCA5_K6	D_REDUCED_NO_POLY_NO_SEX_6	6	0.324	0.990	0.260	0.728	0.551	0.725	0.720	False	False
V2 selected	T05B_C_BINARY_POLYCLINIC_8_KM_STANDARD_PCA7_K3	C_BINARY_POLYCLINIC_8	3	0.176	0.828	0.549	0.622	0.673	0.700	0.677	True	True
V1 baseline	T05B_A_CURRENT_8_KM_STANDARD_PCA6_K6	A_CURRENT_8	6	0.222	0.938	0.270	0.677	0.539	0.425	0.663	False	False

fig, ax = plt.subplots(figsize=(9, 6))
failed = finalists.loc[~finalists.passes_gate]
passed = finalists.loc[finalists.passes_gate]
ax.scatter(failed.silhouette, failed.scaler_robustness, color="#cbd5e1", s=55, label="Tidak lolos gate")
ax.scatter(passed.silhouette, passed.scaler_robustness, color="#2563eb", s=70, label="Lolos gate")
ax.scatter(v2_row.silhouette, v2_row.scaler_robustness, color="#dc2626", marker="*", s=230, label="V2 terpilih")
ax.axhline(.50, color="#475569", linestyle="--", linewidth=1, label="cross-scaler gate")
ax.annotate("V1", (v1_row.silhouette, v1_row.scaler_robustness), xytext=(8, -13), textcoords="offset points")
ax.annotate("V2", (v2_row.silhouette, v2_row.scaler_robustness), xytext=(8, 7), textcoords="offset points")
ax.annotate("raw score tertinggi (gagal gate)", (top_score.silhouette, top_score.scaler_robustness), xytext=(8, 7), textcoords="offset points")
ax.set(
    title="Silhouette yang lebih tinggi tidak otomatis lebih defensible",
    xlabel="Silhouette (lebih tinggi = separation internal lebih baik)",
    ylabel="Cross-scaler ARI (lebih tinggi = lebih robust)",
)
ax.legend(loc="best")
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()

Kandidat raw-score tertinggi menggunakan enam fitur dan k=6, tetapi cross-scaler ARI serta pipeline consensus-nya berada di bawah gate. V2 adalah kandidat ber-score tertinggi setelah syarat berikut diterapkan:

• K-Means dan minimum cluster size ≥15;
• interpretability HIGH atau MEDIUM;
• tidak mempertahankan representasi poliklinik V1;
• cross-scaler ARI ≥0,50;
• pipeline consensus ≥0,60.

Dengan kata lain, pemilihan final adalah constrained model selection, bukan leaderboard satu metrik.

5 — Trade-off V1 versus V2

comparison = pd.DataFrame({
    "metric": [
        "Jumlah cluster", "Jumlah fitur", "Komponen PCA", "Silhouette ↑",
        "Davies–Bouldin ↓", "Calinski–Harabasz ↑", "Bootstrap stability ARI ↑",
        "Cross-scaler ARI ↑", "Feature-ablation robustness ↑", "Pipeline consensus ↑",
        "Cluster terkecil", "Cluster terbesar",
    ],
    "V1": [
        v1_config["k"], len(v1_config["features"]), v1_config["number_of_pca_components"],
        v1_config["metrics"]["silhouette"], v1_config["metrics"]["davies_bouldin"],
        v1_config["metrics"]["calinski_harabasz"], v1_config["metrics"]["bootstrap_mean_ari"],
        v1_row.scaler_robustness, v1_row.feature_ablation_robustness, v1_row.pipeline_consensus,
        v1.cluster_id.value_counts().min(), v1.cluster_id.value_counts().max(),
    ],
    "V2": [
        v2_config["k"], len(v2_config["features"]), v2_config["pca_components"],
        v2_config["metrics"]["silhouette"], v2_config["metrics"]["davies_bouldin"],
        v2_config["metrics"]["calinski_harabasz"], v2_config["metrics"]["bootstrap_stability_ari"],
        v2_config["metrics"]["scaler_robustness"], v2_config["metrics"]["feature_ablation_robustness"],
        v2_config["metrics"]["pipeline_consensus"], v2.cluster_id.value_counts().min(),
        v2.cluster_id.value_counts().max(),
    ],
})
comparison.set_index("metric").round(3)

	V1	V2
metric
Jumlah cluster	6.000	3.000
Jumlah fitur	8.000	8.000
Komponen PCA	6.000	7.000
Silhouette ↑	0.222	0.176
Davies–Bouldin ↓	1.374	1.857
Calinski–Harabasz ↑	86.124	83.430
Bootstrap stability ARI ↑	0.939	0.839
Cross-scaler ARI ↑	0.270	0.549
Feature-ablation robustness ↑	0.677	0.622
Pipeline consensus ↑	0.539	0.673
Cluster terkecil	37.000	60.000
Cluster terbesar	100.000	163.000

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5), sharex=True)
v1_sizes = v1.cluster_id.value_counts().sort_index()
v2_sizes = v2.cluster_label.value_counts().reindex(["Typology A", "Typology B", "Typology C"])

bars1 = axes[0].barh([f"Cluster {i}" for i in v1_sizes.index], v1_sizes.values, color="#94a3b8")
axes[0].bar_label(bars1, padding=3)
axes[0].set(title="V1 — enam cluster", xlabel="Jumlah desa/kelurahan", xlim=(0, 180))

bars2 = axes[1].barh(v2_sizes.index, v2_sizes.values, color=[COLORS["A"], COLORS["B"], COLORS["C"]])
axes[1].bar_label(bars2, padding=3)
axes[1].set(title="V2 — tiga tipologi", xlabel="Jumlah desa/kelurahan", xlim=(0, 180))

for axis in axes:
    axis.spines[["top", "right"]].set_visible(False)
plt.tight_layout()

Trade-off-nya eksplisit: V2 menerima separation internal dan bootstrap stability yang lebih rendah daripada V1, sebagai imbalan atas robustness lintas-scaler yang jauh lebih baik, semantics fitur yang lebih aman, pipeline consensus yang lebih kuat, dan struktur yang lebih mudah dikomunikasikan.

6 — Bagaimana assignment berubah?

merged = v1[["village_id", "cluster_id"]].merge(
    v2[["village_id", "cluster_label"]], on="village_id", validate="one_to_one"
)
transition = pd.crosstab(merged.cluster_id, merged.cluster_label)

fig, ax = plt.subplots(figsize=(8, 5))
sns.heatmap(transition, annot=True, fmt="d", cmap="Blues", linewidths=.5, ax=ax)
ax.set(
    title="V2 bukan sekadar mengganti nama cluster V1",
    xlabel="V2",
    ylabel="V1 cluster",
)
plt.tight_layout()

transition

cluster_label	Typology A	Typology B	Typology C
cluster_id
1	19	0	21
2	0	22	78
3	2	10	60
4	0	35	2
5	39	1	0
6	0	74	2

pd.Series({
    "ARI V1–V2": adjusted_rand_score(merged.cluster_id, merged.cluster_label),
    "NMI V1–V2": normalized_mutual_info_score(merged.cluster_id, merged.cluster_label),
    "changed after maximum-overlap matching": v2_config["comparison_to_v1"]["changed_after_cluster_matching"],
}).round(3).to_frame("nilai")

	nilai
ARI V1–V2	0.337
NMI V1–V2	0.467
changed after maximum-overlap matching	174.000

Heatmap menunjukkan bahwa V2 bukan penggabungan satu-ke-satu dari enam cluster lama. Karena jumlah cluster dan representasi fitur berubah, ARI/NMI dipakai untuk membandingkan struktur partisi. Angka “changed after matching” hanya alat diagnostik berbasis maximum overlap; ia tidak berarti 174 desa sebelumnya “salah”.

7 — Kesimpulan keputusan model

Mengapa V2 dipilih:

1. V1 cukup stabil dalam pengulangan yang sama, tetapi rentan terhadap preprocessing yang wajar.
2. Fitur poliklinik V1 memiliki IQR nol dan 86,85% zero, sehingga laju per kapita memberi presisi yang kurang tepat secara semantik.
3. Sebanyak 504 konfigurasi diuji; 25 finalis menerima evaluasi stability dan ablation tambahan.
4. Hard gate mencegah kandidat dengan silhouette atau raw score tinggi menang ketika cross-scaler robustness/pipeline consensus lemah.
5. V2 k=3 adalah kandidat terkuat yang lolos seluruh gate dan menghasilkan tipologi yang lebih parsimonious.

Yang tidak boleh disimpulkan: V2 bukan bukti adanya tepat tiga “jenis desa alami”, bukan ranking pembangunan, dan bukan klasifikasi resmi. Assignment tetap memiliki 54 unit dengan ambiguitas tinggi; ketidakpastian tersebut dipublikasikan, bukan disembunyikan.