# =========================================================================
# PEMODELAN TINGKAT KEMISKINAN KABUPATEN/KOTA DI INDONESIA MENGGUNAKAN PENDEKATAN REGRESI LOGISTIK ORDINAL
# =========================================================================
# =========================
# LOAD PACKAGE
# =========================
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.3.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.3.3
## corrplot 0.95 loaded
library(car)
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.3
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library(brant)
## Warning: package 'brant' was built under R version 4.3.3
library(scales)
# =========================
# IMPORT DATA
# =========================
data <- read_excel(
"C:/Users/ASUS/OneDrive/Documents/ADK K10/Data tubes ADK.xlsx",
sheet = "Sheet1"
)
# buang kolom kosong tambahan (jika ada) dan ambil 12 kolom utama
data <- data[, 1:12]
# =========================
# RENAME KOLOM
# =========================
colnames(data) <- c(
"Provinsi", # Provinsi
"KabKota", # Kabupaten/Kota
"Kemiskinan", # Persentase Penduduk Miskin (P0) - dasar pembentukan Y
"Sekolah", # Rata-rata Lama Sekolah Penduduk 15+ (Tahun)
"Pengeluaran", # Pengeluaran per Kapita Disesuaikan (Ribu Rupiah/Orang/Tahun)
"IPM", # Indeks Pembangunan Manusia -> DIKELUARKAN dari model (VIF tinggi)
"UHH", # Umur Harapan Hidup (Tahun)
"Sanitasi", # % rumah tangga akses sanitasi layak
"AirMinum", # % rumah tangga akses air minum layak
"Pengangguran", # Tingkat Pengangguran Terbuka
"TPAK", # Tingkat Partisipasi Angkatan Kerja
"PDRB" # PDRB atas Dasar Harga Konstan menurut Pengeluaran (Rupiah)
)
# =========================
# CLEANING DATA
# =========================
# beberapa kolom numerik terbaca sebagai karakter dengan format koma (mis. "18,98") karena format desimal Indonesia, sehingga perlu dikonversi ke numerik
kolom_numerik <- c("Kemiskinan", "Sekolah", "Pengeluaran", "IPM", "UHH",
"Sanitasi", "AirMinum", "Pengangguran", "TPAK", "PDRB")
data <- data %>%
mutate(across(all_of(kolom_numerik),
~ as.numeric(gsub(",", ".", as.character(.)))))
# cek missing value
cat("Jumlah missing value per kolom:\n")
## Jumlah missing value per kolom:
print(colSums(is.na(data)))
## Provinsi KabKota Kemiskinan Sekolah Pengeluaran IPM
## 0 0 0 0 0 0
## UHH Sanitasi AirMinum Pengangguran TPAK PDRB
## 0 0 0 0 0 0
# buang baris dengan missing value (jika ada)
data <- data %>% drop_na(all_of(kolom_numerik))
cat("\nJumlah baris data setelah cleaning:", nrow(data), "\n")
##
## Jumlah baris data setelah cleaning: 514
# =========================
# KLASIFIKASI VARIABEL Y (KEMISKINAN) MENJADI 3 KATEGORI
# =========================
# Y dibentuk dari variabel Kemiskinan (persentase penduduk miskin),
# sehingga Kemiskinan asli TIDAK diikutkan lagi sebagai variabel X.
# Pembagian kategori menggunakan TERSIL (33% & 66%) agar jumlah kab/kota di tiap kategori relatif seimbang.
batas <- quantile(data$Kemiskinan, probs = c(0, 1/3, 2/3, 1))
cat("\nBatas kategori (tersil):\n")
##
## Batas kategori (tersil):
print(batas)
## 0% 33.33333% 66.66667% 100%
## 2.38 7.99 13.41 41.66
data$Kategori <- cut(
data$Kemiskinan,
breaks = batas,
include.lowest = TRUE,
labels = c("Rendah", "Sedang", "Tinggi")
)
data$Kategori <- factor(data$Kategori,
levels = c("Rendah", "Sedang", "Tinggi"),
ordered = TRUE)
cat("\nDistribusi kategori kemiskinan:\n")
##
## Distribusi kategori kemiskinan:
print(table(data$Kategori))
##
## Rendah Sedang Tinggi
## 172 171 171
# =========================
# STATISTIKA DESKRIPTIF
# =========================
str(data)
## tibble [514 × 13] (S3: tbl_df/tbl/data.frame)
## $ Provinsi : chr [1:514] "ACEH" "ACEH" "ACEH" "ACEH" ...
## $ KabKota : chr [1:514] "Simeulue" "Aceh Singkil" "Aceh Selatan" "Aceh Tenggara" ...
## $ Kemiskinan : num [1:514] 19 20.4 13.2 13.4 14.4 ...
## $ Sekolah : num [1:514] 9.48 8.68 8.88 9.67 8.21 ...
## $ Pengeluaran : num [1:514] 7148 8776 8180 8030 8577 ...
## $ IPM : num [1:514] 66.4 69.2 67.4 69.4 67.8 ...
## $ UHH : num [1:514] 65.3 67.4 64.4 68.2 68.7 ...
## $ Sanitasi : num [1:514] 71.6 69.6 62.5 62.7 66.8 ...
## $ AirMinum : num [1:514] 87.5 78.6 79.7 86.7 83.2 ...
## $ Pengangguran: num [1:514] 5.71 8.36 6.46 6.43 7.13 2.61 7.09 7.7 7.28 4.32 ...
## $ TPAK : num [1:514] 71.2 62.9 60.9 69.6 59.5 ...
## $ PDRB : num [1:514] 1648096 1780419 4345784 3487157 8433526 ...
## $ Kategori : Ord.factor w/ 3 levels "Rendah"<"Sedang"<..: 3 3 2 2 3 3 3 3 3 2 ...
summary(data[, c(kolom_numerik, "Kategori")])
## Kemiskinan Sekolah Pengeluaran IPM
## Min. : 2.38 Min. : 1.420 Min. : 3976 Min. :32.84
## 1st Qu.: 7.15 1st Qu.: 7.510 1st Qu.: 8574 1st Qu.:66.64
## Median :10.46 Median : 8.305 Median :10196 Median :69.61
## Mean :12.27 Mean : 8.437 Mean :10325 Mean :69.93
## 3rd Qu.:14.89 3rd Qu.: 9.338 3rd Qu.:11719 3rd Qu.:73.11
## Max. :41.66 Max. :12.830 Max. :23888 Max. :87.18
## UHH Sanitasi AirMinum Pengangguran
## Min. :55.43 Min. : 0.00 Min. : 0.00 Min. : 0.000
## 1st Qu.:67.39 1st Qu.:70.22 1st Qu.: 79.04 1st Qu.: 3.180
## Median :69.97 Median :81.80 Median : 89.80 Median : 4.565
## Mean :69.66 Mean :77.20 Mean : 85.14 Mean : 5.059
## 3rd Qu.:72.04 3rd Qu.:89.88 3rd Qu.: 96.40 3rd Qu.: 6.530
## Max. :77.73 Max. :99.97 Max. :100.00 Max. :13.370
## TPAK PDRB Kategori
## Min. :56.39 Min. : 147485 Rendah:172
## 1st Qu.:65.07 1st Qu.: 3654292 Sedang:171
## Median :68.95 Median : 8814926 Tinggi:171
## Mean :69.46 Mean : 21964077
## 3rd Qu.:72.34 3rd Qu.: 19735101
## Max. :97.93 Max. :460081046
# =========================
# EKSPLORASI DATA (EDA)
# =========================
# variabel X yang dipakai dalam model (semua kecuali Kemiskinan & IPM)
var_x <- c("Sekolah", "Pengeluaran", "UHH", "Sanitasi",
"AirMinum", "Pengangguran", "TPAK", "PDRB")
# 1. distribusi kategori kemiskinan
ggplot(data, aes(x = Kategori, fill = Kategori)) +
geom_bar() +
geom_text(stat = "count", aes(label = after_stat(count)), vjust = -0.5) +
scale_fill_brewer(palette = "YlOrRd") +
labs(title = "Distribusi Kategori Tingkat Kemiskinan Kab/Kota",
x = "Kategori", y = "Jumlah Kab/Kota") +
theme_minimal() +
theme(legend.position = "none")
ggsave("01_distribusi_kategori.png", width = 6, height = 4)
# 2. boxplot tiap variabel X terhadap kategori Y
for (v in var_x) {
p <- ggplot(data, aes(x = Kategori, y = .data[[v]], fill = Kategori)) +
geom_boxplot() +
scale_fill_brewer(palette = "YlOrRd") +
labs(title = paste("Sebaran", v, "Menurut Kategori Kemiskinan"),
x = "Kategori", y = v) +
theme_minimal() +
theme(legend.position = "none")
ggsave(paste0("02_boxplot_", v, ".png"), p, width = 6, height = 4)
}
# 3. matriks korelasi antar variabel X (cek hubungan & gambaran awal multikolinearitas)
png("03_korelasi_antar_X.png", width = 800, height = 800)
corrplot(cor(data[var_x]),
method = "color", type = "upper",
addCoef.col = "black", tl.col = "black", tl.srt = 45,
col = colorRampPalette(c("firebrick", "white", "steelblue"))(200))
dev.off()
## png
## 2
# =========================
# UJI MULTIKOLINEARITAS (VIF)
# =========================
# dicek dulu dengan menyertakan IPM untuk membuktikan IPM punya VIF tinggi
model_cek_ipm <- lm(
Kemiskinan ~ Sekolah + Pengeluaran + IPM + UHH + Sanitasi +
AirMinum + Pengangguran + TPAK + PDRB,
data = data
)
cat("\n=== VIF dengan IPM disertakan ===\n")
##
## === VIF dengan IPM disertakan ===
print(vif(model_cek_ipm))
## Sekolah Pengeluaran IPM UHH Sanitasi AirMinum
## 9.488522 8.580938 36.070465 4.491066 2.277657 1.637052
## Pengangguran TPAK PDRB
## 2.159531 1.951017 1.577570
# VIF setelah IPM dikeluarkan
model_cek_noipm <- lm(
Kemiskinan ~ Sekolah + Pengeluaran + UHH + Sanitasi +
AirMinum + Pengangguran + TPAK + PDRB,
data = data
)
cat("\n=== VIF tanpa IPM (model final) ===\n")
##
## === VIF tanpa IPM (model final) ===
print(vif(model_cek_noipm))
## Sekolah Pengeluaran UHH Sanitasi AirMinum Pengangguran
## 2.365437 2.872825 1.652151 2.206627 1.605275 2.132952
## TPAK PDRB
## 1.808182 1.443317
# VIF IPM jauh di atas 10 (indikasi multikolinearitas tinggi), sedangkan tanpa IPM seluruh VIF variabel turun jauh di bawah 10. Sehingga IPM dikeluarkan dari variabel X.
# =========================
# PEMODELAN REGRESI LOGISTIK ORDINAL
# =========================
# Variabel X memiliki skala yang sangat berbeda (mis. PDRB hingga ratusan juta, sedangkan Sekolah hanya puluhan), sehingga seluruh prediktor distandarisasi (z-score) agar estimasi model lebih stabil.
# Konsekuensinya: odds ratio nanti diinterpretasikan "per kenaikan 1 SD".
data_model <- data
data_model[var_x] <- scale(data_model[var_x])
model <- polr(
Kategori ~ Sekolah + Pengeluaran + UHH + Sanitasi +
AirMinum + Pengangguran + TPAK + PDRB,
data = data_model,
Hess = TRUE,
method = "logistic" # model logit kumulatif (proportional odds)
)
cat("\n=== Ringkasan Model Regresi Logistik Ordinal ===\n")
##
## === Ringkasan Model Regresi Logistik Ordinal ===
print(summary(model))
## Call:
## polr(formula = Kategori ~ Sekolah + Pengeluaran + UHH + Sanitasi +
## AirMinum + Pengangguran + TPAK + PDRB, data = data_model,
## Hess = TRUE, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## Sekolah -0.152124 0.1414 -1.07576
## Pengeluaran -1.204868 0.1729 -6.96781
## UHH -0.541104 0.1214 -4.45557
## Sanitasi -0.179214 0.1451 -1.23482
## AirMinum 0.518370 0.1226 4.22943
## Pengangguran -0.006251 0.1380 -0.04528
## TPAK 0.292977 0.1299 2.25548
## PDRB -0.023296 0.1746 -0.13346
##
## Intercepts:
## Value Std. Error t value
## Rendah|Sedang -1.0052 0.1155 -8.7027
## Sedang|Tinggi 0.9925 0.1161 8.5466
##
## Residual Deviance: 883.5517
## AIC: 903.5517
# =========================
# UJI SIGNIFIKANSI PARAMETER (UJI WALD - PARSIAL)
# =========================
koef <- coef(summary(model))
p_value <- pnorm(abs(koef[, "t value"]), lower.tail = FALSE) * 2
tabel_wald <- cbind(koef, "p_value" = round(p_value, 4))
cat("\n=== Uji Wald (Parsial) ===\n")
##
## === Uji Wald (Parsial) ===
print(round(tabel_wald, 4))
## Value Std. Error t value p_value
## Sekolah -0.1521 0.1414 -1.0758 0.2820
## Pengeluaran -1.2049 0.1729 -6.9678 0.0000
## UHH -0.5411 0.1214 -4.4556 0.0000
## Sanitasi -0.1792 0.1451 -1.2348 0.2169
## AirMinum 0.5184 0.1226 4.2294 0.0000
## Pengangguran -0.0063 0.1380 -0.0453 0.9639
## TPAK 0.2930 0.1299 2.2555 0.0241
## PDRB -0.0233 0.1746 -0.1335 0.8938
## Rendah|Sedang -1.0052 0.1155 -8.7027 0.0000
## Sedang|Tinggi 0.9925 0.1161 8.5466 0.0000
cat("\nVariabel signifikan (p-value < 0.05):\n")
##
## Variabel signifikan (p-value < 0.05):
print(rownames(tabel_wald)[tabel_wald[, "p_value"] < 0.05 &
!rownames(tabel_wald) %in% c("Rendah|Sedang", "Sedang|Tinggi")])
## [1] "Pengeluaran" "UHH" "AirMinum" "TPAK"
# =========================
# UJI SIGNIFIKANSI MODEL (UJI G - SIMULTAN)
# =========================
model_null <- polr(Kategori ~ 1, data = data_model, Hess = TRUE)
G_hitung <- as.numeric(-2 * (logLik(model_null) - logLik(model)))
df_g <- length(coef(model))
p_G <- pchisq(G_hitung, df = df_g, lower.tail = FALSE)
cat("\n=== Uji G (Simultan) ===\n")
##
## === Uji G (Simultan) ===
cat("G hitung :", round(G_hitung, 4), "\n")
## G hitung : 245.8178
cat("df :", df_g, "\n")
## df : 8
cat("p-value :", format.pval(p_G, digits = 4), "\n")
## p-value : < 2.2e-16
cat("Kesimpulan:", ifelse(p_G < 0.05, "Tolak H0, minimal ada satu variabel X yang berpengaruh signifikan terhadap Y", "Gagal tolak H0, tidak ada variabel X yang berpengaruh signifikan terhadap Y"), "\n")
## Kesimpulan: Tolak H0, minimal ada satu variabel X yang berpengaruh signifikan terhadap Y
# =========================
# UJI ASUMSI PROPORTIONAL ODDS (BRANT TEST)
# =========================
cat("\n=== Uji Proportional Odds (Brant Test) ===\n")
##
## === Uji Proportional Odds (Brant Test) ===
brant_test <- brant(model)
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 21.7 8 0.01
## Sekolah 0.28 1 0.6
## Pengeluaran 0.18 1 0.67
## UHH 3.07 1 0.08
## Sanitasi 1.55 1 0.21
## AirMinum 0.33 1 0.56
## Pengangguran 1.97 1 0.16
## TPAK 4.65 1 0.03
## PDRB 1.65 1 0.2
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
print(brant_test)
## X2 df probability
## Omnibus 21.7039635 8 0.005494826
## Sekolah 0.2823169 1 0.595186643
## Pengeluaran 0.1828822 1 0.668907945
## UHH 3.0723144 1 0.079636020
## Sanitasi 1.5452388 1 0.213839698
## AirMinum 0.3321262 1 0.564409805
## Pengangguran 1.9696357 1 0.160486533
## TPAK 4.6494064 1 0.031064287
## PDRB 1.6496895 1 0.199001197
# H0: asumsi proportional odds terpenuhi
# Jika p-value (Omnibus) > 0.05, asumsi proportional odds terpenuhi sehingga model regresi logistik ordinal layak digunakan.
# =========================
# ODDS RATIO DAN INTERPRETASI
# =========================
OR <- exp(coef(model))
CI_OR <- exp(confint(model))
## Waiting for profiling to be done...
tabel_OR <- data.frame(
Variabel = names(OR),
OR = round(OR, 4),
CI_2.5 = round(CI_OR[, 1], 4),
CI_97.5 = round(CI_OR[, 2], 4)
)
cat("\n=== Odds Ratio (per kenaikan 1 SD) & CI 95% ===\n")
##
## === Odds Ratio (per kenaikan 1 SD) & CI 95% ===
print(tabel_OR, row.names = FALSE)
## Variabel OR CI_2.5 CI_97.5
## Sekolah 0.8589 0.6502 1.1327
## Pengeluaran 0.2997 0.2123 0.4185
## UHH 0.5821 0.4573 0.7366
## Sanitasi 0.8359 0.6277 1.1097
## AirMinum 1.6793 1.3229 2.1407
## Pengangguran 0.9938 0.7582 1.3033
## TPAK 1.3404 1.0416 1.7344
## PDRB 0.9770 0.6743 1.3308
# =========================
# EVALUASI MODEL (KLASIFIKASI & AKURASI)
# =========================
prediksi <- predict(model, type = "class")
tabel_klasifikasi <- table(Aktual = data_model$Kategori, Prediksi = prediksi)
cat("\n=== Tabel Klasifikasi (Aktual vs Prediksi) ===\n")
##
## === Tabel Klasifikasi (Aktual vs Prediksi) ===
print(tabel_klasifikasi)
## Prediksi
## Aktual Rendah Sedang Tinggi
## Rendah 107 47 18
## Sedang 46 88 37
## Tinggi 11 46 114
akurasi <- sum(diag(tabel_klasifikasi)) / sum(tabel_klasifikasi)
cat("\nAkurasi model:", percent(akurasi), "\n")
##
## Akurasi model: 60%
# tabel hasil prediksi per kab/kota (opsional, untuk lampiran laporan)
hasil_prediksi <- data.frame(
Provinsi = data$Provinsi,
KabKota = data$KabKota,
Kemiskinan_persen = data$Kemiskinan,
Kategori_Aktual = data_model$Kategori,
Kategori_Prediksi = prediksi
)
write.csv(hasil_prediksi, "hasil_prediksi_kabkota.csv", row.names = FALSE)
cat("\n=== SELESAI ===\n")
##
## === SELESAI ===