LATAR BELAKANG

Stunting merupakan kondisi gagal tumbuh akibat kekurangan gizi, infeksi berulang, dan kondisi lingkungan yang kurang mendukung. Kondisi ini merupakan ancaman serius bagi kualitas sumber daya manusia di Indonesia. Berdasarkan pendataan keluarga tahun 2021, Provinsi Jawa Tengah tercatat memiliki persentase keluarga berisiko stunting yang cukup tinggi, yakni sebesar 50,89%. Sebagai pusat pemerintahan dan ekonomi di Jawa Tengah, Kota Semarang memerlukan strategi pengentasan stunting yang presisi. Penanganan yang optimal membutuhkan pemetaan karakteristik wilayah hingga tingkat kecamatan agar intervensi kebijakan yang dirancang oleh pemerintah daerah dan instansi terkait dapat dialokasikan secara efektif dan tepat sasaran.

Sesuai dengan pedoman Kementerian Kesehatan, intervensi stunting harus dilakukan secara konvergen mencakup intervensi spesifik dan sensitif. Oleh karena itu, pengelompokan wilayah dalam studi ini didasarkan pada empat determinan krusial yang mengawal siklus 1.000 Hari Pertama Kehidupan, yaitu:

  • Morbiditas Maternal (Komplikasi Kehamilan): Komplikasi medis selama masa kehamilan berisiko menghambat fungsi plasenta dan memicu Intrauterine Growth Restriction (IUGR), yang berdampak langsung pada kelahiran bayi dengan panjang dan berat badan di bawah standar.

  • Asupan Gizi Bayi (Cakupan ASI Eksklusif): Ketiadaan praktik ASI eksklusif mengurangi asupan nutrisi dan antibodi alami bayi, sehingga meningkatkan probabilitas terjadinya penyakit infeksi yang menghambat laju pertumbuhan.

  • Kesehatan Lingkungan (Akses Sanitasi Layak): Sebagai faktor tak langsung, praktik sanitasi yang tidak memenuhi syarat (seperti buang air besar sembarangan) meningkatkan risiko infeksi saluran pencernaan kronis, yang merusak kemampuan usus balita dalam menyerap gizi makanan.

  • Beban Malnutrisi Multidimensi (Indeks Komposit Status Gizi): Studi ini mengaplikasikan konsep Composite Index of Anthropometric Failure (CIAF) dengan merata-ratakan indikator balita pendek (TB/U), berat badan kurang (BB/U), dan gizi kurang (BB/TB). Indeks komposit ini memberikan estimasi kerentanan wilayah terhadap malnutrisi.

Mengingat kompleksitas dan keterkaitan antar determinan pemicu stunting, analisis deskriptif univariat tidak memadai untuk memetakan kondisi riil di lapangan. Oleh karena itu, diperlukan pendekatan multivariat berbasis Machine Learning guna membentuk profil kerentanan wilayah secara komprehensif. Penelitian ini mengaplikasikan algoritma klasterisasi Fuzzy Possibilistic C-Means (FPCM) karena keandalannya dalam memproses data kewilayahan.

Algoritma ini dipilih berdasarkan dua keunggulan utama. Pertama, FPCM mampu menangani karakteristik data antar kecamatan yang sering kali memiliki nilai keanggotaan ganda atau batas kondisi yang saling tumpang tindih (overlapping). Kedua, kombinasi fungsi keanggotaan dan kekhasan pada FPCM menjadikannya tangguh (robust) terhadap keberadaan data pencilan. Kemampuan ini sangat krusial mengingat data indikator kesehatan masyarakat sering kali memiliki nilai ekstrem yang dapat mendistorsi pembentukan pusat klaster pada metode tradisional. Melalui pendekatan algoritma ini, hasil pemetaan kecamatan di Kota Semarang diharapkan menjadi lebih presisi dan objektif sebagai landasan penentuan prioritas intervensi kebijakan.

DATA dan VARIABEL

Data yang digunakan dalam penelitian ini merupakan data sekunder yang bersumber dari dokumen resmi Profil Kesehatan Kota Semarang Tahun 2024 yang diterbitkan oleh Dinas Kesehatan Kota Semarang. Penggunaan publikasi resmi dari instansi kesehatan ini memastikan bahwa pemodelan klasterisasi didasarkan pada rekam medis dan hasil penapisan (screening) lapangan yang valid, akurat, dan dapat dipertanggungjawabkan.

Secara struktur, dataset mentah pada laporan tersebut direkapitulasi pada tingkat fasilitas pelayanan kesehatan, yang mencakup 39 Pusat Kesehatan Masyarakat (Puskesmas) di wilayah Kota Semarang. Untuk memenuhi kebutuhan analisis spasial kewilayahan, data dari ke-39 Puskesmas tersebut selanjutnya melalui proses data preprocessing berupa agregasi (penjumlahan) dan standardisasi persentase ke dalam 16 wilayah administrasi tingkat Kecamatan.

Variabel yang digunakan adalah:

Nama Variabel Keterangan 
Asi_Ekslusif
 Persentase jumlah bayi berusia < 6 bulan yang diberi ASI Eksklusif 
Akses_Sanitasi_Layak
 Persentase jumlah kepala keluarga dengan akses terhadap fasilitas sanitasi layak 
Komplikasi_Kehamilan
 Persentase jumlah ibu yang mengalami kasus komplikasi kehamilan 
Indeks_Status_Gizi
 Rata-rata persentase balita pendek, berat badan kurang, dan gizi kurang sebagai indikator asupan gizi balita

TINJAUAN PUSTAKA

1. Standarisasi Data & Deteksi Outlier

Standarisasi Z-Score dilakukan agar tidak ada variabel yang mendominasi perhitungan jarak:

\[Z_{ij} = \frac{X_{ij} - \bar{X}_j}{s_j}\] Pencilan (outlier) dideteksi menggunakan Mahalanobis Distance yang menghasilkan skor anomali: \[D^2(x) = (x - \mu)^T \Sigma^{-1} (x - \mu)\]

2. Uji Multikolinearitas (VIF)

Untuk mendeteksi redundansi antarvariabel, digunakan Variance Inflation Factor (VIF): \[VIF_j = \frac{1}{1 - R_j^2}\]

3. Penentuan Klaster Optimal

Digunakan Modified Partition Coefficient (MPC) untuk menentukan jumlah klaster optimal: \[PC(c) = \frac{1}{n}\sum_{k=1}^{c}\sum_{i=1}^{n}\mu_{ik}^2\]\[MPC(c) = 1 - \frac{c}{c-1}(1 - PC(c))\]

4. Algoritma FPCM

FPCM adalah pengembangan fuzzy yang robust terhadap outlier. Fungsi objektif yang diminimalkan adalah: \[J_{FPCM} = \sum_{i=1}^{n}\sum_{k=1}^{c}(\mu_{ik}^w + t_{ik}^\eta)||X_i - V_k||^2\]Nilai Pusat Klaster (\(V_k\)), Derajat Keanggotaan (\(\mu_{ik}\)), dan Tipikalitas (\(t_{ik}\)) diperbarui melalui: \[V_k = \frac{\sum_{i=1}^{n}(\mu_{ik}^w+t_{ik}^\eta)X_i}{\sum_{i=1}^{n}(\mu_{ik}^w+t_{ik}^\eta)}\]\[\mu_{ik} = \left[ \sum_{j=1}^{c} \left( \frac{||X_i - V_k||^2}{||X_i - V_j||^2} \right)^{\frac{1}{w-1}} \right]^{-1}\]\[t_{ik} = \left[ 1 + \left( \frac{||X_i-V_k||^2}{\eta_k} \right)^{\frac{1}{\eta-1}} \right]^{-1}\]

5. Validasi & Feature Importance

Setelah klaster optimal terbentuk, signifikansi perbedaan karakteristik antar klaster dievaluasi menggunakan uji statistik non-parametrik Kruskal-Wallis. Uji ini digunakan untuk memastikan apakah setiap variabel determinan secara statistik memiliki perbedaan median yang nyata dalam membedakan kelompok klaster:\[H = \frac{12}{N(N+1)} \sum_{i=1}^{c} \frac{R_i^2}{n_i} - 3(N+1)\] Selanjutnya, untuk mengetahui besaran kontribusi tiap variabel dalam mengelompokkan wilayah (Feature Importance), dievaluasi menggunakan algoritma Random Forest. Tingkat kepentingan variabel diukur melalui metrik Mean Decrease Accuracy (MDA), yang menghitung seberapa besar penurunan akurasi prediksi Out-Of-Bag (OOB) model apabila nilai pada variabel tersebut diacak:\[MDA_j = \text{Akurasi}_{OOB}(\text{asli}) - \text{Akurasi}_{OOB}(\text{Variabel } j \text{ diacak})\]

PEMBAHASAN

1. Bagian ini memuat library, membaca data, dan menyesuaikan jumlah variabel secara dinamis.

# ==========================================
# 1. KONFIGURASI AWAL
# ==========================================
file_path  <- "D://GINSA//KULIAH//SMT 6//KP//Data.xlsx"
sheet_name <- "FIX"
col_id     <- "KECAMATAN"

#Input variabel indikator
vars_cluster <- c(
  "Asi_Ekslusif",
  "Akses_Sanitasi_Layak",
  "Komplikasi_Kehamilan",
  "Indeks_Status_Gizi"
)

# ==========================================
# 2. LOAD LIBRARIES & DATA
# ==========================================
library(tidyverse)
library(openxlsx)
library(ppclust)
library(solitude)
library(sf)
library(geodata)
library(ggspatial)
library(MASS)
library(biotools)
library(caret)
library(FSA)
library(randomForest)

df <- read.xlsx(file_path, sheet = sheet_name)

df_clean <- df %>%
  filter(!is.na(!!sym(col_id))) %>%
  mutate(across(all_of(vars_cluster), ~ as.numeric(str_replace(as.character(.), ",", ".")))) %>%
  dplyr::select(all_of(col_id), all_of(vars_cluster))

df_clean
##           KECAMATAN Asi_Ekslusif Akses_Sanitasi_Layak Komplikasi_Kehamilan
## 1        BANYUMANIK     94.97908             98.87479            9.8857143
## 2         CANDISARI     89.98273             96.20045           12.6984127
## 3      GAJAHMUNGKUR     86.90096             88.91024            0.8658009
## 4         GAYAMSARI     87.53181             99.75289           11.9341564
## 5             GENUK     89.82940             99.47856           14.2099682
## 6        GUNUNGPATI     91.08434             99.92974            9.9032018
## 7             MIJEN     99.87374             99.27540           13.7580794
## 8          NGALIYAN     90.81232             99.62207           19.3188163
## 9        PEDURUNGAN     95.14360             99.98215           10.4959076
## 10   SEMARANG BARAT     87.94712             99.91075            9.3278464
## 11 SEMARANG SELATAN     90.08746             94.40715            3.4205231
## 12  SEMARANG TENGAH     89.12281             80.83185           15.1125402
## 13   SEMARANG TIMUR     90.03497             99.29344           10.5169340
## 14   SEMARANG UTARA     82.14592             96.25350           19.0261497
## 15        TEMBALANG     87.34908             99.99821            5.1948052
## 16             TUGU     93.90863             99.44950           13.8755981
##    Indeks_Status_Gizi
## 1           0.9296893
## 2           1.6187433
## 3           1.1363636
## 4           1.4020662
## 5           0.3783805
## 6           0.6434135
## 7           0.7484735
## 8           1.0870939
## 9           0.7321515
## 10          1.6328031
## 11          2.9101778
## 12          3.4515100
## 13          2.6851852
## 14          3.9880041
## 15          0.5011082
## 16          1.2968300

2. Pada bagian ini dilakukan scaling dengan Z-score

df_subset <- df_clean %>% dplyr::select(all_of(vars_cluster))
X_scale   <- scale(df_subset) # Standarisasi Z-Score
print(head(X_scale))
##      Asi_Ekslusif Akses_Sanitasi_Layak Komplikasi_Kehamilan Indeks_Status_Gizi
## [1,]    1.1110981            0.3564432           -0.2656561         -0.5831120
## [2,]   -0.1068010           -0.1549244            0.2937111          0.0430449
## [3,]   -0.8580063           -1.5489059           -2.0594652         -0.3953030
## [4,]   -0.7042322            0.5243489            0.1417218         -0.1538539
## [5,]   -0.1441770            0.4718936            0.5943173         -1.0840972
## [6,]    0.1617246            0.5581637           -0.2621784         -0.8432565

3. Pada bagian ini dilakukan pemeriksaan nilai VIF dan Deteksi Outlier

# ==========================================
# 3. CEK VIF & DETEKSI OUTLIER
# ==========================================
# Uji Multikolinearitas (VIF)
calc_vif <- function(data_sub, vars) {
  vif_df <- data.frame(Variable = vars, VIF = NA)
  for (i in seq_along(vars)) {
    form <- as.formula(paste(vars[i], "~", paste(vars[-i], collapse = " + ")))
    vif_df$VIF[i] <- 1 / (1 - summary(lm(form, data = data_sub))$r.squared)
  }
  return(vif_df)
}

# Deteksi Outlier (Mahalanobis Distance)
# Menghitung pusat (mean) dan matriks kovarians
pusat_data <- colMeans(X_scale)
matriks_kov <- cov(X_scale)

# Menghitung Jarak Mahalanobis untuk setiap observasi
df_clean$Mahalanobis_Dist <- mahalanobis(X_scale, center = pusat_data, cov = matriks_kov)

# Menghitung p-value (menggunakan distribusi Chi-Square, df = jumlah variabel)
df_clean$p_value_MD <- pchisq(df_clean$Mahalanobis_Dist, df = ncol(X_scale), lower.tail = FALSE)

# Melabeli Outlier (batas signifikansi p-value < 0.05 dianggap outlier)
df_clean <- df_clean %>%
  mutate(Status_Outlier = ifelse(p_value_MD < 0.05, "Outlier", "Normal"))

print(calc_vif(df_subset, vars_cluster))
##               Variable      VIF
## 1         Asi_Ekslusif 1.386061
## 2 Akses_Sanitasi_Layak 1.524760
## 3 Komplikasi_Kehamilan 1.180131
## 4   Indeks_Status_Gizi 2.029479
df_clean %>%
  arrange(desc(Mahalanobis_Dist)) %>%
  dplyr::select(all_of(col_id), Mahalanobis_Dist, p_value_MD, Status_Outlier) %>%
  print()
##           KECAMATAN Mahalanobis_Dist p_value_MD Status_Outlier
## 1   SEMARANG TENGAH       11.0918657 0.02555074        Outlier
## 2    SEMARANG UTARA        8.4704985 0.07578652         Normal
## 3      GAJAHMUNGKUR        7.6873741 0.10372513         Normal
## 4  SEMARANG SELATAN        6.4982187 0.16490266         Normal
## 5             MIJEN        5.6073968 0.23044927         Normal
## 6         TEMBALANG        3.6392010 0.45703043         Normal
## 7          NGALIYAN        3.5186979 0.47504099         Normal
## 8    SEMARANG TIMUR        3.2846795 0.51136331         Normal
## 9             GENUK        3.0256370 0.55354431         Normal
## 10       PEDURUNGAN        1.4502981 0.83540754         Normal
## 11       BANYUMANIK        1.3618946 0.85079013         Normal
## 12   SEMARANG BARAT        1.1660537 0.88365541         Normal
## 13        GAYAMSARI        1.1556905 0.88533860         Normal
## 14             TUGU        1.0456970 0.90278947         Normal
## 15       GUNUNGPATI        0.8209346 0.93561966         Normal
## 16        CANDISARI        0.1758623 0.99635338         Normal

Berdasarkan hasil yang diperoleh, seluruh nilai VIF tiap variabel < 10 sehingga dapat disimpulkan tidak terjadi multikolinieritas. Selanjutnya berdasarkan penghitungan outlier dengan Jarak Mahalanobis, diperoleh hasil bahwa terdapat satu wilayah yang dikategorikan sebagai pencilan. Oleh karena itu, metode FPCM akan digunakan untuk proses analisis.

4. Pada bagian ini akan dicari jumlah klaster optimal dengan MPC

# ==========================================
# 4. PENENTUAN K OPTIMAL (MPC)
# ==========================================
calc_mpc <- function(u) {
  n <- nrow(u); k <- ncol(u)
  1 - (k / (k - 1)) * (1 - (sum(u^2) / n))
}

mpc_df <- data.frame(k = 2:10, mpc = NA)
set.seed(123)
for (i in 1:nrow(mpc_df)) {
  res_temp <- fpcm(X_scale, centers = mpc_df$k[i], m = 2)
  mpc_df$mpc[i] <- calc_mpc(res_temp$u)
}

k_optimal <- mpc_df$k[which.max(mpc_df$mpc)]
print(mpc_df)
##    k           mpc
## 1  2  2.698461e-01
## 2  3  4.183771e-01
## 3  4  3.625035e-01
## 4  5  3.073953e-01
## 5  6  1.279917e-01
## 6  7 -2.220446e-16
## 7  8  0.000000e+00
## 8  9  0.000000e+00
## 9 10  0.000000e+00
cat("\n K Optimal terpilih adalah K =", k_optimal, "\n")
## 
##  K Optimal terpilih adalah K = 3
plot(mpc_df$k, mpc_df$mpc, type = "b", pch = 19, col = "darkblue",
     main = "Analisis MPC untuk K Optimal", xlab = "K", ylab = "MPC")

Berdasarkan nilai MPC tertinggi, diperoleh klaster optimal yang terbentuk adalah 3 klaster dengan nilai MPC sebesar 0,418377.

5. Pada bagian ini dilakukan pemodelan FPCM dengan kluster optimal yang telah terpilih

# ==========================================
# 5. MODELING FPCM DENGAN K OPTIMAL
# ==========================================
set.seed(123)
res_final <- fpcm(X_scale, centers = k_optimal, m = 2)

df_clean <- df_clean %>%
  mutate(
    Cluster_Final = apply(res_final$u, 1, which.max),
    Nama_Klaster  = paste("Klaster", Cluster_Final)
  )

df_clean %>%
  group_by(Nama_Klaster) %>%
  summarise(across(all_of(vars_cluster), ~ round(mean(., na.rm = TRUE), 2))) %>%
  print()
## # A tibble: 3 × 5
##   Nama_Klaster Asi_Ekslusif Akses_Sanitasi_Layak Komplikasi_Kehamilan
##   <chr>               <dbl>                <dbl>                <dbl>
## 1 Klaster 1            93.7                 99.5                13.1 
## 2 Klaster 2            88.6                 96.9                 7.71
## 3 Klaster 3            85.6                 88.5                17.1 
## # ℹ 1 more variable: Indeks_Status_Gizi <dbl>
# Menampilkan Daftar Kecamatan per Klaster
df_clean %>%
  dplyr::select(Nama_Klaster, !!sym(col_id)) %>%
  dplyr::arrange(Nama_Klaster, !!sym(col_id)) %>%
  as.data.frame() %>%
  print(row.names = FALSE)
##  Nama_Klaster        KECAMATAN
##     Klaster 1       BANYUMANIK
##     Klaster 1            GENUK
##     Klaster 1       GUNUNGPATI
##     Klaster 1            MIJEN
##     Klaster 1         NGALIYAN
##     Klaster 1       PEDURUNGAN
##     Klaster 1             TUGU
##     Klaster 2        CANDISARI
##     Klaster 2     GAJAHMUNGKUR
##     Klaster 2        GAYAMSARI
##     Klaster 2   SEMARANG BARAT
##     Klaster 2 SEMARANG SELATAN
##     Klaster 2   SEMARANG TIMUR
##     Klaster 2        TEMBALANG
##     Klaster 3  SEMARANG TENGAH
##     Klaster 3   SEMARANG UTARA

Berdasarkan pemodelan, telah diperoleh profiling dan anggota masing-masing klaster dengan deskripsi sebagai berikut:

- Klaster 3: Zona Kritis (Prioritas Intervensi Utama)

Anggota Wilayah: Semarang Tengah, Semarang Utara.

Karakteristik: Merupakan wilayah dengan tingkat kerentanan stunting terburuk. Capaian ASI Eksklusif tercatat paling rendah (85.6%) diiringi dengan akses sanitasi layak yang paling minimum (88.5%). Buruknya asupan bayi dan kesehatan lingkungan ini semakin diperparah oleh tingginya komplikasi kehamilan (17.1%), yang secara akumulatif mendorong Indeks Status Gizi balita (beban malnutrisi) memburuk ke titik tertinggi di angka 3.72%.

- Klaster 2: Zona Transisi (Waspada)

Anggota Wilayah: Candisari, Gajahmungkur, Gayamsari, Semarang Barat, Semarang Selatan, Semarang Timur, Tembalang.

Karakteristik: Merupakan wilayah penyangga dengan tingkat kerentanan menengah. Capaian ASI Eksklusif (88.6%) dan akses sanitasi (96.9%) berada di level moderat. Tingkat komplikasi kehamilan di wilayah ini tercatat paling rendah (7.71%), sehingga berhasil menahan lonjakan Indeks Status Gizi balita di angka pertengahan (1.70%).

- Klaster 1: Zona Ketahanan Kuat (Aman)

Anggota Wilayah: Banyumanik, Genuk, Gunungpati, Mijen, Ngaliyan, Pedurungan, Tugu.

Karakteristik: Merupakan wilayah percontohan dengan pertahanan kesehatan lingkungan dan gizi terbaik. Meskipun angka komplikasi kehamilan berada di taraf moderat (13.1%), capaian ASI Eksklusif yang sangat tinggi (93.7%) dan sanitasi lingkungan yang nyaris sempurna (99.5%) berhasil meredam risiko kelainan tumbuh kembang. Hal ini dibuktikan dengan Indeks Status Gizi yang menempati posisi terbaik (0.83%).

6. Pada bagian ini dilakukan validasi model dan feature importance

for (v in vars_cluster) {
  cat("\n Uji Kruskal-Wallis:", v, "\n")
  form <- as.formula(paste(v, "~ as.factor(Cluster_Final)"))
  kw <- kruskal.test(form, data = df_clean)
  print(kw)
  if (kw$p.value < 0.05) print(dunnTest(form, data = df_clean, method = "bh"))
}
## 
##  Uji Kruskal-Wallis: Asi_Ekslusif 
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Asi_Ekslusif by as.factor(Cluster_Final)
## Kruskal-Wallis chi-squared = 9.4821, df = 2, p-value = 0.008729
##   Comparison        Z    P.unadj      P.adj
## 1      1 - 2 2.638396 0.00832992 0.02498976
## 2      1 - 3 2.376428 0.01748118 0.02622178
## 3      2 - 3 0.617497 0.53690697 0.53690697
## 
##  Uji Kruskal-Wallis: Akses_Sanitasi_Layak 
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Akses_Sanitasi_Layak by as.factor(Cluster_Final)
## Kruskal-Wallis chi-squared = 3.6681, df = 2, p-value = 0.1598
## 
## 
##  Uji Kruskal-Wallis: Komplikasi_Kehamilan 
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Komplikasi_Kehamilan by as.factor(Cluster_Final)
## Kruskal-Wallis chi-squared = 7.062, df = 2, p-value = 0.02928
##   Comparison         Z    P.unadj      P.adj
## 1      1 - 2  1.852491 0.06395536 0.09593304
## 2      1 - 3 -1.178858 0.23845478 0.23845478
## 3      2 - 3 -2.413852 0.01578488 0.04735464
## 
##  Uji Kruskal-Wallis: Indeks_Status_Gizi 
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Indeks_Status_Gizi by as.factor(Cluster_Final)
## Kruskal-Wallis chi-squared = 8.8015, df = 2, p-value = 0.01227
##   Comparison         Z     P.unadj      P.adj
## 1      1 - 2 -1.964763 0.049441661 0.07416249
## 2      1 - 3 -2.750668 0.005947382 0.01784215
## 3      2 - 3 -1.440826 0.149633765 0.14963376
# Validasi Akurasi (LDA LOOCV)
form_lda <- as.formula(paste("Cluster_Final ~", paste(vars_cluster, collapse = "+")))
lda_cv <- lda(form_lda, data = df_clean, CV = TRUE)
akurasi <- sum(diag(table(df_clean$Cluster_Final, lda_cv$class))) / nrow(df_clean)
cat("\nAkurasi Validasi LDA (Hit Ratio):", round(akurasi * 100, 2), "%\n")
## 
## Akurasi Validasi LDA (Hit Ratio): 87.5 %
# Feature Importance (Random Forest)
set.seed(123)
form_rf <- as.formula(paste("as.factor(Cluster_Final) ~", paste(vars_cluster, collapse = "+")))
rf_mod <- randomForest(form_rf, data = df_clean, importance = TRUE)

rf_imp <- as.data.frame(importance(rf_mod)) %>% rownames_to_column("Variable")
ggplot(rf_imp, aes(x = reorder(Variable, MeanDecreaseAccuracy), y = MeanDecreaseAccuracy)) +
  geom_bar(stat = "identity", fill = "#e67e22") +
  coord_flip() + theme_minimal() +
  labs(title = "Feature Importance dalam Pembentukan Klaster", x = "Variabel")

Berdasarkan pengujian Kruskal-Wallis dan evaluasi daya pisah model (Linear Discriminant Analysis), pembentukan klaster kerentanan stunting terbukti valid secara statistik:

Signifikansi Variabel Pembeda (Kruskal-Wallis & Dunn’s Test): Tiga dari empat determinan memiliki p-value < 0.05, yang membuktikan bahwa variabel tersebut secara nyata mampu membedakan karakteristik wilayah:

  • ASI Eksklusif (p = 0.008): Secara signifikan membedakan wilayah dengan ketahanan kuat (Klaster 1) dari wilayah transisi (Klaster 2) maupun kritis (Klaster 3).

  • Indeks Status Gizi (p = 0.012): Menjadi pemisah absolut antara dua kutub wilayah ekstrem, terbukti signifikan membedakan Klaster 1 (kondisi terbaik) dan Klaster 3 (kondisi terburuk).

  • Komplikasi Kehamilan (p = 0.029): Berperan sebagai faktor kunci yang membedakan kerentanan wilayah transisi (Klaster 2) agar tidak terperosok menjadi zona kritis (Klaster 3).

  • Akses Sanitasi Layak (p = 0.159): Menunjukkan hasil yang tidak signifikan secara statistik karena capaian antar wilayah yang sudah relatif tinggi dan merata. Namun, indikator ini tetap dipertahankan dalam model karena urgensi klinisnya sebagai prasyarat dasar intervensi gizi sensitif.

Akurasi Klasifikasi (Validasi LDA): Pengujian Leave-One-Out Cross-Validation (LOOCV) menghasilkan nilai hit ratio sebesar 87.5%. Angka ini mengonfirmasi bahwa algoritma pengelompokan memiliki kinerja prediktif yang baik.

Evaluasi tingkat kepentingan variabel (Feature Importance) menggunakan algoritma Random Forest menunjukkan hasil yang selaras dengan uji Kruskal-Wallis.

PETA KLASTER

# ==========================================
# 6. PEMETAAN GEOSPASIAL (3 KLASTER)
# ==========================================
# Pastikan library ter-load
library(sf)
library(geodata)
library(ggspatial)

# 1. Download & Filter Peta Kota Semarang
idn <- gadm(country = "IDN", level = 3, path = tempdir())
peta_smg <- st_as_sf(idn) %>% filter(NAME_2 == "Kota Semarang")

# 2. Sinkronisasi Nama Kecamatan antara Excel dan Peta
df_map <- df_clean %>%
  mutate(
    Kec_Join = toupper(str_trim(str_squish(!!sym(col_id)))),
    Kec_Join = ifelse(Kec_Join == "GAJAH MUNGKUR", "GAJAHMUNGKUR", Kec_Join)
  )

peta_plot <- peta_smg %>%
  mutate(Kec_Map = toupper(str_trim(str_squish(NAME_3)))) %>%
  inner_join(df_map, by = c("Kec_Map" = "Kec_Join"))

# 3. Plotting Peta
ggplot(data = peta_plot) +
  geom_sf(aes(fill = as.factor(Cluster_Final)), color = "white", size = 0.5) +
  # Atur 3 Warna: 1=Kuning, 2=Merah, 3=Biru
  scale_fill_manual(
    values = c("1" = "green", "2" = "yellow", "3" = "red"),
    name = "Klaster"
  ) +
  geom_sf_text(aes(label = NAME_3), size = 2.5, fontface = "bold", color = "black") +
  # Menambahkan Arah Mata Angin (North Arrow)
  annotation_north_arrow(
    location = "tl", 
    pad_x = unit(0.2, "in"), pad_y = unit(0.2, "in"),
    style = north_arrow_orienteering(text_size = 10)
  ) +
  # Menambahkan Skala Bar (10 km)
  annotation_scale(location = "bl", width_hint = 0.3) +
  # Keterangan Judul
  labs(
    title = "Peta Klaster Kerentanan Stunting Kota Semarang",
    subtitle = "Analisis Fuzzy Possibilistic C-Means (FPCM)"
  ) +
  # Membersihkan background
  theme_void() + 
  theme(
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 12, margin = ggplot2::margin(b = 10)),
    legend.position = "right"
  )
Peta Klaster Kota Semarang
Peta Klaster Kota Semarang