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)binary_effect_interpretation <- make_effect_interpretation(
binary_best_wald,
effect_col="odds_ratio",
effect_label="OR",
response_name="food_crisis"
)
safe_table(binary_effect_interpretation, "Interpretasi model terbaik logistik biner berdasarkan estimate dan odds ratio", digits=3, 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.6.7 Ringkasan Uji Asumsi
binary_assumption_summary <- tibble(
uji=c("VIF", "Likelihood Ratio Test", "Hosmer-Lemeshow", "Box-Tidwell", "Nagelkerke R-square", "Cook's distance"),
angka_hasil=c(
ifelse(all(is.na(binary_vif_table$VIF)), "VIF tidak dihitung", paste0("VIF maksimum = ", round(max(binary_vif_table$VIF, na.rm=TRUE), 3))),
paste0("G² = ", round(binary_best_g2$G2, 3), "; df = ", binary_best_g2$df, "; p-value = ", binary_best_g2$p_value_tampil),
paste0("HL = ", round(binary_hl$test$statistik, 3), "; df = ", binary_hl$test$df, "; p-value = ", binary_hl$test$p_value_tampil),
paste0("Jumlah indikasi non-linearitas = ", sum(binary_box_tidwell$keterangan == "Ada indikasi non-linearitas logit", na.rm=TRUE)),
paste0("Nagelkerke R² = ", round(binary_nagelkerke_table$nilai, 4)),
paste0("Maksimum Cook's distance = ", round(binary_influence_table$nilai[1], 4), "; jumlah Cook > 1 = ", binary_influence_table$nilai[2])
),
interpretasi_hasil=c(
ifelse(all(is.na(binary_vif_table$VIF)) || max(binary_vif_table$VIF, na.rm=TRUE) < 5, "Tidak ada indikasi multikolinearitas berat pada model terbaik.", "Terdapat indikasi multikolinearitas yang perlu diperhatikan."),
ifelse(binary_best_g2$p_value < .05, "Model terbaik signifikan secara simultan dan lebih baik daripada model kosong.", "Model terbaik belum terbukti lebih baik daripada model kosong secara simultan."),
ifelse(binary_hl$test$p_value < .05, "Terdapat indikasi lack of fit berdasarkan Hosmer-Lemeshow.", "Tidak ada bukti kuat lack of fit berdasarkan Hosmer-Lemeshow."),
ifelse(sum(binary_box_tidwell$keterangan == "Ada indikasi non-linearitas logit", na.rm=TRUE) > 0, "Ada prediktor yang berpotensi tidak linear pada skala logit.", "Tidak ada bukti kuat pelanggaran linearitas log-odds."),
"Nilai Nagelkerke R² menunjukkan proporsi peningkatan kecocokan model dibanding model kosong.",
ifelse(binary_influence_table$nilai[2] > 0, "Ada observasi yang berpengaruh kuat dan perlu dicek.", "Tidak ada observasi dengan Cook's distance > 1.")
),
kesimpulan=c(
ifelse(all(is.na(binary_vif_table$VIF)) || max(binary_vif_table$VIF, na.rm=TRUE) < 5, "Model layak digunakan karena tidak ada indikasi multikolinearitas berat.", "Model masih dapat digunakan, tetapi interpretasi koefisien perlu hati-hati karena ada indikasi multikolinearitas."),
ifelse(binary_best_g2$p_value < .05, "Model layak digunakan karena prediktor secara simultan mendukung food_crisis.", "Model dapat dilaporkan dengan catatan dukungan simultan belum kuat."),
ifelse(binary_hl$test$p_value < .05, "Model digunakan dengan catatan terdapat indikasi lack of fit.", "Model layak digunakan karena tidak ada bukti kuat lack of fit."),
ifelse(sum(binary_box_tidwell$keterangan == "Ada indikasi non-linearitas logit", na.rm=TRUE) > 0, "Model masih dapat digunakan, tetapi asumsi linearitas logit perlu diberi catatan.", "Model layak digunakan karena linearitas log-odds relatif terpenuhi."),
"Model memiliki kemampuan penjelasan yang dapat dibaca melalui Nagelkerke R².",
ifelse(binary_influence_table$nilai[2] > 0, "Model dapat digunakan dengan catatan ada observasi berpengaruh.", "Model layak digunakan karena tidak ada Cook's distance > 1.")
)
)
safe_table(binary_assumption_summary, "Ringkasan uji asumsi 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.
Interpretasi hasil prediksi: Pada data testing,
model terbaik menghasilkan accuracy sebesar 0.847, sensitivity sebesar
0.74, specificity sebesar 0.899, F1-score sebesar 0.761, dan AUC sebesar
0.921. Berdasarkan confusion matrix, model memprediksi benar
food_crisis sebanyak 97 observasi dan benar
non_food_crisis sebanyak 241 observasi. Kesalahan prediksi
terdiri dari false positive sebanyak 27 dan false negative sebanyak 34.
Dengan demikian, model dapat digunakan untuk prediksi karena performanya
tidak hanya dilihat dari accuracy, tetapi juga dari kemampuan mengenali
kelas kejadian melalui sensitivity dan kemampuan membedakan dua kelas
melalui AUC.
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)multinom_effect_interpretation <- make_effect_interpretation(
multinom_best_wald,
effect_col="relative_risk_ratio",
effect_label="RRR",
response_name="choice_code",
category_col="kategori_dibanding_referensi"
)
safe_table(multinom_effect_interpretation, "Interpretasi model terbaik multinomial berdasarkan estimate dan relative risk ratio", digits=3, 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.6.7 Ringkasan Uji Asumsi
multinom_assumption_summary <- tibble(
uji=c("Skala respons nominal", "Observasi independen", "VIF", "Perfect separation", "Box-Tidwell", "IIA"),
angka_hasil=c(
paste0("Jumlah kategori = ", n_distinct(multinom_train$choice_code), "; kelas minimum = ", min(table(multinom_train$choice_code))),
paste0("Jumlah observasi training = ", nrow(multinom_train)),
ifelse(all(is.na(multinom_vif_table$VIF)), "VIF tidak dihitung", paste0("VIF maksimum = ", round(max(multinom_vif_table$VIF, na.rm=TRUE), 3))),
paste0("Proporsi max probability > 0,99 = ", round(multinom_separation_table$nilai[1], 4), "; max probability = ", round(multinom_separation_table$nilai[2], 4)),
paste0("Jumlah indikasi non-linearitas = ", sum(multinom_bt_table$keterangan == "Ada indikasi non-linearitas logit", na.rm=TRUE)),
paste0("Kategori alternatif = ", paste(levels(multinom_train$choice_code), collapse=", "))
),
interpretasi_hasil=c(
"Respons memiliki kategori nominal yang saling eksklusif sehingga sesuai untuk multinomial.",
"Setiap baris diperlakukan sebagai satu unit observasi pada data model-ready.",
ifelse(all(is.na(multinom_vif_table$VIF)) || max(multinom_vif_table$VIF, na.rm=TRUE) < 5, "Tidak ada indikasi multikolinearitas berat.", "Terdapat indikasi multikolinearitas yang perlu diperhatikan."),
ifelse(multinom_separation_table$nilai[1] > .25, "Ada cukup banyak probabilitas prediksi ekstrem sehingga perlu waspada separasi.", "Tidak ada indikasi kuat perfect separation berdasarkan probabilitas ekstrem."),
ifelse(sum(multinom_bt_table$keterangan == "Ada indikasi non-linearitas logit", na.rm=TRUE) > 0, "Ada prediktor yang berpotensi tidak linear pada skala logit one-vs-rest.", "Tidak ada bukti kuat pelanggaran linearitas logit one-vs-rest."),
"IIA dievaluasi secara substantif karena uji formal memerlukan struktur/paket khusus."
),
kesimpulan=c(
"Model multinomial layak digunakan karena respons berbentuk nominal lebih dari dua kategori.",
"Model dapat digunakan dengan asumsi setiap observasi merupakan unit yang berbeda pada data model-ready.",
ifelse(all(is.na(multinom_vif_table$VIF)) || max(multinom_vif_table$VIF, na.rm=TRUE) < 5, "Model layak digunakan karena tidak ada indikasi multikolinearitas berat.", "Model masih dapat digunakan, tetapi interpretasi koefisien perlu hati-hati karena ada indikasi multikolinearitas."),
ifelse(multinom_separation_table$nilai[1] > .25, "Model dapat digunakan dengan catatan terdapat probabilitas prediksi ekstrem.", "Model layak digunakan karena tidak ada indikasi kuat perfect separation."),
ifelse(sum(multinom_bt_table$keterangan == "Ada indikasi non-linearitas logit", na.rm=TRUE) > 0, "Model masih dapat digunakan, tetapi asumsi linearitas logit perlu diberi catatan.", "Model layak digunakan karena linearitas logit one-vs-rest relatif terpenuhi."),
"Model dapat digunakan dengan catatan IIA dibaca sebagai asumsi substantif pada pilihan kategori."
)
)
safe_table(multinom_assumption_summary, "Ringkasan uji asumsi model terbaik multinomial", digits=4, 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)
multinom_class_summary <- multinom_metrics$class_metrics |> filter(!is.na(f1_score))
multinom_best_class <- if(nrow(multinom_class_summary)>0) multinom_class_summary$kelas[which.max(multinom_class_summary$f1_score)] else "-"
multinom_worst_class <- if(nrow(multinom_class_summary)>0) multinom_class_summary$kelas[which.min(multinom_class_summary$f1_score)] else "-"
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.
Interpretasi hasil prediksi: Prediksi multinomial dilakukan dengan memilih kategori dengan probabilitas terbesar. Pada data testing, accuracy model terbaik adalah 0.572 dan error rate sebesar 0.428. Berdasarkan metrik per kelas, kelas dengan F1-score terbaik adalah protect_environment, sedangkan kelas dengan F1-score terendah adalah protect_environment. Artinya, kemampuan prediksi model tidak merata pada semua kategori; kategori dengan F1-score lebih tinggi lebih mudah dikenali oleh model, sedangkan kategori dengan F1-score lebih rendah perlu dibaca hati-hati terutama bila jumlah observasinya lebih kecil atau polanya tumpang tindih dengan kategori lain.
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)ordinal_effect_interpretation <- make_effect_interpretation(
ordinal_best_wald,
effect_col="odds_ratio",
effect_label="OR",
response_name="gov_confidence_ord"
)
safe_table(ordinal_effect_interpretation, "Interpretasi model terbaik ordinal berdasarkan estimate dan odds ratio", digits=3, 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.6.3 Ringkasan Uji Asumsi
ordinal_assumption_summary <- tibble(
uji=c("Proportional odds", "VIF", "Korelasi tinggi antarprediktor"),
angka_hasil=c(
ifelse(nrow(ordinal_parallel$ringkas)==0, "Tidak dapat dihitung", paste0("Rentang koefisien maksimum = ", round(max(abs(ordinal_parallel$ringkas$rentang_estimate), na.rm=TRUE), 4))),
ifelse(all(is.na(ordinal_vif_table$VIF)), "VIF tidak dihitung", paste0("VIF maksimum = ", round(max(ordinal_vif_table$VIF, na.rm=TRUE), 3))),
paste0("Jumlah pasangan |r| >= 0,80 = ", sum(!is.na(ordinal_corr_best_table$abs_korelasi)))
),
interpretasi_hasil=c(
ifelse(nrow(ordinal_parallel$ringkas)==0, "Diagnostik proportional odds tidak dapat dihitung.", "Semakin kecil rentang koefisien, semakin stabil efek prediktor pada batas kumulatif kategori."),
ifelse(all(is.na(ordinal_vif_table$VIF)) || max(ordinal_vif_table$VIF, na.rm=TRUE) < 5, "Tidak ada indikasi multikolinearitas berat.", "Terdapat indikasi multikolinearitas yang perlu diperhatikan."),
ifelse(sum(!is.na(ordinal_corr_best_table$abs_korelasi)) == 0, "Tidak ada pasangan prediktor dengan korelasi tinggi.", "Ada pasangan prediktor dengan korelasi tinggi sehingga koefisien perlu ditafsirkan hati-hati.")
),
kesimpulan=c(
ifelse(nrow(ordinal_parallel$ringkas)==0, "Model ordinal digunakan dengan catatan diagnostik proportional odds tidak tersedia.", "Model ordinal dapat digunakan karena stabilitas koefisien kumulatif dapat dievaluasi melalui rentang koefisien."),
ifelse(all(is.na(ordinal_vif_table$VIF)) || max(ordinal_vif_table$VIF, na.rm=TRUE) < 5, "Model layak digunakan karena tidak ada indikasi multikolinearitas berat.", "Model masih dapat digunakan, tetapi interpretasi koefisien perlu hati-hati karena ada indikasi multikolinearitas."),
ifelse(sum(!is.na(ordinal_corr_best_table$abs_korelasi)) == 0, "Model layak digunakan karena tidak ada korelasi tinggi antarprediktor model terbaik.", "Model dapat digunakan dengan catatan ada korelasi tinggi antarprediktor.")
)
)
safe_table(ordinal_assumption_summary, "Ringkasan uji asumsi model terbaik ordinal", digits=4, 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)
ordinal_class_summary <- ordinal_metrics$class_metrics |> filter(!is.na(f1_score))
ordinal_best_class <- if(nrow(ordinal_class_summary)>0) ordinal_class_summary$kelas[which.max(ordinal_class_summary$f1_score)] else "-"
ordinal_worst_class <- if(nrow(ordinal_class_summary)>0) ordinal_class_summary$kelas[which.min(ordinal_class_summary$f1_score)] else "-"
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.
Interpretasi hasil prediksi: Accuracy testing model terbaik ordinal adalah 0.335 dengan error rate sebesar 0.665. Berdasarkan metrik per kelas, kategori dengan F1-score terbaik adalah quite_a_lot, sedangkan kategori dengan F1-score terendah adalah none_at_all. Karena respons bersifat ordinal, hasil prediksi tidak hanya dibaca dari benar-salah secara umum, tetapi juga dari pola kesalahannya: kesalahan ke kategori yang berdekatan lebih ringan secara substantif dibanding kesalahan ke kategori yang jauh.
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)poisson_effect_interpretation <- make_effect_interpretation(
poisson_best_wald,
effect_col="incidence_rate_ratio",
effect_label="IRR",
response_name="case_count"
)
safe_table(poisson_effect_interpretation, "Interpretasi model terbaik Poisson berdasarkan estimate dan incidence rate ratio", digits=3, 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."
),
kesimpulan=case_when(
dispersion_pearson < 1.5 ~ "Model Poisson layak digunakan karena tidak ada indikasi overdispersion berat.",
TRUE ~ "Model Poisson dapat digunakan dengan catatan terdapat indikasi overdispersion."
)
)
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."),
kesimpulan=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.6.7 Ringkasan Uji Asumsi
poisson_gof_p_numeric <- assumption_glm_table(poisson_best_fit, label="Poisson") |>
filter(uji == "Goodness-of-fit deviance") |>
pull(p_value)
poisson_assumption_summary <- tibble(
uji=c("Overdispersi Pearson", "Proporsi nilai nol", "Goodness-of-fit deviance", "VIF", "Likelihood Ratio Test", "Cook's distance"),
angka_hasil=c(
paste0("Dispersion Pearson = ", round(dispersion_pearson, 4), "; p-value = ", format_p(overdispersion_p_value)),
paste0("Proporsi nol = ", round(poisson_zero_table$nilai[3], 4), "; jumlah nol = ", poisson_zero_table$nilai[2]),
paste0("Deviance = ", round(poisson_gof_table$statistik, 3), "; df = ", poisson_gof_table$df, "; p-value = ", poisson_gof_table$p_value),
ifelse(all(is.na(poisson_vif_table$VIF)), "VIF tidak dihitung", paste0("VIF maksimum = ", round(max(poisson_vif_table$VIF, na.rm=TRUE), 3))),
paste0("G² = ", round(poisson_best_g2$G2, 3), "; df = ", poisson_best_g2$df, "; p-value = ", poisson_best_g2$p_value_tampil),
paste0("Maksimum Cook's distance = ", round(poisson_influence_table$nilai[1], 4), "; jumlah Cook > 1 = ", poisson_influence_table$nilai[2])
),
interpretasi_hasil=c(
ifelse(dispersion_pearson < 1.5, "Tidak ada indikasi overdispersi berat.", "Ada indikasi overdispersi sehingga standard error Poisson perlu dibaca hati-hati."),
ifelse(poisson_zero_table$nilai[3] > .50, "Proporsi nilai nol cukup tinggi dan dapat mengindikasikan zero inflation.", "Proporsi nilai nol tidak terlalu dominan berdasarkan batas praktis 50%."),
ifelse(poisson_gof_p_numeric < .05, "Ada indikasi lack of fit berdasarkan deviance.", "Tidak ada bukti kuat lack of fit berdasarkan deviance."),
ifelse(all(is.na(poisson_vif_table$VIF)) || max(poisson_vif_table$VIF, na.rm=TRUE) < 5, "Tidak ada indikasi multikolinearitas berat.", "Terdapat indikasi multikolinearitas yang perlu diperhatikan."),
ifelse(poisson_best_g2$p_value < .05, "Model terbaik signifikan secara simultan dibanding model kosong.", "Model terbaik belum signifikan secara simultan dibanding model kosong."),
ifelse(poisson_influence_table$nilai[2] > 0, "Ada observasi yang berpengaruh kuat dan perlu dicek.", "Tidak ada observasi dengan Cook's distance > 1.")
),
kesimpulan=c(
ifelse(dispersion_pearson < 1.5, "Model Poisson layak digunakan karena tidak ada indikasi overdispersi berat.", "Model Poisson dapat digunakan dengan catatan ada indikasi overdispersi."),
ifelse(poisson_zero_table$nilai[3] > .50, "Model dapat digunakan dengan catatan proporsi nol cukup tinggi.", "Model layak digunakan karena proporsi nol tidak terlalu dominan."),
ifelse(poisson_gof_p_numeric < .05, "Model digunakan dengan catatan terdapat indikasi lack of fit berdasarkan deviance.", "Model layak digunakan karena tidak ada bukti kuat lack of fit berdasarkan deviance."),
ifelse(all(is.na(poisson_vif_table$VIF)) || max(poisson_vif_table$VIF, na.rm=TRUE) < 5, "Model layak digunakan karena tidak ada indikasi multikolinearitas berat.", "Model masih dapat digunakan, tetapi interpretasi koefisien perlu hati-hati karena ada indikasi multikolinearitas."),
ifelse(poisson_best_g2$p_value < .05, "Model layak digunakan karena prediktor secara simultan mendukung expected count case_count.", "Model dapat dilaporkan dengan catatan dukungan simultan belum kuat."),
ifelse(poisson_influence_table$nilai[2] > 0, "Model dapat digunakan dengan catatan ada observasi berpengaruh.", "Model layak digunakan karena tidak ada Cook's distance > 1.")
)
)
safe_table(poisson_assumption_summary, "Ringkasan uji asumsi model terbaik Poisson", 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.
Interpretasi hasil prediksi: Prediksi Poisson menghasilkan expected count, bukan kategori. Pada data testing, MAE sebesar 0.372 menunjukkan bahwa rata-rata kesalahan absolut prediksi sekitar 0.372 kejadian. RMSE sebesar 0.784 menunjukkan besarnya kesalahan prediksi setelah memberi penalti lebih besar pada error yang ekstrem. Mean aktual adalah 0.42, sedangkan mean prediksi adalah 0.387. Korelasi aktual-prediksi sebesar 0.666 menunjukkan tingkat kesesuaian pola antara jumlah kejadian aktual dan expected count hasil model. Jika mean prediksi mendekati mean aktual dan korelasi bernilai positif, model dapat dipakai sebagai pendekatan prediksi rata-rata jumlah kejadian.
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.