Analisis Regresi Biner, Multinomial, Ordinal, dan Poisson
Food Crisis, Environment Priority, Government Confidence, dan Violence Case Count
Febrianggi Caesar Immanuel
01 June 2026
1 Pendahuluan
1.1 Gambaran Umum Analisis
Laporan ini menerapkan empat model regresi berdasarkan skala variabel respons: regresi logistik biner untuk respons dua kategori, regresi logistik multinomial untuk respons nominal lebih dari dua kategori, regresi logistik ordinal untuk respons berurutan, dan regresi Poisson untuk respons hitungan. Alur analisis mengikuti materi regresi kategorik dan Poisson: pemahaman data, estimasi model, uji simultan \(G^2\), uji Wald, pemilihan model terbaik, interpretasi variabel yang mendukung respons, dan prediksi menggunakan model terbaik.
1.2 Library yang Digunakan
Library berikut digunakan untuk membaca data, mengolah data, membangun model regresi, membuat visualisasi, dan menyajikan tabel pada dokumen HTML.
library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(nnet)
library(MASS)
library(broom)
library(knitr)
library(kableExtra)
library(scales)
library(tibble)
library(purrr)
library(stringr)
library(pROC)1.3 Catatan Khusus
Catatan 1 — Kriteria pemilihan model terbaik.
Kriteria pemilihan model terbaik dalam laporan ini mengutamakan
jumlah variabel X yang signifikan terlebih dahulu,
kemudian mempertimbangkan AIC terkecil. Artinya,
apabila satu model memiliki lebih banyak variabel signifikan meskipun
AIC-nya sedikit lebih besar, model tersebut tetap lebih diprioritaskan
daripada model lain yang AIC-nya lebih kecil tetapi kehilangan banyak
variabel signifikan.
Catatan 2 — Interpretasi hasil model.
Interpretasi koefisien, odds ratio, relative risk ratio, dan incidence
rate ratio dalam laporan ini bersifat asosiasi statistik, bukan klaim
sebab-akibat. Dengan demikian, variabel X yang signifikan dipahami
sebagai variabel yang mendukung atau berkaitan dengan variasi variabel Y
dalam model, bukan sebagai penyebab langsung.
2 Sumber Data
data_sources <- tibble(
bagian=c("Regresi logistik biner","Regresi logistik multinomial","Regresi logistik ordinal","Regresi Poisson tahunan"),
topik=c("Food Crisis","Environment Priority","Government Confidence","Violence Case Count"),
sumber=c("World Bank - World Development Indicators","World Values Survey Wave 7","World Values Survey Wave 7","UCDP Georeferenced Event Dataset Global"),
tautan=c("https://databank.worldbank.org/source/world-development-indicators","https://www.worldvaluessurvey.org/WVSDocumentationWV7.jsp","https://www.worldvaluessurvey.org/WVSDocumentationWV7.jsp","https://ucdp.uu.se/downloads/#ged_global"),
file_model_ready=c("wdi_binary_logit_model_ready.xlsx","multinomial_model_ready_sample.xlsx","ordinal_model_ready_sample.xlsx","poisson_model_ready.xlsx")
)
safe_table(data_sources, "Sumber data dan file model-ready yang digunakan", full_width=TRUE)3 Bagian I — Regresi Logistik Biner
3.1 Konsep model
Regresi logistik biner digunakan ketika respons hanya memiliki dua kemungkinan hasil. Model logit menghubungkan prediktor dengan log odds kejadian:
\[ \log\left(\frac{p_i}{1-p_i}\right)=\beta_0+\beta_1X_{1i}+\cdots+\beta_kX_{ki} \]
Uji simultan model dilakukan dengan statistik likelihood ratio:
\[ G^2=-2(\ell_0-\ell_1)=2(\ell_1-\ell_0) \]
Jika p-value dari \(G^2\) kecil, model penuh lebih baik daripada model kosong.
3.2 Import, struktur data, dan distribusi respons
binary_data <- readxl::read_excel(resolve_path("wdi_binary_logit_model_ready.xlsx")) |>
clean_model_data() |>
mutate(food_crisis=as.integer(food_crisis), food_crisis_label=factor(food_crisis, levels=c(0,1), labels=c("non_food_crisis","food_crisis")))
binary_dict <- make_dictionary(binary_data)
binary_summary <- binary_data |> count(food_crisis, food_crisis_label, name="n") |> mutate(proporsi=n/sum(n))
safe_table(binary_dict, "Kamus data biner setelah seleksi variabel analisis")safe_table(binary_summary, "Distribusi respons biner food_crisis", digits=3)binary_data |> count(food_crisis_label) |>
ggplot(aes(x=food_crisis_label, y=n, fill=food_crisis_label))+
geom_col(width=.68)+geom_text(aes(label=n), vjust=-.45, fontface="bold")+
scale_fill_manual(values=binary_palette)+
labs(title="Distribusi Respons Food Crisis", subtitle="Respons biner: 0 = non-food crisis dan 1 = food crisis.", x="Kategori respons", y="Jumlah observasi", fill="Kategori")+
theme_report()
Gambar 1. Distribusi kategori respons food_crisis.
Variabel food_crisis sudah sesuai untuk regresi logistik
biner karena hanya memiliki dua kategori. Distribusi kelas dibaca
sebelum pemodelan karena ketidakseimbangan kelas dapat membuat accuracy
terlihat baik meskipun model kurang sensitif dalam mengenali kelas
kejadian.
3.3 Statistik Deskriptif dan Eksplorasi Data
binary_descriptive <- numeric_descriptive(binary_data, response="food_crisis")
binary_group_summary <- binary_data |>
select(food_crisis_label, where(is.numeric)) |>
select(-food_crisis) |>
group_by(food_crisis_label) |>
summarise(across(everything(), ~mean(.x, na.rm=TRUE)), .groups="drop") |>
pivot_longer(-food_crisis_label, names_to="variabel", values_to="mean") |>
pivot_wider(names_from=food_crisis_label, values_from=mean)
safe_table(binary_descriptive, "Statistik deskriptif prediktor numerik pada data logistik biner", digits=3, full_width=TRUE)safe_table(binary_group_summary, "Rata-rata prediktor menurut status food_crisis", digits=3, full_width=TRUE)binary_plot_vars <- binary_x <- model_predictors(binary_data, response="food_crisis")
binary_plot_vars <- binary_plot_vars[seq_len(min(4, length(binary_plot_vars)))]
binary_data |>
select(food_crisis_label, all_of(binary_plot_vars)) |>
pivot_longer(-food_crisis_label, names_to="variabel", values_to="nilai") |>
ggplot(aes(x=food_crisis_label, y=nilai, fill=food_crisis_label))+
geom_boxplot(alpha=.75, outlier.alpha=.25)+
facet_wrap(~variabel, scales="free_y")+
scale_fill_manual(values=binary_palette)+
labs(title="Eksplorasi Prediktor terhadap Status Food Crisis", x="Status food crisis", y="Nilai prediktor", fill="Kategori")+
theme_report()
Gambar 2. Eksplorasi prediktor numerik utama menurut status food_crisis.
3.4 Pembagian data training dan testing
binary_split <- split_train_test(binary_data, strata="food_crisis", prop=.80)
binary_train <- binary_split$train; binary_test <- binary_split$test
binary_split_table <- bind_rows(binary_train |> count(food_crisis_label, name="n") |> mutate(subset="Training"), binary_test |> count(food_crisis_label, name="n") |> mutate(subset="Testing")) |> group_by(subset) |> mutate(proporsi=n/sum(n)) |> ungroup() |> select(subset, food_crisis_label, n, proporsi)
safe_table(binary_split_table, "Komposisi kelas food_crisis pada data training dan testing", digits=3)3.5 Kandidat model dan pemilihan model terbaik
binary_x <- model_predictors(binary_train, response="food_crisis")
binary_selection <- select_best_glm(binary_train, response="food_crisis", predictors=binary_x, family=binomial(link="logit"), effect_name="odds_ratio")
binary_best_fit <- binary_selection$best$fit
binary_best_null <- glm(food_crisis ~ 1, data=binary_train, family=binomial(link="logit"))
binary_best_wald <- binary_selection$best$wald
binary_best_g2 <- binary_selection$best$g2
binary_selected_vars <- binary_best_wald |> filter(p_value < 0.05) |> pull(variabel)
binary_model_selection_table <- binary_selection$summary |>
mutate(p_value_G2 = p_value_G2_tampil) |>
select(model_ke, formula, jumlah_prediktor, jumlah_variabel_signifikan, AIC, G2, p_value_G2)
binary_fit_table <- tibble(keterangan=c("Jumlah observasi training","Jumlah prediktor model terbaik","AIC model terbaik","Log-likelihood model terbaik","Log-likelihood model kosong","McFadden pseudo R-square"), nilai=c(nrow(binary_train), length(binary_selection$best$vars), AIC(binary_best_fit), as.numeric(logLik(binary_best_fit)), as.numeric(logLik(binary_best_null)), mcfadden_r2(binary_best_fit,binary_best_null)))
binary_best_wald_show <- binary_best_wald |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil)
binary_best_g2_show <- binary_best_g2 |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil)
safe_table(binary_model_selection_table, "Perbandingan kandidat model logistik biner berdasarkan jumlah variabel signifikan dan AIC", digits=4, full_width=TRUE)safe_table(binary_fit_table, "Ringkasan kecocokan model terbaik logistik biner", digits=4)safe_table(binary_best_g2_show, "Uji simultan G² model terbaik logistik biner", digits=4)safe_table(binary_best_wald_show, "Uji Wald parsial dan odds ratio model terbaik logistik biner", digits=4, full_width=TRUE)Model terbaik logistik biner adalah kandidat ke-4 karena memiliki 4
variabel signifikan dan AIC sebesar 997.248. Uji simultan \(G^2\) menghasilkan p-value < 0.001,
sehingga model terbaik memberikan peningkatan kecocokan dibanding model
kosong. Berdasarkan uji Wald, variabel X yang signifikan dan mendukung
food_crisis adalah: log_gdp_per_capita, population_growth,
unemployment, year. Dengan demikian, prediksi selanjutnya dilakukan
menggunakan model terbaik tersebut, bukan model penuh awal.
3.6 Uji Asumsi Model Terbaik
3.6.1 Multikolinearitas — VIF
Hipotesis/Kriteria:
\(H_0\): Tidak terdapat
multikolinearitas berat antarprediktor.
\(H_1\): Terdapat multikolinearitas
berat antarprediktor.
Kriteria praktis: VIF < 5 menunjukkan tidak ada indikasi
multikolinearitas berat; VIF 5–10 perlu perhatian; VIF > 10
mengindikasikan multikolinearitas berat.
binary_vif_table <- vif_table(binary_train, binary_selection$best$vars)
safe_table(binary_vif_table, "Uji multikolinearitas model terbaik logistik biner menggunakan VIF", digits=3, full_width=TRUE)3.6.2 Likelihood Ratio Test (LRT)
Hipotesis:
\(H_0\): Semua koefisien prediktor sama
dengan nol, sehingga model dengan prediktor tidak lebih baik daripada
model kosong.
\(H_1\): Minimal ada satu koefisien
prediktor tidak sama dengan nol, sehingga model dengan prediktor lebih
baik daripada model kosong.
binary_lrt_show <- binary_best_g2 |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil)
safe_table(binary_lrt_show, "Likelihood Ratio Test atau uji G² model terbaik logistik biner", digits=4, full_width=TRUE)3.6.3 Hosmer-Lemeshow Goodness-of-Fit
Hipotesis:
\(H_0\): Model memiliki kecocokan yang
baik; tidak ada perbedaan bermakna antara frekuensi aktual dan frekuensi
harapan.
\(H_1\): Model tidak memiliki kecocokan
yang baik; terdapat perbedaan bermakna antara frekuensi aktual dan
frekuensi harapan.
binary_train_prob <- predict(binary_best_fit, newdata=binary_train, type="response")
binary_hl <- hosmer_lemeshow_table(actual=binary_train$food_crisis, predicted=binary_train_prob, g=10)
binary_hl_test <- binary_hl$test |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil)
safe_table(binary_hl$detail, "Detail pengelompokan Hosmer-Lemeshow model terbaik logistik biner", digits=3, full_width=TRUE)safe_table(binary_hl_test, "Uji Hosmer-Lemeshow Goodness-of-Fit model terbaik logistik biner", digits=4, full_width=TRUE)3.6.4 Linearitas Log-Odds (Box-Tidwell)
Hipotesis:
\(H_0\): Hubungan prediktor numerik
dengan log-odds bersifat linear.
\(H_1\): Hubungan prediktor numerik
dengan log-odds tidak linear.
binary_box_tidwell <- box_tidwell_table(binary_train, response="food_crisis", predictors=binary_selection$best$vars, family=binomial(link="logit")) |>
mutate(p_value=p_value_tampil) |>
select(-p_value_tampil)
safe_table(binary_box_tidwell, "Uji linearitas log-odds Box-Tidwell model terbaik logistik biner", digits=4, full_width=TRUE)3.6.5 Nagelkerke \(R^2\)
Nagelkerke \(R^2\) bukan uji hipotesis, sehingga tidak memiliki \(H_0\) dan \(H_1\). Ukuran ini digunakan sebagai indeks kekuatan penjelasan model. Nilai yang lebih besar menunjukkan model memberikan peningkatan kecocokan yang lebih baik dibanding model kosong.
binary_nagelkerke_table <- tibble(
ukuran="Nagelkerke R-square",
nilai=nagelkerke_r2(binary_best_fit, binary_best_null),
keterangan="Semakin besar nilainya, semakin besar peningkatan kecocokan model dibanding model kosong."
)
safe_table(binary_nagelkerke_table, "Nagelkerke R-square model terbaik logistik biner", digits=4, full_width=TRUE)3.6.6 Diagnostik Pengaruh Observasi
Diagnostik pengaruh observasi bukan uji hipotesis formal. Pemeriksaan ini digunakan untuk melihat apakah ada observasi yang terlalu memengaruhi estimasi model. Kriteria praktis: Cook’s distance > 1 perlu diperiksa, dan leverage tinggi dilihat dengan batas \(2p/n\).
binary_influence_table <- influence_check_table(binary_best_fit)
safe_table(binary_influence_table, "Diagnostik pengaruh observasi model terbaik logistik biner", digits=4, full_width=TRUE)3.7 Prediksi dengan model terbaik
binary_prob <- predict(binary_best_fit, newdata=binary_test, type="response")
binary_pred <- ifelse(binary_prob >= .50, 1, 0)
binary_cm <- table(Aktual=binary_test$food_crisis, Prediksi=binary_pred)
binary_metric_table <- binary_metrics(binary_test$food_crisis, binary_pred)
binary_cm_table <- cm_labeled(binary_cm)
names(binary_cm_table) <- c("Aktual \\ Prediksi", "Prediksi: non_food_crisis (0)", "Prediksi: food_crisis (1)")
binary_cm_table[["Aktual \\ Prediksi"]] <- c("Aktual: non_food_crisis (0)", "Aktual: food_crisis (1)")
safe_table(binary_cm_table, "Confusion matrix logistik biner model terbaik dengan keterangan aktual dan prediksi")safe_table(binary_metric_table, "Metrik evaluasi testing model terbaik logistik biner pada threshold 0,50", digits=3)binary_roc <- pROC::roc(response=binary_test$food_crisis, predictor=binary_prob, quiet=TRUE)
binary_auc <- as.numeric(pROC::auc(binary_roc))
roc_df <- tibble(specificity=binary_roc$specificities, sensitivity=binary_roc$sensitivities) |> mutate(false_positive_rate=1-specificity)
roc_df |> ggplot(aes(x=false_positive_rate,y=sensitivity))+geom_line(linewidth=1.15,color="#6d28d9")+geom_abline(slope=1,intercept=0,linetype="dashed",color="#e11d48")+coord_equal()+labs(title="Kurva ROC Model Terbaik Logistik Biner", subtitle=paste0("AUC = ",round(binary_auc,3)), x="False positive rate", y="True positive rate / sensitivity")+theme_report()
Gambar 2. Kurva ROC model terbaik logistik biner.
Pada data testing, model terbaik menghasilkan accuracy sebesar 0.847,
sensitivity sebesar 0.74, specificity sebesar 0.899, dan AUC sebesar
0.921. Sensitivity penting karena menunjukkan kemampuan model mengenali
observasi yang benar-benar mengalami food_crisis.
3.8 Korelasi prediktor
binary_corr_top <- correlation_pairs(binary_train |> select(all_of(binary_x)), top_n=10)
binary_corr_high <- correlation_pairs(binary_train |> select(all_of(binary_x)), top_n=100, threshold=.80)
safe_table(binary_corr_top, "Sepuluh pasangan prediktor biner dengan korelasi absolut terbesar", digits=3)safe_table(binary_corr_high, "Pasangan prediktor biner dengan korelasi tinggi atau |r| >= 0,80", digits=3)4 Bagian II — Regresi Logistik Multinomial
4.1 Konsep model
Regresi logistik multinomial digunakan ketika respons bersifat nominal dan memiliki lebih dari dua kategori. Dengan kategori referensi \(Y=1\), model untuk kategori \(j\) adalah:
\[ \log\left(\frac{P(Y_i=j)}{P(Y_i=1)}\right)=\alpha_j+\beta_{j1}X_{1i}+\cdots+\beta_{jk}X_{ki} \]
Ukuran efek yang digunakan adalah relative risk ratio (RRR), yaitu \(\exp(\beta)\).
4.2 Import, struktur data, dan distribusi respons
multinom_data <- readxl::read_excel(resolve_path("multinomial_model_ready_sample.xlsx")) |>
clean_model_data() |>
mutate(choice_code=factor(choice_code, levels=sort(unique(choice_code)), labels=c("protect_environment","economic_growth_jobs","other_answer")[seq_along(sort(unique(choice_code)))]))
multinom_dict <- make_dictionary(multinom_data)
multinom_summary <- multinom_data |> count(choice_code, name="n") |> mutate(proporsi=n/sum(n))
safe_table(multinom_dict, "Kamus data multinomial setelah seleksi variabel analisis")safe_table(multinom_summary, "Distribusi respons multinomial choice_code", digits=3)multinom_summary |> ggplot(aes(x=reorder(choice_code,n), y=n, fill=choice_code))+geom_col(width=.70)+geom_text(aes(label=n), hjust=-.15, fontface="bold")+coord_flip()+scale_fill_manual(values=report_palette)+labs(title="Distribusi Environment Priority", subtitle="Respons multinomial dengan kategori nominal.", x="Kategori respons", y="Jumlah observasi", fill="Kategori")+theme_report()
Gambar 3. Distribusi kategori respons multinomial.
Respons choice_code bersifat nominal, sehingga tidak
tepat diperlakukan sebagai urutan rendah-tinggi. Ketimpangan jumlah
observasi antarkategori juga perlu dibaca karena kategori dengan
frekuensi kecil dapat lebih sulit diprediksi secara stabil.
4.3 Statistik Deskriptif dan Eksplorasi Data
multinom_descriptive <- numeric_descriptive(multinom_data, response="choice_code")
multinom_x_preview <- model_predictors(multinom_data, response="choice_code")
multinom_group_summary <- multinom_data |>
select(choice_code, all_of(multinom_x_preview)) |>
group_by(choice_code) |>
summarise(across(everything(), ~mean(.x, na.rm=TRUE)), .groups="drop") |>
pivot_longer(-choice_code, names_to="variabel", values_to="mean") |>
pivot_wider(names_from=choice_code, values_from=mean)
safe_table(multinom_descriptive, "Statistik deskriptif prediktor numerik pada data multinomial", digits=3, full_width=TRUE)safe_table(multinom_group_summary, "Rata-rata prediktor menurut kategori choice_code", digits=3, full_width=TRUE)multinom_plot_vars <- multinom_x_preview[seq_len(min(4, length(multinom_x_preview)))]
multinom_data |>
select(choice_code, all_of(multinom_plot_vars)) |>
pivot_longer(-choice_code, names_to="variabel", values_to="nilai") |>
ggplot(aes(x=choice_code, y=nilai, fill=choice_code))+
geom_boxplot(alpha=.75, outlier.alpha=.25)+
facet_wrap(~variabel, scales="free_y")+
scale_fill_manual(values=report_palette)+
labs(title="Eksplorasi Prediktor terhadap Environment Priority", x="Kategori respons", y="Nilai prediktor", fill="Kategori")+
theme_report()
Gambar 4. Eksplorasi prediktor numerik utama menurut kategori choice_code.
4.4 Pembagian data training dan testing
multinom_split <- split_train_test(multinom_data, strata="choice_code", prop=.80)
multinom_train <- multinom_split$train; multinom_test <- multinom_split$test
multinom_split_table <- bind_rows(multinom_train |> count(choice_code, name="n") |> mutate(subset="Training"), multinom_test |> count(choice_code, name="n") |> mutate(subset="Testing")) |> group_by(subset) |> mutate(proporsi=n/sum(n)) |> ungroup() |> select(subset, choice_code, n, proporsi)
safe_table(multinom_split_table, "Komposisi respons multinomial pada data training dan testing", digits=3)4.5 Kandidat model dan pemilihan model terbaik
multinom_x <- model_predictors(multinom_train, response="choice_code")
multinom_selection <- select_best_multinom(multinom_train, response="choice_code", predictors=multinom_x)
multinom_best_fit <- multinom_selection$best$fit
multinom_best_null <- nnet::multinom(choice_code ~ 1, data=multinom_train, trace=FALSE, Hess=TRUE, model=TRUE)
multinom_best_wald <- multinom_selection$best$wald
multinom_best_g2 <- multinom_selection$best$g2
multinom_selected_vars <- multinom_best_wald |> group_by(variabel) |> summarise(min_p=min(p_value, na.rm=TRUE), .groups="drop") |> filter(min_p < 0.05) |> pull(variabel)
multinom_model_selection_table <- multinom_selection$summary |>
mutate(p_value_G2=p_value_G2_tampil) |>
select(model_ke, formula, jumlah_prediktor, jumlah_variabel_signifikan, AIC, G2, p_value_G2)
multinom_fit_table <- tibble(keterangan=c("Jumlah observasi training","Jumlah prediktor model terbaik","AIC model terbaik","Log-likelihood model terbaik","Log-likelihood model kosong","McFadden pseudo R-square"), nilai=c(nrow(multinom_train), length(multinom_selection$best$vars), AIC(multinom_best_fit), as.numeric(logLik(multinom_best_fit)), as.numeric(logLik(multinom_best_null)), mcfadden_r2(multinom_best_fit,multinom_best_null)))
multinom_best_wald_show <- multinom_best_wald |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil)
multinom_best_g2_show <- multinom_best_g2 |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil)
safe_table(multinom_model_selection_table, "Perbandingan kandidat model multinomial berdasarkan jumlah variabel signifikan dan AIC", digits=4, full_width=TRUE)safe_table(multinom_fit_table, "Ringkasan kecocokan model terbaik multinomial", digits=4)safe_table(multinom_best_g2_show, "Uji simultan G² model terbaik multinomial", digits=4)safe_table(multinom_best_wald_show |> slice_head(n=30), "Uji Wald parsial dan RRR model terbaik multinomial", digits=4, full_width=TRUE)Model terbaik multinomial adalah kandidat ke-9. Model ini dipilih
karena mempertahankan 3 variabel signifikan, kemudian dibandingkan
berdasarkan AIC. Uji \(G^2\) model
terbaik memiliki p-value < 0.001, sehingga model dengan prediktor
lebih informatif daripada model kosong. Variabel X yang mendukung
perbedaan pilihan choice_code adalah: education_code,
employment_label_other, employment_label_student. Karena respons
multinomial memiliki beberapa kategori, satu variabel dapat signifikan
pada salah satu perbandingan kategori meskipun tidak signifikan pada
perbandingan lain.
4.6 Uji Asumsi Model Terbaik
4.6.1 Asumsi 1 — Respons Bersifat Nominal
Hipotesis/Kriteria:
\(H_0\): Variabel respons bersifat
nominal dan kategori respons saling eksklusif.
\(H_1\): Variabel respons tidak
bersifat nominal atau kategori tidak saling eksklusif.
Kriteria: regresi multinomial sesuai jika respons memiliki lebih dari
dua kategori nominal tanpa urutan alami.
multinom_nominal_table <- multinom_train |>
count(choice_code, name="n") |>
mutate(proporsi=n/sum(n), skala_respons="Nominal")
safe_table(multinom_nominal_table, "Pemeriksaan skala respons nominal pada model multinomial", digits=3, full_width=TRUE)4.6.2 Asumsi 2 — Observasi Independen
Hipotesis/Kriteria:
\(H_0\): Observasi saling
independen.
\(H_1\): Observasi tidak saling
independen atau terdapat struktur pengelompokan/pengulangan.
Kriteria: setiap baris data harus merepresentasikan unit observasi
berbeda dan tidak berulang tanpa koreksi.
multinom_independent_table <- tibble(
asumsi="Observasi independen",
pemeriksaan="Setiap baris data diperlakukan sebagai satu unit observasi yang berbeda.",
hasil="Tidak ada struktur pengulangan responden yang digunakan dalam model-ready.",
catatan="Apabila terdapat clustering survei atau panel, koreksi standard error dapat dipertimbangkan."
)
safe_table(multinom_independent_table, "Pemeriksaan asumsi observasi independen pada model multinomial", full_width=TRUE)4.6.3 Asumsi 3 — Tidak Ada Multikolinearitas Berat (VIF)
Hipotesis/Kriteria:
\(H_0\): Tidak terdapat
multikolinearitas berat antarprediktor.
\(H_1\): Terdapat multikolinearitas
berat antarprediktor.
Kriteria praktis: VIF < 5 aman; VIF 5–10 perlu perhatian; VIF > 10
menunjukkan multikolinearitas berat.
multinom_vif_table <- vif_table(multinom_train, multinom_selection$best$vars)
safe_table(multinom_vif_table, "Uji multikolinearitas model terbaik multinomial menggunakan VIF", digits=3, full_width=TRUE)4.6.4 Asumsi 4 — Tidak Ada Perfect Separation
Hipotesis/Kriteria:
\(H_0\): Tidak terdapat perfect
separation atau prediksi yang terlalu ekstrem.
\(H_1\): Terdapat indikasi perfect
separation.
Kriteria praktis: proporsi probabilitas maksimum > 0,99 yang sangat
tinggi dapat menjadi sinyal separasi atau prediksi ekstrem.
multinom_prob_train <- predict(multinom_best_fit, newdata=multinom_train, type="probs")
if(is.null(dim(multinom_prob_train))){multinom_prob_train <- cbind(multinom_prob_train, 1-multinom_prob_train)}
multinom_separation_table <- tibble(
diagnostik=c("Proporsi probabilitas maksimum > 0,99", "Probabilitas maksimum terbesar", "Probabilitas minimum terkecil"),
nilai=c(mean(apply(multinom_prob_train, 1, max) > .99), max(multinom_prob_train), min(multinom_prob_train)),
keterangan=c(
"Nilai tinggi dapat mengindikasikan prediksi ekstrem.",
"Jika mendekati 1 untuk banyak observasi, perlu waspada separasi.",
"Jika mendekati 0 untuk banyak observasi, perlu waspada separasi."
)
)
safe_table(multinom_separation_table, "Diagnostik perfect separation model terbaik multinomial", digits=4, full_width=TRUE)4.6.5 Asumsi 5 — Linearitas pada Skala Logit (Box-Tidwell)
Hipotesis:
\(H_0\): Hubungan prediktor numerik
dengan skala logit bersifat linear.
\(H_1\): Hubungan prediktor numerik
dengan skala logit tidak linear.
Karena model multinomial memiliki beberapa kategori, pemeriksaan
dilakukan secara one-vs-rest.
multinom_bt_table <- purrr::map_dfr(levels(multinom_train$choice_code), function(k){
tmp <- multinom_train |> mutate(.y_bin=as.integer(choice_code == k))
box_tidwell_table(tmp, response=".y_bin", predictors=multinom_selection$best$vars, family=binomial(link="logit")) |>
mutate(kategori_one_vs_rest=k)
}) |>
mutate(p_value=p_value_tampil) |>
select(kategori_one_vs_rest, everything(), -p_value_tampil)
safe_table(multinom_bt_table, "Diagnostik Box-Tidwell one-vs-rest untuk model multinomial", digits=4, full_width=TRUE)4.6.6 Asumsi 6 — Independence of Irrelevant Alternatives (IIA)
Hipotesis/Kriteria:
\(H_0\): Perbandingan odds antara dua
kategori tidak dipengaruhi secara tidak wajar oleh keberadaan kategori
lain.
\(H_1\): Perbandingan odds antara dua
kategori dipengaruhi oleh keberadaan kategori lain.
Pada laporan ini, IIA dievaluasi secara substantif karena uji formal
memerlukan struktur data dan paket khusus.
multinom_iia_table <- tibble(
asumsi="Independence of Irrelevant Alternatives (IIA)",
pemeriksaan="Evaluasi substantif berdasarkan struktur pilihan kategori.",
hasil="Kategori respons diperlakukan sebagai alternatif nominal yang saling eksklusif.",
catatan="Uji formal IIA seperti Hausman-McFadden memerlukan struktur data khusus; pada laporan ini IIA dibahas sebagai asumsi substantif model multinomial."
)
safe_table(multinom_iia_table, "Pemeriksaan asumsi IIA pada model multinomial", full_width=TRUE)4.7 Prediksi dengan model terbaik
multinom_pred <- predict(multinom_best_fit, newdata=multinom_test, type="class")
multinom_prob <- predict(multinom_best_fit, newdata=multinom_test, type="probs")
multinom_metrics <- multi_metrics(multinom_test$choice_code, multinom_pred)
safe_table(cm_labeled(multinom_metrics$confusion_matrix), "Confusion matrix model terbaik multinomial dengan keterangan aktual dan prediksi")safe_table(multinom_metrics$accuracy_table, "Metrik agregat testing model terbaik multinomial", digits=3)safe_table(multinom_metrics$class_metrics, "Metrik evaluasi per kelas model terbaik multinomial", digits=3)multinom_prob_df <- as.data.frame(multinom_prob) |> mutate(observasi=row_number()) |> pivot_longer(-observasi, names_to="kategori", values_to="probabilitas")
multinom_prob_df |> group_by(kategori) |> summarise(rata_rata_probabilitas=mean(probabilitas), .groups="drop") |> ggplot(aes(x=reorder(kategori,rata_rata_probabilitas), y=rata_rata_probabilitas, fill=kategori))+geom_col(width=.70)+geom_text(aes(label=round(rata_rata_probabilitas,3)), hjust=-.12, fontface="bold")+coord_flip()+scale_fill_manual(values=report_palette)+scale_y_continuous(labels=percent_format(accuracy=1))+labs(title="Rata-rata Probabilitas Prediksi Model Terbaik Multinomial", subtitle="Diringkas dari data testing.", x="Kategori", y="Rata-rata probabilitas prediksi", fill="Kategori")+theme_report()
Gambar 4. Rata-rata probabilitas prediksi model terbaik multinomial.
Prediksi dilakukan dengan memilih kategori yang memiliki probabilitas terbesar. Accuracy model terbaik pada data testing adalah 0.572. Karena respons memiliki beberapa kategori, evaluasi model juga perlu membaca sensitivity, specificity, precision, dan F1-score pada setiap kelas.
4.8 Korelasi prediktor multinomial
multinom_corr_top <- correlation_pairs(multinom_train |> select(all_of(multinom_x)), top_n=10)
safe_table(multinom_corr_top, "Sepuluh pasangan prediktor multinomial dengan korelasi absolut terbesar", digits=3)5 Bagian III — Regresi Logistik Ordinal
5.1 Konsep model
Regresi logistik ordinal digunakan ketika respons memiliki urutan kategori. Pada model cumulative logit:
\[ \log\left(\frac{P(Y_i\le j)}{P(Y_i>j)}\right)=\theta_j-\left(\beta_1X_{1i}+\cdots+\beta_kX_{ki}\right) \]
Model ini mempertahankan informasi urutan tanpa menganggap jarak antar kategori harus sama.
5.2 Import, struktur data, dan distribusi respons
ordinal_data <- readxl::read_excel(resolve_path("ordinal_model_ready_sample.xlsx")) |>
clean_model_data() |>
mutate(gov_confidence_ord=ordered(gov_confidence_ord, levels=sort(unique(gov_confidence_ord)), labels=c("none_at_all","not_very_much","quite_a_lot","a_great_deal")[seq_along(sort(unique(gov_confidence_ord)))]))
ordinal_dict <- make_dictionary(ordinal_data)
ordinal_summary <- ordinal_data |> count(gov_confidence_ord, name="n") |> mutate(proporsi=n/sum(n))
safe_table(ordinal_dict, "Kamus data ordinal setelah seleksi variabel analisis")safe_table(ordinal_summary, "Distribusi respons ordinal gov_confidence_ord", digits=3)ordinal_summary |> ggplot(aes(x=gov_confidence_ord, y=n, fill=gov_confidence_ord))+geom_col(width=.70)+geom_text(aes(label=n), vjust=-.45, fontface="bold")+scale_fill_manual(values=report_palette)+labs(title="Distribusi Government Confidence", subtitle="Kategori respons memiliki urutan alami.", x="Tingkat kepercayaan terhadap pemerintah", y="Jumlah observasi", fill="Kategori")+theme_report()
Gambar 5. Distribusi kategori respons ordinal.
Variabel gov_confidence_ord memiliki urutan substantif
dari kepercayaan rendah ke tinggi. Karena jarak antar kategori tidak
harus sama, regresi logistik ordinal lebih tepat dibanding model linear
biasa.
5.3 Statistik Deskriptif dan Eksplorasi Data
ordinal_descriptive <- numeric_descriptive(ordinal_data, response="gov_confidence_ord")
ordinal_x_preview <- model_predictors(ordinal_data, response="gov_confidence_ord")
ordinal_group_summary <- ordinal_data |>
select(gov_confidence_ord, all_of(ordinal_x_preview)) |>
group_by(gov_confidence_ord) |>
summarise(across(everything(), ~mean(.x, na.rm=TRUE)), .groups="drop") |>
pivot_longer(-gov_confidence_ord, names_to="variabel", values_to="mean") |>
pivot_wider(names_from=gov_confidence_ord, values_from=mean)
safe_table(ordinal_descriptive, "Statistik deskriptif prediktor numerik pada data ordinal", digits=3, full_width=TRUE)safe_table(ordinal_group_summary, "Rata-rata prediktor menurut tingkat gov_confidence_ord", digits=3, full_width=TRUE)ordinal_plot_vars <- ordinal_x_preview[seq_len(min(4, length(ordinal_x_preview)))]
ordinal_data |>
select(gov_confidence_ord, all_of(ordinal_plot_vars)) |>
pivot_longer(-gov_confidence_ord, names_to="variabel", values_to="nilai") |>
ggplot(aes(x=gov_confidence_ord, y=nilai, fill=gov_confidence_ord))+
geom_boxplot(alpha=.75, outlier.alpha=.25)+
facet_wrap(~variabel, scales="free_y")+
scale_fill_manual(values=report_palette)+
labs(title="Eksplorasi Prediktor terhadap Government Confidence", x="Tingkat respons ordinal", y="Nilai prediktor", fill="Kategori")+
theme_report()
Gambar 6. Eksplorasi prediktor numerik utama menurut tingkat gov_confidence_ord.
5.4 Pembagian data training dan testing
ordinal_split <- split_train_test(ordinal_data, strata="gov_confidence_ord", prop=.80)
ordinal_train <- ordinal_split$train; ordinal_test <- ordinal_split$test
ordinal_split_table <- bind_rows(ordinal_train |> count(gov_confidence_ord, name="n") |> mutate(subset="Training"), ordinal_test |> count(gov_confidence_ord, name="n") |> mutate(subset="Testing")) |> group_by(subset) |> mutate(proporsi=n/sum(n)) |> ungroup() |> select(subset, gov_confidence_ord, n, proporsi)
safe_table(ordinal_split_table, "Komposisi respons ordinal pada data training dan testing", digits=3)5.5 Kandidat model dan pemilihan model terbaik
ordinal_x <- model_predictors(ordinal_train, response="gov_confidence_ord")
ordinal_selection <- select_best_polr(ordinal_train, response="gov_confidence_ord", predictors=ordinal_x)
ordinal_best_fit <- ordinal_selection$best$fit
ordinal_best_null <- MASS::polr(gov_confidence_ord ~ 1, data=ordinal_train, Hess=TRUE, method="logistic")
ordinal_best_wald <- ordinal_selection$best$wald
ordinal_best_g2 <- ordinal_selection$best$g2
ordinal_selected_vars <- ordinal_best_wald |> filter(p_value < 0.05) |> pull(variabel)
ordinal_model_selection_table <- ordinal_selection$summary |>
mutate(p_value_G2=p_value_G2_tampil) |>
select(model_ke, formula, jumlah_prediktor, jumlah_variabel_signifikan, AIC, G2, p_value_G2)
ordinal_fit_table <- tibble(keterangan=c("Jumlah observasi training","Jumlah prediktor model terbaik","AIC model terbaik","Log-likelihood model terbaik","Log-likelihood model kosong","McFadden pseudo R-square"), nilai=c(nrow(ordinal_train), length(ordinal_selection$best$vars), AIC(ordinal_best_fit), as.numeric(logLik(ordinal_best_fit)), as.numeric(logLik(ordinal_best_null)), mcfadden_r2(ordinal_best_fit,ordinal_best_null)))
ordinal_best_wald_show <- ordinal_best_wald |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil)
ordinal_best_g2_show <- ordinal_best_g2 |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil)
safe_table(ordinal_model_selection_table, "Perbandingan kandidat model ordinal berdasarkan jumlah variabel signifikan dan AIC", digits=4, full_width=TRUE)safe_table(ordinal_fit_table, "Ringkasan kecocokan model terbaik ordinal", digits=4)safe_table(ordinal_best_g2_show, "Uji simultan G² model terbaik ordinal", digits=4)safe_table(ordinal_best_wald_show |> slice_head(n=30), "Uji Wald parsial dan odds ratio model terbaik ordinal", digits=4, full_width=TRUE)Model terbaik ordinal adalah kandidat ke-9, dengan 3 variabel
signifikan dan AIC sebesar 1.5160099^{4}. Uji \(G^2\) model terbaik memiliki p-value <
0.001, sehingga model terpilih memberikan peningkatan informasi
dibanding model kosong. Variabel X yang mendukung variasi tingkat
gov_confidence_ord adalah: education_code, income_scale,
employment_label_unemployed. Dengan demikian, prediksi ordinal berikut
dilakukan menggunakan model terbaik tersebut.
5.6 Uji Asumsi Model Terbaik
5.6.1 Asumsi 1 — Asumsi Proportional Odds
Hipotesis/Kriteria:
\(H_0\): Efek prediktor relatif sama
pada setiap batas kumulatif kategori respons.
\(H_1\): Efek prediktor berbeda pada
beberapa batas kumulatif kategori respons.
Kriteria praktis: rentang koefisien logit kumulatif yang kecil
menunjukkan efek prediktor relatif stabil; rentang yang besar
mengindikasikan potensi pelanggaran proportional odds.
ordinal_parallel <- ordinal_parallel_diagnostic(
data=ordinal_train,
response="gov_confidence_ord",
predictors=ordinal_selection$best$vars
)
safe_table(ordinal_parallel$ringkas, "Ringkasan diagnostik proportional odds berdasarkan rentang koefisien logit kumulatif", digits=4, full_width=TRUE)safe_table(ordinal_parallel$detail |> slice_head(n=30) |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil), "Detail uji Wald pada logit kumulatif sebagai diagnostik proportional odds", digits=4, full_width=TRUE)5.6.2 Asumsi 2 — Tidak Ada Multikolinearitas Berat
Hipotesis/Kriteria:
\(H_0\): Tidak terdapat
multikolinearitas berat antarprediktor.
\(H_1\): Terdapat multikolinearitas
berat antarprediktor.
Kriteria praktis: VIF < 5 aman; VIF 5–10 perlu perhatian; VIF > 10
menunjukkan multikolinearitas berat.
ordinal_vif_table <- vif_table(ordinal_train, ordinal_selection$best$vars)
ordinal_corr_best_table <- selected_correlation_check(ordinal_train, ordinal_selection$best$vars)
safe_table(ordinal_vif_table, "Uji multikolinearitas model terbaik ordinal menggunakan VIF", digits=3, full_width=TRUE)safe_table(ordinal_corr_best_table, "Pemeriksaan korelasi tinggi antarprediktor model terbaik ordinal", digits=3, full_width=TRUE)5.7 Prediksi dengan model terbaik
ordinal_pred <- predict(ordinal_best_fit, newdata=ordinal_test, type="class")
ordinal_prob <- predict(ordinal_best_fit, newdata=ordinal_test, type="probs")
ordinal_metrics <- multi_metrics(ordinal_test$gov_confidence_ord, ordinal_pred)
safe_table(cm_labeled(ordinal_metrics$confusion_matrix), "Confusion matrix model terbaik ordinal dengan keterangan aktual dan prediksi")safe_table(ordinal_metrics$accuracy_table, "Metrik agregat testing model terbaik ordinal", digits=3)safe_table(ordinal_metrics$class_metrics, "Metrik evaluasi per kelas model terbaik ordinal", digits=3)ordinal_prob_df <- as.data.frame(ordinal_prob) |> mutate(observasi=row_number()) |> pivot_longer(-observasi, names_to="kategori", values_to="probabilitas")
ordinal_prob_df |> group_by(kategori) |> summarise(rata_rata_probabilitas=mean(probabilitas), .groups="drop") |> ggplot(aes(x=kategori, y=rata_rata_probabilitas, fill=kategori))+geom_col(width=.70)+geom_text(aes(label=round(rata_rata_probabilitas,3)), vjust=-.45, fontface="bold")+scale_fill_manual(values=report_palette)+scale_y_continuous(labels=percent_format(accuracy=1))+labs(title="Rata-rata Probabilitas Prediksi Model Terbaik Ordinal", subtitle="Diringkas dari data testing.", x="Kategori ordinal", y="Rata-rata probabilitas prediksi", fill="Kategori")+theme_report()
Gambar 6. Rata-rata probabilitas prediksi model terbaik ordinal.
Accuracy testing model terbaik ordinal adalah 0.335. Pada respons ordinal, pola kesalahan juga perlu dibaca karena kesalahan ke kategori yang berdekatan lebih ringan secara substantif dibanding kesalahan ke kategori yang jauh.
5.8 Korelasi prediktor ordinal
ordinal_corr_top <- correlation_pairs(ordinal_train |> select(all_of(ordinal_x)), top_n=10)
safe_table(ordinal_corr_top, "Sepuluh pasangan prediktor ordinal dengan korelasi absolut terbesar", digits=3)6 Bagian IV — Regresi Poisson Tahunan
6.1 Konsep model
Regresi Poisson digunakan untuk respons berupa hitungan non-negatif. Model menggunakan link log:
\[ \log(\mu_i)=\beta_0+\beta_1X_{1i}+\cdots+\beta_kX_{ki} \]
Ukuran efek yang digunakan adalah incidence rate ratio (IRR), yaitu \(\exp(\beta)\).
6.2 Import, struktur data, dan distribusi respons
poisson_data <- readxl::read_excel(resolve_path("poisson_model_ready.xlsx")) |>
clean_model_data()
poisson_overdispersion_file <- readxl::read_excel(resolve_path("poisson_overdispersion_check.xlsx"))
poisson_duplicate_check <- tibble(variabel_1="lag_conflict_count", variabel_2="lag_actor_count", identik=if(all(c("lag_conflict_count","lag_actor_count") %in% names(poisson_data))) isTRUE(all.equal(poisson_data$lag_conflict_count, poisson_data$lag_actor_count)) else NA, keputusan="lag_actor_count tidak digunakan pada model akhir")
poisson_data <- poisson_data |> select(-any_of("lag_actor_count"))
poisson_dict <- make_dictionary(poisson_data)
poisson_summary <- tibble(n_observasi=nrow(poisson_data), min_case_count=min(poisson_data$case_count,na.rm=TRUE), mean_case_count=mean(poisson_data$case_count,na.rm=TRUE), variance_case_count=var(poisson_data$case_count,na.rm=TRUE), median_case_count=median(poisson_data$case_count,na.rm=TRUE), max_case_count=max(poisson_data$case_count,na.rm=TRUE), proporsi_nol=mean(poisson_data$case_count==0,na.rm=TRUE))
safe_table(poisson_dict, "Kamus data Poisson tahunan setelah seleksi variabel analisis")safe_table(poisson_duplicate_check, "Cek duplikasi prediktor lag pada data Poisson tahunan")safe_table(poisson_overdispersion_file, "Ringkasan awal cek overdispersion Poisson tahunan")safe_table(poisson_summary, "Ringkasan deskriptif respons case_count pada data Poisson tahunan", digits=3)poisson_data |> ggplot(aes(x=case_count))+geom_histogram(binwidth=1,boundary=-.5,fill="#0891b2",color="white")+scale_x_continuous(breaks=scales::pretty_breaks())+labs(title="Distribusi Violence Case Count Tahunan", subtitle="Respons berupa data hitungan non-negatif.", x="Jumlah kejadian kekerasan", y="Frekuensi")+theme_report()
Gambar 7. Distribusi respons case_count pada data Poisson tahunan.
case_count memenuhi bentuk dasar respons Poisson karena
berupa hitungan non-negatif. Rata-rata dan varians respons perlu
diperhatikan karena model Poisson klasik mengasumsikan equidispersion,
yaitu rata-rata dan varians kondisional relatif sama.
6.3 Statistik Deskriptif dan Eksplorasi Data
poisson_descriptive <- numeric_descriptive(poisson_data, response="case_count")
poisson_x_preview <- model_predictors(poisson_data, response="case_count")
poisson_response_exploration <- tibble(
ukuran=c("Mean case_count", "Variance case_count", "Variance / Mean", "Proporsi nilai nol"),
nilai=c(
mean(poisson_data$case_count, na.rm=TRUE),
var(poisson_data$case_count, na.rm=TRUE),
var(poisson_data$case_count, na.rm=TRUE)/mean(poisson_data$case_count, na.rm=TRUE),
mean(poisson_data$case_count == 0, na.rm=TRUE)
)
)
safe_table(poisson_descriptive, "Statistik deskriptif prediktor numerik pada data Poisson tahunan", digits=3, full_width=TRUE)safe_table(poisson_response_exploration, "Eksplorasi awal respons case_count untuk indikasi overdispersi dan zero inflation", digits=4, full_width=TRUE)poisson_plot_vars <- poisson_x_preview[seq_len(min(4, length(poisson_x_preview)))]
poisson_data |>
select(case_count, all_of(poisson_plot_vars)) |>
pivot_longer(-case_count, names_to="variabel", values_to="nilai") |>
ggplot(aes(x=nilai, y=case_count))+
geom_point(alpha=.35, color="#172033")+
geom_smooth(method="glm", method.args=list(family=poisson(link="log")), se=TRUE, color="#6d28d9", fill="#6d28d9", alpha=.18)+
facet_wrap(~variabel, scales="free_x")+
labs(title="Eksplorasi Hubungan Prediktor dengan Violence Case Count", x="Nilai prediktor", y="case_count")+
theme_report()
Gambar 8. Hubungan prediktor numerik utama dengan case_count.
6.4 Pembagian data training dan testing
poisson_split <- split_train_test(poisson_data, prop=.80)
poisson_train <- poisson_split$train; poisson_test <- poisson_split$test
poisson_split_table <- tibble(subset=c("Training","Testing"), n=c(nrow(poisson_train),nrow(poisson_test)), mean_case_count=c(mean(poisson_train$case_count),mean(poisson_test$case_count)), proporsi_nol=c(mean(poisson_train$case_count==0),mean(poisson_test$case_count==0)))
safe_table(poisson_split_table, "Komposisi training dan testing pada data Poisson tahunan", digits=3)6.5 Kandidat model dan pemilihan model terbaik
poisson_x <- model_predictors(poisson_train, response="case_count")
poisson_selection <- select_best_glm(poisson_train, response="case_count", predictors=poisson_x, family=poisson(link="log"), effect_name="incidence_rate_ratio")
poisson_best_fit <- poisson_selection$best$fit
poisson_best_null <- glm(case_count ~ 1, data=poisson_train, family=poisson(link="log"))
poisson_best_wald <- poisson_selection$best$wald
poisson_best_g2 <- poisson_selection$best$g2
poisson_selected_vars <- poisson_best_wald |> filter(p_value < 0.05) |> pull(variabel)
poisson_model_selection_table <- poisson_selection$summary |>
mutate(p_value_G2=p_value_G2_tampil) |>
select(model_ke, formula, jumlah_prediktor, jumlah_variabel_signifikan, AIC, G2, p_value_G2)
poisson_fit_table <- tibble(keterangan=c("Jumlah observasi training","Jumlah prediktor model terbaik","Null deviance model terbaik","Residual deviance model terbaik","AIC model terbaik","Log-likelihood model terbaik","Log-likelihood model kosong","McFadden pseudo R-square"), nilai=c(nrow(poisson_train), length(poisson_selection$best$vars), poisson_best_fit$null.deviance, poisson_best_fit$deviance, AIC(poisson_best_fit), as.numeric(logLik(poisson_best_fit)), as.numeric(logLik(poisson_best_null)), mcfadden_r2(poisson_best_fit,poisson_best_null)))
poisson_best_wald_show <- poisson_best_wald |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil)
poisson_best_g2_show <- poisson_best_g2 |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil)
safe_table(poisson_model_selection_table, "Perbandingan kandidat model Poisson berdasarkan jumlah variabel signifikan dan AIC", digits=4, full_width=TRUE)safe_table(poisson_fit_table, "Ringkasan kecocokan model terbaik Poisson tahunan", digits=4)safe_table(poisson_best_g2_show, "Uji simultan G² model terbaik Poisson tahunan", digits=4)safe_table(poisson_best_wald_show, "Uji Wald parsial dan IRR model terbaik Poisson tahunan", digits=4, full_width=TRUE)Model terbaik Poisson adalah kandidat ke-4, dengan 5 variabel
signifikan dan AIC sebesar 2628.673. Uji \(G^2\) model terbaik memiliki p-value <
0.001, sehingga prediktor pada model terbaik secara simultan
berkontribusi terhadap expected count case_count. Variabel
X yang mendukung case_count adalah: log1p_lag_case_count,
region_label_Europe, lag_conflict_count, log1p_lag_total_fatalities,
region_label_Middle_East. Prediksi berikutnya dilakukan menggunakan
model terbaik tersebut.
6.6 Uji Asumsi Model Terbaik
6.6.1 Asumsi 1 — Overdispersi (Pearson Dispersion)
Hipotesis:
\(H_0\): Tidak terdapat overdispersi;
varians kondisional relatif sesuai dengan mean kondisional.
\(H_1\): Terdapat overdispersi; varians
kondisional lebih besar daripada mean kondisional.
Kriteria praktis: dispersion Pearson mendekati 1 menunjukkan tidak ada
overdispersi berat.
pearson_chi_square <- sum(residuals(poisson_best_fit, type="pearson")^2)
df_residual <- df.residual(poisson_best_fit)
dispersion_pearson <- pearson_chi_square/df_residual
overdispersion_p_value <- pchisq(pearson_chi_square, df=df_residual, lower.tail=FALSE)
poisson_overdispersion_table <- tibble(
uji="Pearson chi-square overdispersion",
pearson_chi_square=pearson_chi_square,
df_residual=df_residual,
dispersion_pearson=dispersion_pearson,
p_value=format_p(overdispersion_p_value),
keputusan=case_when(
dispersion_pearson < 1.5 ~ "Tidak ada indikasi overdispersion berat",
dispersion_pearson < 2.5 ~ "Ada indikasi overdispersion ringan-sedang",
TRUE ~ "Ada indikasi overdispersion kuat"
),
penjelasan_hasil=case_when(
dispersion_pearson < 1.5 ~ "Asumsi equidispersion relatif terpenuhi.",
dispersion_pearson < 2.5 ~ "Terdapat indikasi overdispersi ringan-sedang.",
TRUE ~ "Terdapat indikasi overdispersi kuat."
),
penanganan_jika_tidak_memenuhi=case_when(
dispersion_pearson < 1.5 ~ "Tidak diperlukan penanganan khusus.",
TRUE ~ "Gunakan quasi-Poisson, negative binomial, atau robust standard error."
)
)
safe_table(poisson_overdispersion_table, "Uji overdispersi Pearson model terbaik Poisson tahunan", digits=4, full_width=TRUE)6.6.2 Asumsi 2 — Proporsi Nilai Nol
Hipotesis/Kriteria:
\(H_0\): Proporsi nilai nol masih wajar
untuk model Poisson.
\(H_1\): Proporsi nilai nol terlalu
tinggi dan mengindikasikan potensi zero inflation.
Kriteria praktis: proporsi nol yang sangat besar perlu dibandingkan
dengan konteks data dan hasil prediksi model.
poisson_zero_table <- tibble(
ukuran=c("Jumlah observasi", "Jumlah case_count = 0", "Proporsi nilai nol"),
nilai=c(nrow(poisson_train), sum(poisson_train$case_count == 0, na.rm=TRUE), mean(poisson_train$case_count == 0, na.rm=TRUE)),
keterangan=c("Total observasi training", "Banyaknya observasi tanpa kejadian", "Proporsi nol yang tinggi dapat mengindikasikan zero inflation"),
penjelasan_hasil=c("Total observasi yang digunakan pada model terbaik.", "Banyaknya observasi tanpa kejadian pada data training.", "Jika proporsi nilai nol sangat tinggi, terdapat potensi zero inflation."),
penanganan_jika_tidak_memenuhi=c("Tidak diperlukan.", "Tidak diperlukan.", "Pertimbangkan zero-inflated Poisson, hurdle model, atau negative binomial jika juga terjadi overdispersi.")
)
safe_table(poisson_zero_table, "Pemeriksaan proporsi nilai nol pada data training Poisson", digits=4, full_width=TRUE)6.6.3 Asumsi 3 — Goodness-of-Fit (Deviance Test)
Hipotesis:
\(H_0\): Model Poisson memiliki
kecocokan yang memadai terhadap data.
\(H_1\): Model Poisson tidak memiliki
kecocokan yang memadai terhadap data.
Kriteria: p-value kecil menunjukkan indikasi lack of fit.
poisson_gof_table <- assumption_glm_table(poisson_best_fit, label="Poisson") |>
filter(uji == "Goodness-of-fit deviance") |>
mutate(p_value=p_value_tampil) |>
select(-p_value_tampil)
safe_table(poisson_gof_table, "Goodness-of-Fit Deviance Test model terbaik Poisson", digits=4, full_width=TRUE)6.6.4 Asumsi 4 — Multikolinearitas (VIF)
Hipotesis/Kriteria:
\(H_0\): Tidak terdapat
multikolinearitas berat antarprediktor.
\(H_1\): Terdapat multikolinearitas
berat antarprediktor.
Kriteria praktis: VIF < 5 aman; VIF 5–10 perlu perhatian; VIF > 10
menunjukkan multikolinearitas berat.
poisson_vif_table <- vif_table(poisson_train, poisson_selection$best$vars)
poisson_corr_best_table <- selected_correlation_check(poisson_train, poisson_selection$best$vars)
safe_table(poisson_vif_table, "Uji multikolinearitas model terbaik Poisson menggunakan VIF", digits=3, full_width=TRUE)safe_table(poisson_corr_best_table, "Pemeriksaan korelasi tinggi antarprediktor model terbaik Poisson", digits=3, full_width=TRUE)6.6.5 Asumsi 5 — Likelihood Ratio Test
Hipotesis:
\(H_0\): Semua koefisien prediktor sama
dengan nol; model dengan prediktor tidak lebih baik daripada model
kosong.
\(H_1\): Minimal ada satu koefisien
prediktor tidak sama dengan nol; model dengan prediktor lebih baik
daripada model kosong.
poisson_lrt_show <- poisson_best_g2 |> mutate(p_value=p_value_tampil) |> select(-p_value_tampil)
safe_table(poisson_lrt_show, "Likelihood Ratio Test atau uji G² model terbaik Poisson", digits=4, full_width=TRUE)6.6.6 Diagnostik Pengaruh Observasi
Diagnostik pengaruh observasi bukan uji hipotesis formal. Pemeriksaan ini digunakan untuk melihat apakah ada observasi yang terlalu memengaruhi estimasi model. Kriteria praktis: Cook’s distance > 1 perlu diperiksa, dan leverage tinggi dilihat dengan batas \(2p/n\).
poisson_influence_table <- influence_check_table(poisson_best_fit)
safe_table(poisson_influence_table, "Diagnostik pengaruh observasi model terbaik Poisson tahunan", digits=4, full_width=TRUE)6.7 Prediksi dengan model terbaik
poisson_pred_mu <- predict(poisson_best_fit, newdata=poisson_test, type="response")
poisson_pred_table <- tibble(actual=poisson_test$case_count, predicted_mu=poisson_pred_mu, residual=actual-predicted_mu, abs_error=abs(residual), squared_error=residual^2)
poisson_metrics <- poisson_pred_table |> summarise(MAE=mean(abs_error), RMSE=sqrt(mean(squared_error)), mean_actual=mean(actual), mean_predicted=mean(predicted_mu), correlation_actual_predicted=cor(actual,predicted_mu))
safe_table(poisson_metrics, "Metrik prediksi testing model terbaik Poisson tahunan", digits=4)poisson_pred_table |> ggplot(aes(x=predicted_mu, y=actual))+geom_point(alpha=.45,color="#172033")+geom_abline(slope=1,intercept=0,linetype="dashed",color="#e11d48")+labs(title="Aktual vs Prediksi Expected Count", subtitle="Titik yang dekat dengan garis diagonal menunjukkan prediksi yang lebih akurat.", x="Predicted expected count", y="Actual count")+theme_report()
Gambar 8. Scatterplot aktual versus prediksi expected count model terbaik Poisson.
poisson_main_x <- poisson_selected_vars[1]
if(is.na(poisson_main_x) || length(poisson_main_x)==0){poisson_main_x <- poisson_selection$best$vars[1]}
poisson_profile <- make_profile_grid(poisson_train, poisson_selection$best$vars, xvar=poisson_main_x, n=120)
poisson_profile_pred <- predict(poisson_best_fit, newdata=poisson_profile, type="link", se.fit=TRUE)
poisson_profile <- poisson_profile |>
mutate(
fit_link = poisson_profile_pred$fit,
se_link = poisson_profile_pred$se.fit,
predicted_mu = exp(fit_link),
ci_low = exp(fit_link - 1.96*se_link),
ci_high = exp(fit_link + 1.96*se_link)
)
poisson_profile |>
ggplot(aes(x=.data[[poisson_main_x]], y=predicted_mu))+
geom_ribbon(aes(ymin=ci_low, ymax=ci_high), fill="#6d28d9", alpha=.18)+
geom_line(linewidth=1.25, color="#6d28d9")+
labs(
title="Kurva Prediksi Expected Violence Case Count",
subtitle="Garis menunjukkan expected count; arsiran menunjukkan confidence interval 95%. Prediktor lain ditahan pada nilai median.",
x=poisson_main_x,
y="Predicted expected count"
)+
theme_report()
Gambar 9. Kurva prediksi expected count model terbaik Poisson dengan confidence interval 95 persen.
Prediksi Poisson menghasilkan expected count. MAE sebesar 0.372 berarti rata-rata kesalahan absolut prediksi sekitar nilai tersebut dalam satuan jumlah kejadian. Perbandingan mean aktual (0.42) dan mean prediksi (0.387) membantu membaca apakah model cenderung underpredict atau overpredict secara rata-rata.
6.8 Korelasi prediktor Poisson
poisson_corr_top <- correlation_pairs(poisson_train |> select(all_of(poisson_x)), top_n=10)
safe_table(poisson_corr_top, "Sepuluh pasangan prediktor Poisson dengan korelasi absolut terbesar", digits=3)7 Perbandingan Empat Model
comparison_table <- tibble(
model=c("Logistik biner","Logistik multinomial","Logistik ordinal","Poisson tahunan"),
respons=c("food_crisis","choice_code","gov_confidence_ord","case_count"),
sifat_respons=c("Biner","Nominal lebih dari dua kategori","Ordinal","Hitungan non-negatif"),
kriteria_model_terbaik=c("Jumlah variabel signifikan, lalu AIC","Jumlah variabel signifikan, lalu AIC","Jumlah variabel signifikan, lalu AIC","Jumlah variabel signifikan, lalu AIC"),
uji_model=c("G² likelihood ratio","G² likelihood ratio","G² likelihood ratio","G² likelihood ratio"),
uji_koefisien=c("Wald z-test","Wald z-test","Wald z-test","Wald z-test"),
ukuran_efek=c("Odds ratio","Relative risk ratio","Odds ratio","Incidence rate ratio"),
bentuk_prediksi=c("Probabilitas dan kelas 0/1","Probabilitas tiap kategori","Probabilitas tiap tingkat ordinal","Expected count")
)
safe_table(comparison_table, "Perbandingan model berdasarkan skala respons, pemilihan model, uji model, ukuran efek, dan bentuk prediksi", full_width=TRUE)Secara keseluruhan, model dipilih berdasarkan skala respons dan kemudian diseleksi berdasarkan jumlah variabel signifikan terlebih dahulu, baru AIC. Pendekatan ini menjaga agar model akhir tidak semata-mata mengejar AIC kecil, tetapi tetap mempertahankan variabel X yang secara statistik mendukung variabel Y. Setelah model terbaik diperoleh, interpretasi diarahkan pada variabel X yang signifikan, kemudian prediksi dilakukan memakai model terbaik pada data testing.