Nama: Stefanie Aurelia Clarissa Harisusanti
NPM: 140610240096
Analisis ini mengevaluasi hubungan IPM, TPT, dan kemiskinan menggunakan empat pendekatan regresi kategorik. Hasil menunjukkan model biner dan multinomial mengalami perfect separation/perfect classification sehingga sangat baik untuk klasifikasi tetapi lemah untuk inferensi. Model ordinal menghasilkan akurasi yang lebih realistis (53,125%) sehingga paling sesuai untuk menjelaskan hubungan substantif antar variabel. Model Poisson menunjukkan overdispersion ekstrem (11382,76) sehingga Negative Binomial direkomendasikan.
Temuan utama menunjukkan pembangunan manusia, kemiskinan, dan pengangguran memiliki keterkaitan yang kompleks dalam pembangunan wilayah Indonesia.
Data BPS 2024 digunakan untuk menggambarkan kondisi sosial ekonomi provinsi di Indonesia. Variabel yang digunakan meliputi Indeks Pembangunan Manusia (IPM), Tingkat Pengangguran Terbuka (TPT), persentase penduduk miskin, dan jumlah penduduk miskin. Analisis data kategorik sangat relevan karena banyak fenomena pembangunan direpresentasikan dalam bentuk kategori.
Model digunakan ketika respons terdiri dari dua kategori.
Model digunakan untuk respons nominal lebih dari dua kategori.
Model digunakan untuk kategori yang memiliki urutan alami.
Model digunakan untuk data cacah.
Data berasal dari publikasi BPS tahun 2024 yang memuat indikator pembangunan provinsi di Indonesia.
Pendekatan yang digunakan adalah regresi logistik biner, multinomial, ordinal, dan Poisson.
# ============================================================
# Analisis Regresi pada Data BPS 2024 (35 Provinsi)
# ============================================================
library(tidyverse)
library(knitr)
library(kableExtra)
library(nnet)
library(MASS)
library(scales)
theme_brown <- function() {
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", color = "#243142", size = 13),
plot.subtitle = element_text(color = "#5c6b7a", size = 10),
axis.title = element_text(color = "#243142"),
axis.text = element_text(color = "#3d4f5c"),
legend.title = element_text(face = "bold", color = "#243142"),
panel.grid.minor = element_blank()
)
}
# ============================================================
# DATA BPS 2024 (35 provinsi)
# ============================================================
data_bps <- tibble::tribble(
~provinsi, ~IPM, ~TPT, ~persen_miskin, ~jumlah_penduduk_miskin,
"Aceh", 75.36, 5.75, 12.64, 718.96,
"Sumatera Utara", 75.76, 5.60, 7.19, 1110.92,
"Sumatera Barat", 76.43, 5.75, 5.42, 315.43,
"Riau", 75.67, 3.70, 6.36, 473.04,
"Jambi", 74.36, 4.48, 7.26, 272.70,
"Sumatera Selatan", 73.84, 3.86, 10.51, 948.84,
"Bengkulu", 74.91, 3.11, 12.52, 261.15,
"Lampung", 73.13, 4.19, 10.62, 939.30,
"Jawa Barat", 74.92, 6.75, 7.08, 3668.35,
"Jawa Tengah", 73.87, 4.78, 9.58, 3396.34,
"Jawa Timur", 75.35, 4.19, 9.56, 3893.82,
"Banten", 76.35, 6.68, 5.70, 777.49,
"Bali", 78.63, 1.79, 3.80, 176.21,
"Nusa Tenggara Barat", 73.10, 2.73, 11.91, 658.60,
"Nusa Tenggara Timur", 69.14, 3.02, 19.02, 1107.94,
"Kalimantan Barat", 71.19, 4.86, 6.25, 333.99,
"Kalimantan Tengah", 74.28, 4.01, 5.26, 149.24,
"Kalimantan Selatan", 75.19, 4.20, 4.02, 180.20,
"Kalimantan Timur", 78.79, 5.14, 5.51, 211.88,
"Kalimantan Utara", 73.41, 3.90, 5.38, 41.11,
"Sulawesi Utara", 75.68, 5.85, 6.70, 173.30,
"Sulawesi Tengah", 72.24, 2.94, 11.04, 358.33,
"Sulawesi Selatan", 75.18, 4.19, 7.77, 711.77,
"Sulawesi Tenggara", 73.62, 3.09, 10.63, 305.27,
"Gorontalo", 72.01, 3.13, 13.87, 170.03,
"Sulawesi Barat", 70.46, 2.68, 10.71, 155.91,
"Maluku", 73.40, 6.11, 15.78, 293.99,
"Maluku Utara", 71.84, 4.03, 6.03, 79.69,
"Papua Barat", 67.69, 4.13, 21.09, 108.28,
"Papua Barat Daya", 69.65, 6.48, 16.95, 96.81,
"Papua", 73.83, 6.48, 18.09, 161.07,
"Indonesia", 75.02, 4.91, 8.57, 24054.72
)
dplyr::glimpse(data_bps)
# ============================================================
# BAGIAN 1: REGRESI LOGISTIK BINER
# ============================================================
# ---- Membuat variabel respons biner -------------------------
data_bin <- data_bps %>%
mutate(
miskin_tinggi = factor(
ifelse(persen_miskin > median(persen_miskin), "Tinggi", "Rendah"),
levels = c("Rendah", "Tinggi")
)
)
dplyr::glimpse(data_bin)
# ---- Distribusi variabel respons ----------------------------
data_bin %>%
count(miskin_tinggi) %>%
mutate(proporsi = n / sum(n)) %>%
kable(
digits = 3,
caption = "Distribusi kategori kemiskinan provinsi (Biner)"
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
)
data_bin %>%
count(miskin_tinggi) %>%
mutate(proporsi = n / sum(n)) %>%
ggplot(aes(x = miskin_tinggi, y = proporsi, fill = miskin_tinggi)) +
geom_col(width = 0.65, alpha = 0.94) +
geom_text(aes(label = percent(proporsi, accuracy = 0.1)),
vjust = -0.4, fontface = "bold", color = "#243142") +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.75)) +
scale_fill_manual(values = c("#2f7f73", "#d18b2f")) +
labs(
title = "Distribusi Kategori Kemiskinan",
subtitle = "Respons biner: Rendah vs Tinggi berdasarkan median persen miskin.",
x = NULL,
y = "Proporsi"
) +
theme_brown() +
theme(legend.position = "none")
# ---- Fitting model regresi logistik biner -------------------
fit_bin <- glm(
miskin_tinggi ~ IPM + TPT + persen_miskin,
data = data_bin,
family = binomial(link = "logit")
)
summary(fit_bin)
# ---- Tabel ringkasan koefisien ------------------------------
bin_coef <- as.data.frame(coef(summary(fit_bin))) %>%
tibble::rownames_to_column("parameter") %>%
dplyr::rename(
estimate = Estimate,
std_error = `Std. Error`,
z_value = `z value`,
p_value = `Pr(>|z|)`
) %>%
mutate(
OR = exp(estimate),
CI_low = exp(estimate - 1.96 * std_error),
CI_high = exp(estimate + 1.96 * std_error)
)
bin_coef %>%
mutate(
across(
c(estimate, std_error, z_value, p_value, OR, CI_low, CI_high),
~ round(.x, 4)
)
) %>%
kable(
caption = "Ringkasan hasil regresi logistik biner",
col.names = c(
"Parameter", "Estimate", "SE", "z-value", "p-value",
"OR", "CI 2.5%", "CI 97.5%"
)
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
)
# ---- Prediksi & evaluasi model ------------------------------
pred_prob_bin <- predict(fit_bin, type = "response")
data_bin_pred <- data_bin %>%
mutate(
prob_tinggi = pred_prob_bin,
prediksi = factor(
ifelse(prob_tinggi >= 0.5, "Tinggi", "Rendah"),
levels = c("Rendah", "Tinggi")
)
)
head(data_bin_pred)
conf_bin <- table(
Aktual = data_bin_pred$miskin_tinggi,
Prediksi = data_bin_pred$prediksi
)
conf_bin
accuracy_bin <- sum(diag(conf_bin)) / sum(conf_bin)
accuracy_bin
# ---- Visualisasi probabilitas prediksi vs IPM ---------------
grid_bin <- expand.grid(
IPM = seq(min(data_bin$IPM), max(data_bin$IPM), length.out = 120),
TPT = mean(data_bin$TPT),
persen_miskin = mean(data_bin$persen_miskin)
)
pred_bin_link <- predict(fit_bin, newdata = grid_bin, type = "link", se.fit = TRUE)
grid_bin_plot <- grid_bin %>%
mutate(
fit_link = pred_bin_link$fit,
se_link = pred_bin_link$se.fit,
prob = plogis(fit_link),
prob_low = plogis(fit_link - 1.96 * se_link),
prob_high = plogis(fit_link + 1.96 * se_link)
)
ggplot(grid_bin_plot, aes(x = IPM, y = prob)) +
geom_ribbon(aes(ymin = prob_low, ymax = prob_high), fill = "#2f7f73", alpha = 0.16) +
geom_line(color = "#2f7f73", linewidth = 1.35) +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "#d18b2f", linewidth = 0.8) +
scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
labs(
title = "Prediksi Probabilitas Kemiskinan Tinggi vs IPM",
subtitle = "TPT dan persen miskin ditahan pada nilai rata-rata; garis putus = ambang 50%.",
x = "Indeks Pembangunan Manusia (IPM)",
y = "P(Kemiskinan Tinggi)"
) +
theme_brown()
# ============================================================
# BAGIAN 2: REGRESI LOGISTIK MULTINOMIAL
# ============================================================
# ---- Membuat variabel respons multinomial -------------------
data_multi <- data_bps %>%
mutate(
status_wilayah = factor(
case_when(
IPM >= 78 & persen_miskin <= 8 ~ "Maju",
IPM < 72 & persen_miskin >= 12 ~ "Tertinggal",
TRUE ~ "Berkembang"
),
levels = c("Berkembang", "Maju", "Tertinggal")
)
)
dplyr::glimpse(data_multi)
# ---- Distribusi variabel respons ----------------------------
data_multi %>%
count(status_wilayah) %>%
mutate(proporsi = n / sum(n)) %>%
kable(
digits = 3,
caption = "Distribusi status wilayah provinsi (Multinomial)"
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
)
data_multi %>%
count(status_wilayah) %>%
mutate(proporsi = n / sum(n)) %>%
ggplot(aes(x = status_wilayah, y = proporsi, fill = status_wilayah)) +
geom_col(width = 0.65, alpha = 0.94) +
geom_text(aes(label = percent(proporsi, accuracy = 0.1)),
vjust = -0.4, fontface = "bold", color = "#243142") +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.85)) +
scale_fill_manual(values = c("#2f7f73", "#5568b8", "#d18b2f")) +
labs(
title = "Distribusi Status Wilayah",
subtitle = "Respons multinomial: kategori tidak memiliki urutan alami.",
x = NULL,
y = "Proporsi"
) +
theme_brown() +
theme(legend.position = "none")
# ---- Fitting model regresi logistik multinomial -------------
fit_multi <- nnet::multinom(
status_wilayah ~ IPM + TPT + persen_miskin,
data = data_multi,
trace = FALSE
)
summary(fit_multi)
# ---- Tabel ringkasan koefisien ------------------------------
multi_sum <- summary(fit_multi)
coef_multi <- as.data.frame(multi_sum$coefficients)
se_multi <- as.data.frame(multi_sum$standard.errors)
coef_long <- coef_multi %>%
tibble::rownames_to_column("kategori") %>%
tidyr::pivot_longer(
cols = -kategori,
names_to = "variabel",
values_to = "estimate"
)
se_long <- se_multi %>%
tibble::rownames_to_column("kategori") %>%
tidyr::pivot_longer(
cols = -kategori,
names_to = "variabel",
values_to = "std_error"
)
result_multi <- coef_long %>%
left_join(se_long, by = c("kategori", "variabel")) %>%
mutate(
z_value = estimate / std_error,
p_value = 2 * (1 - pnorm(abs(z_value))),
RRR = exp(estimate),
CI_low = exp(estimate - 1.96 * std_error),
CI_high = exp(estimate + 1.96 * std_error)
)
result_multi %>%
mutate(
across(
c(estimate, std_error, z_value, p_value, RRR, CI_low, CI_high),
~ round(.x, 4)
)
) %>%
kable(
caption = "Ringkasan hasil regresi logistik multinomial",
col.names = c(
"Kategori", "Variabel", "Estimate", "SE", "z-value",
"p-value", "RRR", "CI 2.5%", "CI 97.5%"
)
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
)
# ---- Prediksi & evaluasi model ------------------------------
pred_prob_multi <- predict(fit_multi, type = "probs")
head(pred_prob_multi)
data_multi_pred <- data_multi %>%
bind_cols(as.data.frame(pred_prob_multi)) %>%
mutate(
prediksi = predict(fit_multi, type = "class")
)
head(data_multi_pred)
conf_multi <- table(
Aktual = data_multi_pred$status_wilayah,
Prediksi = data_multi_pred$prediksi
)
conf_multi
accuracy_multi <- sum(diag(conf_multi)) / sum(conf_multi)
accuracy_multi
# ---- Visualisasi probabilitas prediksi vs IPM ---------------
grid_multi <- expand.grid(
IPM = seq(min(data_multi$IPM), max(data_multi$IPM), length.out = 100),
TPT = mean(data_multi$TPT),
persen_miskin = mean(data_multi$persen_miskin)
)
grid_prob_multi <- predict(fit_multi, newdata = grid_multi, type = "probs")
grid_multi_plot <- grid_multi %>%
bind_cols(as.data.frame(grid_prob_multi)) %>%
tidyr::pivot_longer(
cols = c("Berkembang", "Maju", "Tertinggal"),
names_to = "status_wilayah",
values_to = "probabilitas"
)
ggplot(
grid_multi_plot,
aes(x = IPM, y = probabilitas, color = status_wilayah)
) +
geom_line(linewidth = 1.35) +
scale_y_continuous(labels = percent_format()) +
scale_color_manual(values = c("#2f7f73", "#5568b8", "#d18b2f")) +
labs(
title = "Prediksi Probabilitas Status Wilayah",
subtitle = "TPT dan persen miskin ditahan pada nilai rata-rata.",
x = "Indeks Pembangunan Manusia (IPM)",
y = "Probabilitas prediksi",
color = "Status wilayah"
) +
theme_brown()
# ============================================================
# BAGIAN 3: REGRESI LOGISTIK ORDINAL
# ============================================================
# ---- Membuat variabel respons ordinal -----------------------
data_ord <- data_bps %>%
mutate(
kategori_ipm = cut(
IPM,
breaks = quantile(IPM, probs = c(0, .25, .5, .75, 1)),
include.lowest = TRUE,
labels = c("Rendah", "Sedang", "Tinggi", "Sangat Tinggi"),
ordered_result = TRUE
)
)
dplyr::glimpse(data_ord)
# ---- Distribusi variabel respons ----------------------------
data_ord %>%
count(kategori_ipm) %>%
mutate(proporsi = n / sum(n)) %>%
kable(
digits = 3,
caption = "Distribusi kategori IPM provinsi (Ordinal)"
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
)
data_ord %>%
count(kategori_ipm) %>%
mutate(proporsi = n / sum(n)) %>%
ggplot(aes(x = kategori_ipm, y = proporsi, fill = kategori_ipm)) +
geom_col(width = 0.65, alpha = 0.94) +
geom_text(aes(label = percent(proporsi, accuracy = 0.1)),
vjust = -0.4, fontface = "bold", color = "#243142") +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.55)) +
scale_fill_manual(values = c("#2f7f73", "#5568b8", "#d18b2f", "#b84f5a")) +
labs(
title = "Distribusi Kategori IPM",
subtitle = "Respons ordinal: kategori memiliki urutan alami.",
x = NULL,
y = "Proporsi"
) +
theme_brown() +
theme(legend.position = "none")
# ---- Fitting model regresi logistik ordinal -----------------
fit_ord <- MASS::polr(
kategori_ipm ~ TPT + persen_miskin,
data = data_ord,
method = "logistic",
Hess = TRUE
)
summary(fit_ord)
# ---- Tabel ringkasan koefisien ------------------------------
ord_coef <- coef(summary(fit_ord))
result_ord <- as.data.frame(ord_coef) %>%
tibble::rownames_to_column("parameter") %>%
dplyr::rename(
estimate = Value,
std_error = `Std. Error`,
t_value = `t value`
) %>%
mutate(
p_value = 2 * (1 - pnorm(abs(t_value))),
jenis = ifelse(grepl("\\|", parameter), "Cutpoint", "Koefisien"),
estimate_slide = ifelse(jenis == "Koefisien", -estimate, estimate),
OR_polr = ifelse(jenis == "Koefisien", exp(estimate), NA_real_),
OR_slide = ifelse(jenis == "Koefisien", exp(estimate_slide), NA_real_),
CI_low_polr = ifelse(jenis == "Koefisien", exp(estimate - 1.96 * std_error), NA_real_),
CI_high_polr = ifelse(jenis == "Koefisien", exp(estimate + 1.96 * std_error), NA_real_)
)
result_ord %>%
mutate(
across(
c(estimate, estimate_slide, std_error, t_value, p_value,
OR_polr, OR_slide, CI_low_polr, CI_high_polr),
~ round(.x, 4)
)
) %>%
kable(
caption = "Ringkasan hasil regresi logistik ordinal",
col.names = c(
"Parameter", "Estimate polr", "SE", "z/t-value", "p-value",
"Jenis", "Estimate slide", "OR polr", "OR slide",
"CI polr 2.5%", "CI polr 97.5%"
)
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
)
# ---- Prediksi & evaluasi model ------------------------------
pred_prob_ord <- predict(fit_ord, type = "probs")
head(pred_prob_ord)
data_ord_pred <- data_ord %>%
bind_cols(as.data.frame(pred_prob_ord)) %>%
mutate(
prediksi = predict(fit_ord, type = "class")
)
head(data_ord_pred)
conf_ord <- table(
Aktual = data_ord_pred$kategori_ipm,
Prediksi = data_ord_pred$prediksi
)
conf_ord
accuracy_ord <- sum(diag(conf_ord)) / sum(conf_ord)
accuracy_ord
# ---- Visualisasi probabilitas prediksi vs persen_miskin -----
grid_ord <- expand.grid(
TPT = mean(data_ord$TPT),
persen_miskin = seq(
min(data_ord$persen_miskin),
max(data_ord$persen_miskin),
length.out = 120
)
)
grid_prob_ord <- predict(fit_ord, newdata = grid_ord, type = "probs")
grid_ord_plot <- grid_ord %>%
bind_cols(as.data.frame(grid_prob_ord)) %>%
tidyr::pivot_longer(
cols = c("Rendah", "Sedang", "Tinggi", "Sangat Tinggi"),
names_to = "kategori_ipm",
values_to = "probabilitas"
) %>%
mutate(
kategori_ipm = factor(
kategori_ipm,
levels = c("Rendah", "Sedang", "Tinggi", "Sangat Tinggi")
)
)
ggplot(
grid_ord_plot,
aes(x = persen_miskin, y = probabilitas, color = kategori_ipm)
) +
geom_line(linewidth = 1.25) +
scale_y_continuous(labels = percent_format()) +
scale_color_manual(values = c("#2f7f73", "#5568b8", "#d18b2f", "#b84f5a")) +
labs(
title = "Prediksi Probabilitas Kategori IPM",
subtitle = "TPT ditahan pada nilai rata-rata.",
x = "Persentase penduduk miskin (%)",
y = "Probabilitas prediksi",
color = "Kategori IPM"
) +
theme_brown()
# ---- Plot parallel lines (cumulative logit) -----------------
eps <- 1e-6
grid_ord_cumlogit <- grid_ord %>%
bind_cols(as.data.frame(grid_prob_ord)) %>%
mutate(
`P(Y <= Rendah)` = Rendah,
`P(Y <= Sedang)` = Rendah + Sedang,
`P(Y <= Tinggi)` = Rendah + Sedang + Tinggi
) %>%
tidyr::pivot_longer(
cols = c(`P(Y <= Rendah)`, `P(Y <= Sedang)`, `P(Y <= Tinggi)`),
names_to = "batas_kumulatif",
values_to = "prob_kumulatif"
) %>%
mutate(
prob_kumulatif = pmin(pmax(prob_kumulatif, eps), 1 - eps),
cumulative_logit = qlogis(prob_kumulatif)
)
ggplot(
grid_ord_cumlogit,
aes(x = persen_miskin, y = cumulative_logit, color = batas_kumulatif)
) +
geom_line(linewidth = 1.25) +
scale_color_manual(values = c(
"P(Y <= Rendah)" = "#2f7f73",
"P(Y <= Sedang)" = "#d18b2f",
"P(Y <= Tinggi)" = "#5568b8"
)) +
labs(
title = "Parallel Lines pada Model Ordinal",
subtitle = "Garis yang sejajar adalah cumulative logit, bukan probabilitas kategori tunggal.",
x = "Persentase penduduk miskin (%)",
y = "Logit peluang kumulatif",
color = "Batas kumulatif"
) +
theme_brown()
# ============================================================
# BAGIAN 4: REGRESI POISSON
# ============================================================
dplyr::glimpse(data_bps)
# ---- Distribusi variabel respons ----------------------------
data_bps %>%
mutate(bin_penduduk_miskin = round(jumlah_penduduk_miskin / 500) * 500) %>%
count(bin_penduduk_miskin) %>%
ggplot(aes(x = bin_penduduk_miskin, y = n)) +
geom_col(width = 450, fill = "#2f7f73", alpha = 0.92) +
scale_x_continuous(breaks = scales::pretty_breaks()) +
labs(
title = "Distribusi Jumlah Penduduk Miskin",
subtitle = "Respons Poisson berupa hitungan (dalam ribuan jiwa).",
x = "Jumlah penduduk miskin (ribu jiwa)",
y = "Frekuensi"
) +
theme_brown()
# ---- Fitting model regresi Poisson --------------------------
fit_pois <- glm(
jumlah_penduduk_miskin ~ IPM + TPT,
data = data_bps,
family = poisson(link = "log")
)
summary(fit_pois)
# ---- Tabel ringkasan koefisien ------------------------------
pois_coef <- as.data.frame(coef(summary(fit_pois))) %>%
tibble::rownames_to_column("parameter") %>%
dplyr::rename(
estimate = Estimate,
std_error = `Std. Error`,
z_value = `z value`,
p_value = `Pr(>|z|)`
) %>%
mutate(
IRR = exp(estimate),
CI_low = exp(estimate - 1.96 * std_error),
CI_high = exp(estimate + 1.96 * std_error),
perubahan_persen = 100 * (IRR - 1)
)
pois_coef %>%
mutate(
across(
c(estimate, std_error, z_value, p_value, IRR, CI_low, CI_high, perubahan_persen),
~ round(.x, 4)
)
) %>%
kable(
caption = "Ringkasan hasil regresi Poisson",
col.names = c(
"Parameter", "Estimate", "SE", "z-value", "p-value",
"IRR", "CI 2.5%", "CI 97.5%", "Perubahan (%)"
)
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
)
# ---- Visualisasi prediksi rate vs IPM -----------------------
grid_pois <- expand.grid(
IPM = seq(min(data_bps$IPM), max(data_bps$IPM), length.out = 120),
TPT = mean(data_bps$TPT)
)
pred_pois <- predict(fit_pois, newdata = grid_pois, type = "link", se.fit = TRUE)
grid_pois_plot <- grid_pois %>%
mutate(
fit_link = pred_pois$fit,
se_link = pred_pois$se.fit,
rate = exp(fit_link),
rate_low = exp(fit_link - 1.96 * se_link),
rate_high = exp(fit_link + 1.96 * se_link)
)
ggplot(grid_pois_plot, aes(x = IPM, y = rate)) +
geom_ribbon(aes(ymin = rate_low, ymax = rate_high), fill = "#2f7f73", alpha = 0.16) +
geom_line(color = "#2f7f73", linewidth = 1.35) +
labs(
title = "Prediksi Jumlah Penduduk Miskin vs IPM",
subtitle = "TPT ditahan pada nilai rata-rata; pita = interval kepercayaan 95%.",
x = "Indeks Pembangunan Manusia (IPM)",
y = "Prediksi jumlah penduduk miskin (ribu jiwa)"
) +
theme_brown()
# ---- Cek overdispersion -------------------------------------
dispersion_pois <- sum(residuals(fit_pois, type = "pearson")^2) / df.residual(fit_pois)
tibble::tibble(
`Dispersion Pearson` = dispersion_pois,
`Interpretasi ringkas` = dplyr::case_when(
dispersion_pois < 1.5 ~ "Tidak ada indikasi overdispersion berat",
dispersion_pois < 2.5 ~ "Ada indikasi overdispersion sedang",
TRUE ~ "Ada indikasi overdispersion kuat"
)
) %>%
mutate(`Dispersion Pearson` = round(`Dispersion Pearson`, 3)) %>%
kable(
caption = "Indikasi overdispersion pada model Poisson"
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
)
Residual Deviance = 3.3934e-09
AIC = 8
Accuracy = 1.00
Nilai residual deviance yang hampir nol menunjukkan model memisahkan kategori secara sempurna. Akurasi 100% mengindikasikan perfect classification. Dalam literatur regresi logistik kondisi ini dikenal sebagai perfect separation.
Odds Ratio IPM = exp(-11.13), menunjukkan kenaikan IPM menurunkan peluang masuk kategori kemiskinan tinggi. TPT juga memiliki arah negatif. Sebaliknya persentase miskin berasosiasi positif terhadap peluang masuk kelompok kemiskinan tinggi.
Hasil menunjukkan daerah dengan kualitas pembangunan manusia lebih baik cenderung memiliki probabilitas lebih rendah untuk masuk kategori kemiskinan tinggi.
Residual Deviance = 0.07875826
AIC = 16.07876
Accuracy = 1.00
RRR IPM > 1 sehingga peningkatan IPM meningkatkan peluang wilayah masuk kategori maju dibanding berkembang.
Kenaikan IPM secara drastis menurunkan peluang wilayah menjadi tertinggal. Sebaliknya peningkatan TPT dan kemiskinan meningkatkan peluang wilayah masuk kategori tertinggal.
Akurasi 100% menunjukkan perfect classification. Model sangat baik mengklasifikasikan data yang tersedia namun berpotensi terlalu menyesuaikan data.
Accuracy = 0.53125
OR(TPT)=exp(0.6417) ≈ 1.90.
Artinya kenaikan satu satuan TPT meningkatkan odds berada pada kategori IPM yang lebih rendah sekitar 90%.
OR(persen_miskin)=exp(-0.3497) ≈ 0.705.
Artinya peningkatan kemiskinan menurunkan peluang berada pada kategori IPM yang lebih tinggi.
Cutpoint menggambarkan batas laten antar kategori IPM.
Walaupun akurasi tidak setinggi model lain, model ordinal lebih realistis karena mempertahankan struktur urutan kategori. Oleh karena itu model ini paling layak digunakan untuk inferensi kebijakan.
Intercept = -3.395164
IPM = 0.130984
TPT = 0.201902
Seluruh p-value < 0.001
IRR(IPM)=exp(0.130984) ≈ 1.14.
Setiap kenaikan satu poin IPM diasosiasikan dengan kenaikan ekspektasi jumlah penduduk miskin sebesar sekitar 14%. Hasil ini perlu dibaca hati-hati karena dipengaruhi ukuran populasi provinsi.
IRR(TPT)=exp(0.201902) ≈ 1.22.
Setiap kenaikan satu poin TPT meningkatkan ekspektasi jumlah penduduk miskin sekitar 22%.
Dispersion Pearson = 11382.76.
Nilai ini jauh melebihi 1 sehingga terjadi overdispersion ekstrem. Asumsi Poisson tidak terpenuhi.
Model Negative Binomial direkomendasikan karena lebih mampu menangani varians yang jauh lebih besar daripada rata-rata.
| Model | Kelebihan | Keterbatasan |
|---|---|---|
| Biner | Akurasi 100% | Perfect separation |
| Multinomial | Akurasi 100% | Perfect classification |
| Ordinal | Inferensi realistis | Akurasi sedang |
| Poisson | Data cacah | Overdispersion ekstrem |
Model ordinal dinilai paling tepat untuk penjelasan substantif, sedangkan model Poisson memerlukan pengembangan ke Negative Binomial.
Peningkatan IPM merupakan strategi utama pembangunan daerah. Pengurangan kemiskinan dan pengangguran harus dilakukan secara simultan. Daerah tertinggal memerlukan intervensi pendidikan, kesehatan, dan penciptaan lapangan kerja yang lebih kuat.
Analisis menunjukkan bahwa hubungan IPM, kemiskinan, dan pengangguran sangat erat dalam pembangunan daerah Indonesia. Model biner dan multinomial memberikan klasifikasi sempurna namun kurang ideal untuk inferensi. Model ordinal merupakan model yang paling realistis. Model Poisson menunjukkan overdispersion ekstrem sehingga Negative Binomial direkomendasikan.
Badan Pusat Statistik. (2024). Indikator kesejahteraan rakyat Indonesia.
Agresti, A. (2018). An Introduction to Categorical Data Analysis.