Analisis Regresi Logistik
Regresi Logistik Biner, Multinomial, Ordinal, dan Poisson
David Christian
01 June 2026
1 Bagian I: Regresi Logistik Biner (Heart Disease)
required_packages <- c("dplyr", "ggplot2", "broom", "knitr", "scales", "tidyr")
missing_packages <- required_packages[
!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages) > 0) {
stop(
"Paket berikut belum tersedia: ",
paste(missing_packages, collapse = ", "),
". Silakan install terlebih dahulu."
)
}
invisible(lapply(required_packages, library, character.only = TRUE))2 Ringkasan Masalah
Penyakit jantung (heart disease) merupakan salah satu penyebab kematian tertinggi di dunia. Kemampuan untuk memprediksi apakah seseorang menderita penyakit jantung berdasarkan karakteristik klinis dan demografis memiliki nilai diagnostik yang tinggi.
Tujuan analisis ini adalah membangun model regresi logistik biner untuk memprediksi apakah seorang pasien terdiagnosis penyakit jantung, berdasarkan sejumlah prediktor klinis.
Kasus ini sesuai dianalisis dengan regresi logistik biner karena variabel respons bersifat biner:
- \(Y_i = 1\) : Pasien terdiagnosis penyakit jantung (heart disease)
- \(Y_i = 0\) : Pasien tidak terdiagnosis penyakit jantung
Model yang digunakan adalah:
\[ \log\left(\frac{p_i}{1 - p_i}\right) = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \cdots + \beta_k X_{ki} \]
di mana \(p_i = P(Y_i = 1 \mid X_i)\) adalah peluang pasien ke-\(i\) terdiagnosis penyakit jantung.
3 Sumber Data
Dataset yang digunakan adalah Heart Disease UCI yang tersedia di Kaggle:
Sumber: https://www.kaggle.com/datasets/redwankarimsony/heart-disease-data
Dataset ini merupakan gabungan dari empat repositori UCI: Cleveland, Hungary, Switzerland, dan VA Long Beach, dengan total 920 observasi dan 16 variabel.
# Sesuaikan path ke lokasi file CSV Anda
raw_heart <- read.csv(
file.path(dirname(knitr::current_input(dir = TRUE)), "heart_disease_uci.csv"),
sep = ";",
stringsAsFactors = FALSE
)
ringkasan_data <- data.frame(
Keterangan = c("Jumlah observasi", "Jumlah variabel"),
Nilai = c(nrow(raw_heart), ncol(raw_heart))
)
knitr::kable(ringkasan_data, caption = "Ukuran dataset Heart Disease UCI")| Keterangan | Nilai |
|---|---|
| Jumlah observasi | 920 |
| Jumlah variabel | 16 |
Interpretasi output: Dataset berisi 920 observasi
pasien dari empat pusat medis dengan 16 variabel. Variabel respons
adalah num, yaitu tingkat keparahan penyakit jantung (0–4),
yang akan dibinarisasi menjadi ada/tidaknya penyakit jantung.
4 Persiapan Data
4.1 Kamus Variabel
kamus <- data.frame(
`Nama variabel` = c("id", "age", "sex", "dataset", "cp", "trestbps",
"chol", "fbs", "restecg", "thalch", "exang",
"oldpeak", "slope", "ca", "thal", "num"),
`Keterangan` = c(
"Nomor identitas pasien",
"Usia pasien (tahun)",
"Jenis kelamin pasien",
"Asal pusat data (Cleveland, Hungary, Switzerland, VA Long Beach)",
"Tipe nyeri dada (chest pain type)",
"Tekanan darah istirahat saat masuk RS (mmHg)",
"Kolesterol serum (mg/dl)",
"Gula darah puasa > 120 mg/dl (TRUE/FALSE)",
"Hasil elektrokardiogram istirahat",
"Detak jantung maksimum yang tercapai",
"Angina akibat olahraga (TRUE/FALSE)",
"Depresi ST akibat olahraga relatif terhadap istirahat",
"Kemiringan segmen ST puncak olahraga",
"Jumlah pembuluh darah utama berwarna fluoroskopi (0–3)",
"Hasil thalassemia",
"Tingkat keparahan penyakit jantung (0 = tidak ada, 1–4 = ada)"
),
`Tipe` = c(
"ID", "Numerik", "Kategorik", "Kategorik", "Kategorik", "Numerik",
"Numerik", "Kategorik", "Kategorik", "Numerik", "Kategorik",
"Numerik", "Kategorik", "Numerik", "Kategorik", "Respons biner"
),
check.names = FALSE
)
knitr::kable(kamus, caption = "Kamus variabel dataset Heart Disease UCI")| Nama variabel | Keterangan | Tipe |
|---|---|---|
| id | Nomor identitas pasien | ID |
| age | Usia pasien (tahun) | Numerik |
| sex | Jenis kelamin pasien | Kategorik |
| dataset | Asal pusat data (Cleveland, Hungary, Switzerland, VA Long Beach) | Kategorik |
| cp | Tipe nyeri dada (chest pain type) | Kategorik |
| trestbps | Tekanan darah istirahat saat masuk RS (mmHg) | Numerik |
| chol | Kolesterol serum (mg/dl) | Numerik |
| fbs | Gula darah puasa > 120 mg/dl (TRUE/FALSE) | Kategorik |
| restecg | Hasil elektrokardiogram istirahat | Kategorik |
| thalch | Detak jantung maksimum yang tercapai | Numerik |
| exang | Angina akibat olahraga (TRUE/FALSE) | Kategorik |
| oldpeak | Depresi ST akibat olahraga relatif terhadap istirahat | Numerik |
| slope | Kemiringan segmen ST puncak olahraga | Kategorik |
| ca | Jumlah pembuluh darah utama berwarna fluoroskopi (0–3) | Numerik |
| thal | Hasil thalassemia | Kategorik |
| num | Tingkat keparahan penyakit jantung (0 = tidak ada, 1–4 = ada) | Respons biner |
4.2 Transformasi dan Pembersihan Data
Variabel respons num memiliki nilai 0–4. Nilai 0
menunjukkan tidak ada penyakit jantung, sedangkan nilai 1–4 menunjukkan
adanya penyakit jantung. Variabel ini dibinarisasi menjadi
heart_disease (0/1). Variabel id dan
dataset tidak digunakan sebagai prediktor model.
heart <- raw_heart %>%
transmute(
# Respons biner
heart_disease = as.integer(num > 0),
hd_label = factor(
ifelse(num > 0, "Penyakit Jantung", "Normal"),
levels = c("Normal", "Penyakit Jantung")
),
# Prediktor numerik
age = age,
trestbps = trestbps,
chol = chol,
thalch = thalch,
oldpeak = oldpeak,
ca = ca,
# Prediktor kategorik
sex = factor(sex, levels = c("Female", "Male")),
cp = factor(cp, levels = c("typical angina", "atypical angina",
"non-anginal", "asymptomatic")),
fbs = factor(ifelse(fbs == "TRUE" | fbs == TRUE, "Ya", "Tidak"),
levels = c("Tidak", "Ya")),
restecg = factor(restecg,
levels = c("normal", "st-t abnormality", "lv hypertrophy")),
exang = factor(ifelse(exang == "TRUE" | exang == TRUE, "Ya", "Tidak"),
levels = c("Tidak", "Ya")),
slope = factor(slope,
levels = c("upsloping", "flat", "downsloping")),
thal = factor(thal,
levels = c("normal", "fixed defect", "reversable defect"))
)
# Ringkasan missing values
missing_summary <- data.frame(
Variabel = names(heart),
`Jumlah NA` = sapply(heart, function(x) sum(is.na(x))),
`Persen NA` = sapply(heart, function(x) scales::percent(mean(is.na(x)), accuracy = 0.1)),
check.names = FALSE
) %>% filter(`Jumlah NA` > 0)
if (nrow(missing_summary) > 0) {
knitr::kable(missing_summary, caption = "Variabel dengan nilai hilang (missing values)")
} else {
cat("Tidak ada missing values pada variabel yang digunakan.\n")
}| Variabel | Jumlah NA | Persen NA | |
|---|---|---|---|
| trestbps | trestbps | 59 | 6.4% |
| chol | chol | 30 | 3.3% |
| thalch | thalch | 55 | 6.0% |
| oldpeak | oldpeak | 62 | 6.7% |
| ca | ca | 611 | 66.4% |
| fbs | fbs | 90 | 9.8% |
| restecg | restecg | 2 | 0.2% |
| exang | exang | 55 | 6.0% |
| slope | slope | 309 | 33.6% |
| thal | thal | 486 | 52.8% |
# Hapus baris dengan missing values pada variabel yang digunakan dalam model
heart_clean <- heart %>%
select(heart_disease, hd_label, age, trestbps, chol, thalch, oldpeak,
sex, cp, fbs, restecg, exang, slope, thal) %>%
na.omit()
cat(sprintf(
"Observasi sebelum na.omit: %d | Setelah na.omit: %d | Dihapus: %d\n",
nrow(heart), nrow(heart_clean), nrow(heart) - nrow(heart_clean)
))## Observasi sebelum na.omit: 920 | Setelah na.omit: 371 | Dihapus: 549Interpretasi output: Beberapa variabel memiliki
nilai hilang yang cukup banyak, terutama slope,
thal, dan ca. Setelah penghapusan baris dengan
nilai hilang, dataset final yang digunakan dalam pemodelan lebih bersih
dan siap dianalisis.
5 Eksplorasi Singkat
5.1 Distribusi Kelas Respons
class_summary <- heart_clean %>%
count(hd_label, name = "Jumlah") %>%
mutate(Proporsi = scales::percent(Jumlah / sum(Jumlah), accuracy = 0.1)) %>%
rename(`Status` = hd_label)
knitr::kable(class_summary, caption = "Distribusi kelas respons")| Status | Jumlah | Proporsi |
|---|---|---|
| Normal | 171 | 46.1% |
| Penyakit Jantung | 200 | 53.9% |
ggplot(heart_clean, aes(x = hd_label, fill = hd_label)) +
geom_bar(width = 0.6, color = "white", linewidth = 0.8) +
geom_text(
stat = "count",
aes(label = after_stat(count)),
vjust = -0.4, fontface = "bold"
) +
scale_fill_manual(values = c("Normal" = "#2a9d8f", "Penyakit Jantung" = "#e76f51")) +
labs(
title = "Distribusi Status Penyakit Jantung",
x = NULL,
y = "Jumlah pasien"
) +
theme_minimal(base_size = 12) +
theme(legend.position = "none")
Interpretasi output: Grafik batang memperlihatkan distribusi pasien berdasarkan status diagnosis. Proporsi masing-masing kelas menentukan apakah data cukup seimbang untuk pemodelan klasifikasi biner.
5.3 Scatter Plot: Usia vs Detak Jantung Maksimum
ggplot(heart_clean, aes(x = age, y = thalch, color = hd_label)) +
geom_point(alpha = 0.65, size = 2) +
scale_color_manual(values = c("Normal" = "#2a9d8f", "Penyakit Jantung" = "#e76f51")) +
labs(
title = "Usia vs Detak Jantung Maksimum",
subtitle = "Titik oranye menunjukkan pasien dengan penyakit jantung",
x = "Usia (tahun)",
y = "Detak jantung maksimum",
color = "Status"
) +
theme_minimal(base_size = 12)
Interpretasi output: Scatter plot memperlihatkan hubungan antara usia pasien dan detak jantung maksimum yang dicapai. Pasien dengan penyakit jantung cenderung memiliki detak jantung maksimum yang lebih rendah pada usia yang lebih tua, yang mengindikasikan bahwa kombinasi kedua variabel ini relevan sebagai prediktor.
6 Pembagian Data Training dan Testing
Data dibagi secara stratified dengan proporsi 80% training dan 20% testing agar komposisi kelas respons tetap seimbang di kedua subset.
stratified_split <- function(y, prop = 0.8) {
idx_by_class <- split(seq_along(y), y)
train_idx <- lapply(
idx_by_class,
function(idx) sample(idx, size = floor(length(idx) * prop))
)
unlist(train_idx, use.names = FALSE)
}
set.seed(42)
train_id <- stratified_split(heart_clean$heart_disease, prop = 0.8)
train_data <- heart_clean[train_id, ]
test_data <- heart_clean[-train_id, ]
split_summary <- bind_rows(
train_data %>% count(hd_label) %>% mutate(Data = "Training"),
test_data %>% count(hd_label) %>% mutate(Data = "Testing")
) %>%
group_by(Data) %>%
mutate(Proporsi = scales::percent(n / sum(n), accuracy = 0.1)) %>%
ungroup() %>%
rename(Status = hd_label, Jumlah = n) %>%
select(Data, Status, Jumlah, Proporsi)
knitr::kable(split_summary, caption = "Distribusi kelas pada training dan testing")| Data | Status | Jumlah | Proporsi |
|---|---|---|---|
| Training | Normal | 136 | 45.9% |
| Training | Penyakit Jantung | 160 | 54.1% |
| Testing | Normal | 35 | 46.7% |
| Testing | Penyakit Jantung | 40 | 53.3% |
Interpretasi output: Data training berisi 80% observasi dan data testing berisi 20% observasi. Karena split dilakukan secara stratified, proporsi kelas pada training dan testing mendekati proporsi data keseluruhan, sehingga evaluasi tidak bias akibat perbedaan komposisi kelas.
7 Model Regresi Logistik Biner
7.1 Rumus Model Logistik
Misalkan \(Y_i\) adalah status penyakit jantung pasien ke-\(i\), dengan:
\[ Y_i = \begin{cases} 1, & \text{jika terdiagnosis penyakit jantung},\\ 0, & \text{jika tidak terdiagnosis penyakit jantung}. \end{cases} \]
Peluang terdiagnosis penyakit jantung dinotasikan sebagai:
\[p_i = P(Y_i = 1 \mid X_i).\]
Pada regresi logistik biner, peluang dimodelkan melalui fungsi logit:
\[ \text{logit}(p_i) = \ln\left(\frac{p_i}{1-p_i}\right) = \eta_i = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \cdots + \beta_kX_{ki}. \]
Bentuk peluang prediksi diperoleh dengan transformasi balik:
\[ \hat{p}_i = \frac{\exp(\hat{\eta}_i)}{1 + \exp(\hat{\eta}_i)} = \frac{1}{1 + \exp(-\hat{\eta}_i)}. \]
Interpretasi odds ratio: Untuk prediktor numerik \(X_j\), odds ratio dihitung dengan:
\[OR_j = \exp(\beta_j).\]
- Jika \(OR_j > 1\) : kenaikan \(X_j\) meningkatkan odds penyakit jantung.
- Jika \(OR_j < 1\) : kenaikan \(X_j\) menurunkan odds penyakit jantung.
- Untuk prediktor kategorik, odds ratio dibandingkan terhadap kategori referensi.
7.2 Estimasi Model pada Data Training
heart_fit <- glm(
heart_disease ~ age + trestbps + chol + thalch + oldpeak +
sex + cp + fbs + restecg + exang + slope + thal,
data = train_data,
family = binomial(link = "logit")
)
ringkasan_model <- data.frame(
Keterangan = c(
"Jumlah observasi training",
"Null deviance",
"Residual deviance",
"Derajat bebas residual",
"AIC"
),
Nilai = c(
nobs(heart_fit),
round(heart_fit$null.deviance, 3),
round(heart_fit$deviance, 3),
heart_fit$df.residual,
round(AIC(heart_fit), 3)
)
)
knitr::kable(ringkasan_model, caption = "Ringkasan kecocokan model regresi logistik")| Keterangan | Nilai |
|---|---|
| Jumlah observasi training | 296.000 |
| Null deviance | 408.395 |
| Residual deviance | 210.101 |
| Derajat bebas residual | 278.000 |
| AIC | 246.101 |
Interpretasi output: Nilai residual deviance yang lebih kecil dari null deviance menunjukkan bahwa prediktor yang dimasukkan memberikan perbaikan nyata dibandingkan model kosong. Nilai AIC digunakan sebagai tolok ukur perbandingan antar model; semakin kecil AIC, semakin baik keseimbangan antara kecocokan model dan jumlah parameter.
7.3 Tabel Odds Ratio
coef_table <- broom::tidy(heart_fit) %>%
filter(term != "(Intercept)") %>%
mutate(
odds_ratio = exp(estimate),
ci_low = exp(estimate - 1.96 * std.error),
ci_high = exp(estimate + 1.96 * std.error)
) %>%
arrange(p.value) %>%
transmute(
`Variabel/level` = term,
`Odds ratio` = round(odds_ratio, 3),
`Interval kepercayaan 95 persen` = paste0(round(ci_low, 3), " – ", round(ci_high, 3)),
`p-value` = signif(p.value, 3)
)
knitr::kable(
head(coef_table, 20),
caption = "Dua puluh koefisien paling menonjol berdasarkan p-value"
)| Variabel/level | Odds ratio | Interval kepercayaan 95 persen | p-value |
|---|---|---|---|
| cpasymptomatic | 14.072 | 3.555 – 55.699 | 0.000165 |
| thalreversable defect | 4.261 | 1.898 – 9.57 | 0.000445 |
| sexMale | 2.755 | 1.119 – 6.781 | 0.027400 |
| oldpeak | 1.645 | 1.04 – 2.6 | 0.033200 |
| thalch | 0.981 | 0.963 – 1 | 0.047900 |
| trestbps | 1.020 | 0.999 – 1.042 | 0.058800 |
| chol | 0.995 | 0.99 – 1 | 0.059100 |
| cpatypical angina | 3.906 | 0.827 – 18.458 | 0.085500 |
| thalfixed defect | 3.512 | 0.775 – 15.927 | 0.103000 |
| slopeflat | 2.040 | 0.85 – 4.899 | 0.111000 |
| restecglv hypertrophy | 1.649 | 0.781 – 3.479 | 0.189000 |
| slopedownsloping | 0.383 | 0.062 – 2.349 | 0.300000 |
| age | 1.023 | 0.98 – 1.068 | 0.306000 |
| cpnon-anginal | 1.597 | 0.4 – 6.381 | 0.508000 |
| restecgst-t abnormality | 1.908 | 0.275 – 13.252 | 0.513000 |
| fbsYa | 0.756 | 0.259 – 2.209 | 0.609000 |
| exangYa | 0.946 | 0.413 – 2.167 | 0.895000 |
Interpretasi output: Tabel odds ratio diurutkan dari p-value terkecil. Prediktor dengan odds ratio jauh di atas 1 atau jauh di bawah 1 dan interval kepercayaan yang tidak melewati angka 1 menunjukkan pengaruh yang signifikan secara statistik terhadap probabilitas diagnosis penyakit jantung.
8 Prediksi dan Evaluasi Klasifikasi
8.1 Rumus Prediksi Kelas
Setelah model menghasilkan peluang prediksi \(\hat{p}_i\), status diagnosis ditentukan menggunakan threshold \(c\):
\[ \hat{Y}_i = \begin{cases} 1, & \text{jika } \hat{p}_i \ge c,\\ 0, & \text{jika } \hat{p}_i < c. \end{cases} \]
Jika \(c = 0{,}50\), pasien diklasifikasikan menderita penyakit jantung ketika peluang prediksinya minimal 50%.
p_train <- predict(heart_fit, newdata = train_data, type = "response")
p_test <- predict(heart_fit, newdata = test_data, type = "response")
prediction_preview <- head(
data.frame(
`Status aktual` = test_data$hd_label,
`Peluang prediksi (HD)` = round(p_test, 4),
check.names = FALSE
),
6
)
knitr::kable(prediction_preview, caption = "Contoh peluang prediksi penyakit jantung pada data testing")| Status aktual | Peluang prediksi (HD) | |
|---|---|---|
| 9 | Penyakit Jantung | 0.9597 |
| 12 | Normal | 0.2920 |
| 13 | Penyakit Jantung | 0.4985 |
| 14 | Normal | 0.2348 |
| 18 | Normal | 0.5576 |
| 23 | Penyakit Jantung | 0.4845 |
Interpretasi output: Semakin besar nilai peluang prediksi, semakin tinggi risiko pasien diklasifikasikan sebagai penderita penyakit jantung. Nilai ini belum menjadi kelas akhir hingga dibandingkan dengan threshold.
8.2 Rumus Confusion Matrix dan Metrik Evaluasi
Untuk klasifikasi biner penyakit jantung, notasi confusion matrix yang digunakan:
- TP (true positive): aktual penyakit jantung, diprediksi penyakit jantung.
- TN (true negative): aktual normal, diprediksi normal.
- FP (false positive): aktual normal, diprediksi penyakit jantung.
- FN (false negative): aktual penyakit jantung, diprediksi normal.
\[\text{Akurasi} = \frac{TP + TN}{TP + TN + FP + FN}\]
\[\text{Sensitivity} = \text{Recall} = \frac{TP}{TP + FN}\]
\[\text{Specificity} = \frac{TN}{TN + FP}\]
\[\text{Presisi} = \frac{TP}{TP + FP}\]
\[\text{F1-score} = 2 \times \frac{\text{Presisi} \times \text{Sensitivity}}{\text{Presisi} + \text{Sensitivity}}\]
\[\text{Balanced accuracy} = \frac{\text{Sensitivity} + \text{Specificity}}{2}\]
Dalam konteks penyakit jantung, sensitivity sangat penting karena menunjukkan kemampuan model mendeteksi pasien yang benar-benar menderita penyakit jantung sehingga tidak terlewat untuk mendapat penanganan.
safe_div <- function(num, den) {
ifelse(den == 0, NA_real_, num / den)
}
classification_metrics <- function(actual, prob, threshold = 0.5) {
pred <- as.integer(prob >= threshold)
tp <- sum(pred == 1 & actual == 1)
tn <- sum(pred == 0 & actual == 0)
fp <- sum(pred == 1 & actual == 0)
fn <- sum(pred == 0 & actual == 1)
sensitivity <- safe_div(tp, tp + fn)
specificity <- safe_div(tn, tn + fp)
precision <- safe_div(tp, tp + fp)
negative_predictive_value <- safe_div(tn, tn + fn)
accuracy <- safe_div(tp + tn, tp + tn + fp + fn)
data.frame(
threshold = threshold,
accuracy = accuracy,
error_rate = 1 - accuracy,
sensitivity = sensitivity,
specificity = specificity,
precision = precision,
negative_predictive_value = negative_predictive_value,
f1_score = safe_div(2 * precision * sensitivity,
precision + sensitivity),
balanced_accuracy = (sensitivity + specificity) / 2,
false_positive_rate = 1 - specificity,
false_negative_rate = 1 - sensitivity
)
}
format_metrics_indonesia <- function(x) {
x %>%
transmute(
Threshold = threshold,
Akurasi = accuracy,
`Error rate` = error_rate,
Sensitivity = sensitivity,
Specificity = specificity,
Presisi = precision,
NPV = negative_predictive_value,
`F1-score` = f1_score,
`Balanced accuracy` = balanced_accuracy,
FPR = false_positive_rate,
FNR = false_negative_rate
)
}
confusion_matrix_tbl <- function(actual, prob, threshold = 0.5) {
pred_label <- factor(
ifelse(prob >= threshold, "Prediksi HD", "Prediksi Normal"),
levels = c("Prediksi HD", "Prediksi Normal")
)
actual_label <- factor(
ifelse(actual == 1, "Aktual HD", "Aktual Normal"),
levels = c("Aktual HD", "Aktual Normal")
)
addmargins(table(actual_label, pred_label))
}8.3 Evaluasi pada Threshold 0,50
confusion_default <- confusion_matrix_tbl(test_data$heart_disease, p_test, threshold = 0.5)
metrics_default <- classification_metrics(test_data$heart_disease, p_test, threshold = 0.5) %>%
format_metrics_indonesia() %>%
mutate(across(where(is.numeric), round, 3))
knitr::kable(confusion_default, caption = "Confusion matrix testing pada threshold 0,50")| Prediksi HD | Prediksi Normal | Sum | |
|---|---|---|---|
| Aktual HD | 34 | 6 | 40 |
| Aktual Normal | 12 | 23 | 35 |
| Sum | 46 | 29 | 75 |
| Threshold | Akurasi | Error rate | Sensitivity | Specificity | Presisi | NPV | F1-score | Balanced accuracy | FPR | FNR |
|---|---|---|---|---|---|---|---|---|---|---|
| 0.5 | 0.76 | 0.24 | 0.85 | 0.657 | 0.739 | 0.793 | 0.791 | 0.754 | 0.343 | 0.15 |
Interpretasi output: Pada threshold 0,50, confusion matrix menunjukkan berapa pasien penyakit jantung yang berhasil dideteksi dan berapa pasien normal yang tetap dikenali sebagai normal. Sensitivity mengukur kemampuan model mendeteksi penyakit jantung; jika nilai ini rendah, banyak pasien yang benar-benar sakit akan terlewat oleh model.
9 Kurva ROC dan Threshold Optimal
Kurva ROC mengevaluasi trade-off antara:
- Sensitivity / True Positive Rate: proporsi pasien penyakit jantung yang berhasil dideteksi.
- False Positive Rate = \(1 - \text{Specificity}\): proporsi pasien normal yang keliru ditandai sakit.
Setiap titik pada kurva ROC berasal dari satu nilai threshold \(c\). Untuk setiap threshold dihitung:
\[TPR(c) = \frac{TP(c)}{TP(c) + FN(c)}, \quad FPR(c) = \frac{FP(c)}{FP(c) + TN(c)}.\]
Nilai AUC (area under the curve):
\[AUC = \int_0^1 TPR(FPR)\, d(FPR).\]
Semakin dekat AUC ke 1, semakin baik kemampuan model membedakan pasien sakit dan normal.
Threshold optimal dipilih menggunakan indeks Youden:
\[J = \text{Sensitivity} + \text{Specificity} - 1.\]
\[c_{\text{optimal}} = \arg\max_c \left\{ \text{Sensitivity}(c) + \text{Specificity}(c) - 1 \right\}.\]
roc_points <- function(actual, prob) {
thresholds <- c(Inf, sort(unique(prob), decreasing = TRUE), -Inf)
out <- lapply(thresholds, function(th) {
pred <- as.integer(prob >= th)
tp <- sum(pred == 1 & actual == 1)
tn <- sum(pred == 0 & actual == 0)
fp <- sum(pred == 1 & actual == 0)
fn <- sum(pred == 0 & actual == 1)
sensitivity <- safe_div(tp, tp + fn)
specificity <- safe_div(tn, tn + fp)
data.frame(
threshold = th,
sensitivity = sensitivity,
specificity = specificity,
fpr = 1 - specificity,
youden = sensitivity + specificity - 1
)
})
bind_rows(out)
}
auc_value <- function(roc_df) {
roc_sorted <- roc_df %>% arrange(fpr, sensitivity)
sum(
diff(roc_sorted$fpr) *
(head(roc_sorted$sensitivity, -1) + tail(roc_sorted$sensitivity, -1)) / 2
)
}
roc_train <- roc_points(train_data$heart_disease, p_train) %>% mutate(Data = "Training")
roc_test <- roc_points(test_data$heart_disease, p_test) %>% mutate(Data = "Testing")
auc_train <- auc_value(roc_train)
auc_test <- auc_value(roc_test)
optimal_train <- roc_train %>%
filter(is.finite(threshold)) %>%
arrange(desc(youden), desc(sensitivity)) %>%
slice(1)
threshold_opt <- optimal_train$threshold[1]
test_at_opt <- roc_points(test_data$heart_disease, p_test) %>%
filter(is.finite(threshold)) %>%
slice_min(abs(threshold - threshold_opt), n = 1, with_ties = FALSE) %>%
mutate(Data = "Testing — threshold optimal")
auc_table <- data.frame(
Data = c("Training", "Testing"),
AUC = round(c(auc_train, auc_test), 3)
)
threshold_table <- optimal_train %>%
transmute(
`Threshold optimal` = round(threshold, 3),
Sensitivity = round(sensitivity, 3),
Specificity = round(specificity, 3),
`Indeks Youden` = round(youden, 3)
)
knitr::kable(auc_table, caption = "Nilai AUC model regresi logistik")| Data | AUC |
|---|---|
| Training | 0.921 |
| Testing | 0.900 |
knitr::kable(threshold_table, caption = "Threshold optimal berdasarkan ROC training (Indeks Youden)")| Threshold optimal | Sensitivity | Specificity | Indeks Youden |
|---|---|---|---|
| 0.487 | 0.881 | 0.853 | 0.734 |
roc_plot <- bind_rows(roc_train, roc_test)
ggplot(roc_plot, aes(x = fpr, y = sensitivity, color = Data)) +
geom_path(linewidth = 1.1) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "#6c757d") +
geom_point(
data = optimal_train,
aes(x = fpr, y = sensitivity),
inherit.aes = FALSE,
color = "#ffb703", fill = "#fb8500",
shape = 21, size = 4, stroke = 1.2
) +
geom_point(
data = test_at_opt,
aes(x = fpr, y = sensitivity),
inherit.aes = FALSE,
color = "#8338ec", fill = "#3a86ff",
shape = 24, size = 4, stroke = 1.2
) +
coord_equal() +
scale_color_manual(values = c("Training" = "#0077b6", "Testing" = "#e76f51")) +
labs(
title = "Kurva ROC Model Regresi Logistik — Heart Disease",
subtitle = paste0(
"AUC training = ", round(auc_train, 3),
"; AUC testing = ", round(auc_test, 3),
"; threshold optimal = ", round(threshold_opt, 3)
),
x = "False positive rate (1 – Specificity)",
y = "Sensitivity / True positive rate",
color = "Data"
) +
theme_minimal(base_size = 12) +
theme(legend.position = "bottom")
Interpretasi output: AUC testing mengukur kemampuan model membedakan pasien penyakit jantung dan pasien normal pada data baru. Nilai AUC yang mendekati 1 menunjukkan performa yang baik. Jika kurva testing mendekati kurva training, model relatif stabil dan tidak mengalami overfitting yang parah. Titik yang ditandai pada kurva menunjukkan posisi threshold optimal berdasarkan indeks Youden.
10 Evaluasi dengan Threshold Optimal
Threshold optimal dari ROC training diterapkan pada data testing. Pendekatan ini lebih tepat karena threshold dipilih tanpa melihat data testing, sehingga evaluasi tetap bersifat out-of-sample.
metrics_compare <- bind_rows(
classification_metrics(test_data$heart_disease, p_test, threshold = 0.5) %>%
mutate(aturan = "Threshold 0,50"),
classification_metrics(test_data$heart_disease, p_test, threshold = threshold_opt) %>%
mutate(aturan = "Threshold optimal ROC")
) %>%
format_metrics_indonesia() %>%
bind_cols(`Aturan klasifikasi` = c("Threshold 0,50", "Threshold optimal ROC"), .) %>%
select(`Aturan klasifikasi`, everything()) %>%
mutate(across(where(is.numeric), round, 3))
confusion_opt <- confusion_matrix_tbl(test_data$heart_disease, p_test, threshold = threshold_opt)
knitr::kable(metrics_compare, caption = "Perbandingan metrik testing: threshold 0,50 vs threshold optimal")| Aturan klasifikasi | Threshold | Akurasi | Error rate | Sensitivity | Specificity | Presisi | NPV | F1-score | Balanced accuracy | FPR | FNR |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Threshold 0,50 | 0.500 | 0.760 | 0.240 | 0.850 | 0.657 | 0.739 | 0.793 | 0.791 | 0.754 | 0.343 | 0.150 |
| Threshold optimal ROC | 0.487 | 0.773 | 0.227 | 0.875 | 0.657 | 0.745 | 0.821 | 0.805 | 0.766 | 0.343 | 0.125 |
| Prediksi HD | Prediksi Normal | Sum | |
|---|---|---|---|
| Aktual HD | 35 | 5 | 40 |
| Aktual Normal | 12 | 23 | 35 |
| Sum | 47 | 28 | 75 |
Interpretasi output: Tabel perbandingan memperlihatkan dampak perubahan threshold terhadap sensitivity dan specificity. Pada threshold optimal, sensitivity meningkat sehingga model lebih agresif mendeteksi penyakit jantung; konsekuensinya specificity dapat turun karena sebagian pasien normal mungkin ikut diklasifikasikan sebagai sakit. Dalam konteks medis, pilihan threshold perlu mempertimbangkan biaya melewatkan pasien sakit (false negative) versus biaya salah mendiagnosis pasien sehat (false positive).
10.1 Distribusi Peluang Prediksi
test_prob_plot <- test_data %>% mutate(peluang_hd = p_test)
ggplot(test_prob_plot, aes(x = peluang_hd, fill = hd_label)) +
geom_density(alpha = 0.55, color = "white", linewidth = 0.7) +
geom_vline(
xintercept = threshold_opt,
color = "#fb8500", linewidth = 1.2, linetype = "dashed"
) +
annotate(
"label",
x = threshold_opt, y = Inf,
label = paste0("threshold = ", round(threshold_opt, 3)),
vjust = 1.4, fill = "#fff3b0", color = "#5f370e", label.size = 0
) +
scale_fill_manual(values = c("Normal" = "#2a9d8f", "Penyakit Jantung" = "#e76f51")) +
labs(
title = "Distribusi Peluang Prediksi Penyakit Jantung — Data Testing",
x = "Peluang prediksi penyakit jantung",
y = "Kepadatan",
fill = "Status aktual"
) +
theme_minimal(base_size = 12)
Interpretasi output: Grafik kepadatan memperlihatkan sebaran peluang prediksi untuk dua kelas. Semakin terpisah dua kurva, semakin baik model membedakan pasien sakit dan normal. Garis vertikal menunjukkan threshold optimal; observasi di sebelah kanannya diklasifikasikan sebagai penyakit jantung. Area tumpang tindih kedua kurva menunjukkan bagian data yang sulit dibedakan oleh model.
11 Interpretasi Hasil
Ringkasan hasil utama model regresi logistik biner untuk prediksi penyakit jantung:
ringkasan_akhir <- data.frame(
Metrik = c(
"AUC Testing",
"Threshold optimal (Youden)",
"Sensitivity pada threshold optimal",
"Specificity pada threshold optimal",
"Akurasi pada threshold optimal",
"Balanced accuracy pada threshold optimal"
),
Nilai = c(
round(auc_test, 3),
round(threshold_opt, 3),
round(classification_metrics(test_data$heart_disease, p_test, threshold_opt)$sensitivity, 3),
round(classification_metrics(test_data$heart_disease, p_test, threshold_opt)$specificity, 3),
round(classification_metrics(test_data$heart_disease, p_test, threshold_opt)$accuracy, 3),
round(classification_metrics(test_data$heart_disease, p_test, threshold_opt)$balanced_accuracy, 3)
)
)
knitr::kable(ringkasan_akhir, caption = "Ringkasan performa model pada data testing")| Metrik | Nilai |
|---|---|
| AUC Testing | 0.900 |
| Threshold optimal (Youden) | 0.487 |
| Sensitivity pada threshold optimal | 0.875 |
| Specificity pada threshold optimal | 0.657 |
| Akurasi pada threshold optimal | 0.773 |
| Balanced accuracy pada threshold optimal | 0.766 |
Secara substantif, model memprediksi peluang seorang pasien menderita penyakit jantung berdasarkan variabel klinis. Jika peluang prediksi berada di atas threshold optimal, pasien diklasifikasikan sebagai berisiko tinggi dan perlu penanganan lebih lanjut.
Dalam konteks klinis, threshold tidak harus selalu 0,50. Dokter atau sistem diagnostik yang lebih konservatif dapat memilih threshold lebih rendah agar lebih banyak kasus positif terdeteksi, meskipun konsekuensinya jumlah pasien sehat yang keliru diklasifikasikan sebagai sakit akan meningkat.
12 Kesimpulan
Regresi logistik biner dapat digunakan untuk membangun model prediksi penyakit jantung berbasis variabel klinis dan demografis. Evaluasi dengan pembagian training-testing memberikan gambaran kinerja model pada data baru yang belum pernah dilihat model sebelumnya, sedangkan kurva ROC dan indeks Youden membantu menentukan threshold klasifikasi yang lebih informatif dibandingkan menggunakan angka default 0,50 secara otomatis.
13 Bagian II: Regresi Logistik Multinomial (TravelMode)
# ── Paket yang digunakan ──────────────────────────────────────────────────────
library(AER) # data TravelMode
library(nnet) # multinom()
library(tidyverse) # manipulasi & visualisasi data
library(kableExtra) # tabel HTML
library(scales) # format label
library(GGally) # ggpairs
library(car) # vif
library(broom) # augment / tidy14 Pendahuluan
14.1 Latar Belakang
Pemilihan moda transportasi merupakan salah satu keputusan perilaku yang paling banyak diteliti dalam ekonomi transportasi. Individu memilih moda perjalanan berdasarkan pertimbangan biaya, waktu, kenyamanan, dan karakteristik rumah tangga. Karena pilihan moda terdiri atas lebih dari dua kategori dan tidak memiliki urutan alami, regresi logistik multinomial merupakan kerangka pemodelan yang paling sesuai.
14.2 Tujuan Analisis
Analisis ini bertujuan untuk:
- Mendeskripsikan distribusi pilihan moda transportasi beserta karakteristik prediktornya.
- Membangun model regresi logistik multinomial untuk menjelaskan faktor-faktor yang memengaruhi pilihan moda.
- Menginterpretasikan koefisien model dalam bentuk relative risk ratio (RRR).
- Mengevaluasi performa model dan menyajikan visualisasi probabilitas prediksi.
14.3 Data
Data yang digunakan adalah TravelMode dari paket
AER (Greene, 2003). Data ini memuat observasi pilihan
moda perjalanan antara Sydney dan Melbourne, Australia, untuk 210
individu terhadap empat moda: mobil (car),
pesawat (air), kereta
(train), dan bus (bus).
Format asli data berbentuk long (840 baris = 210 individu × 4 moda). Untuk keperluan regresi logistik multinomial, data akan diubah ke format wide sehingga setiap baris mewakili satu individu dengan pilihan yang dilakukan.
15 Struktur Data
## Rows: 840
## Columns: 9
## $ individual <fct> 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5,…
## $ mode <fct> air, train, bus, car, air, train, bus, car, air, train, bus…
## $ choice <fct> no, no, no, yes, no, no, no, yes, no, no, no, yes, no, no, …
## $ wait <int> 69, 34, 35, 0, 64, 44, 53, 0, 69, 34, 35, 0, 64, 44, 53, 0,…
## $ vcost <int> 59, 31, 25, 10, 58, 31, 25, 11, 115, 98, 53, 23, 49, 26, 21…
## $ travel <int> 100, 372, 417, 180, 68, 354, 399, 255, 125, 892, 882, 720, …
## $ gcost <int> 70, 71, 70, 30, 68, 84, 85, 50, 129, 195, 149, 101, 59, 79,…
## $ income <int> 35, 35, 35, 35, 30, 30, 30, 30, 40, 40, 40, 40, 70, 70, 70,…
## $ size <int> 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 3, 3, 3, 3, 2, 2, 2, 2,…# Contoh baris pertama
head(TravelMode, 12) %>%
kable(caption = "Contoh 12 baris pertama data TravelMode (format long)") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE)| individual | mode | choice | wait | vcost | travel | gcost | income | size |
|---|---|---|---|---|---|---|---|---|
| 1 | air | no | 69 | 59 | 100 | 70 | 35 | 1 |
| 1 | train | no | 34 | 31 | 372 | 71 | 35 | 1 |
| 1 | bus | no | 35 | 25 | 417 | 70 | 35 | 1 |
| 1 | car | yes | 0 | 10 | 180 | 30 | 35 | 1 |
| 2 | air | no | 64 | 58 | 68 | 68 | 30 | 2 |
| 2 | train | no | 44 | 31 | 354 | 84 | 30 | 2 |
| 2 | bus | no | 53 | 25 | 399 | 85 | 30 | 2 |
| 2 | car | yes | 0 | 11 | 255 | 50 | 30 | 2 |
| 3 | air | no | 69 | 115 | 125 | 129 | 40 | 1 |
| 3 | train | no | 34 | 98 | 892 | 195 | 40 | 1 |
| 3 | bus | no | 35 | 53 | 882 | 149 | 40 | 1 |
| 3 | car | yes | 0 | 23 | 720 | 101 | 40 | 1 |
Keterangan variabel:
| Variabel | Keterangan |
|---|---|
individual |
Identitas individu (1–210) |
mode |
Moda transportasi: air, bus,
car, train |
choice |
Apakah moda ini yang dipilih? (yes /
no) |
wait |
Waktu tunggu di terminal (menit); 0 untuk mobil |
vcost |
Komponen biaya kendaraan |
travel |
Waktu perjalanan dalam kendaraan (menit) |
gcost |
Ukuran biaya umum (generalized cost) |
income |
Pendapatan rumah tangga |
size |
Ukuran rombongan perjalanan (party size) |
16 Persiapan Data
Karena regresi logistik multinomial membutuhkan satu baris per
individu, data diubah dari format long ke format wide.
Variabel yang bersifat individual-specific (sama untuk
semua moda per individu) — yaitu income dan
size — diambil langsung.
Untuk variabel alternative-specific (berbeda per
moda), diambil nilai dari moda yang benar-benar dipilih
oleh masing-masing individu, karena dalam regresi logistik multinomial
standar (nnet::multinom), prediktor yang dimasukkan
bersifat individual-level.
Catatan: Pendekatan ini menggunakan atribut dari moda yang dipilih sebagai prediktor individual. Untuk model conditional logit (yang mempertimbangkan atribut semua alternatif), diperlukan pendekatan berbeda seperti
mlogit. Analisis ini berfokus pada regresi logistik multinomial standar sesuai landasan materi.
# ── Pilihan aktual tiap individu ─────────────────────────────────────────────
chosen <- TravelMode %>%
filter(choice == "yes") %>%
select(individual, mode, wait, vcost, travel, gcost)
# ── Atribut individual-specific ──────────────────────────────────────────────
indiv_attr <- TravelMode %>%
distinct(individual, income, size)
# ── Gabungkan ─────────────────────────────────────────────────────────────────
data_wide <- chosen %>%
left_join(indiv_attr, by = "individual") %>%
mutate(
mode = factor(mode, levels = c("car", "air", "train", "bus"))
)
glimpse(data_wide)## Rows: 210
## Columns: 8
## $ individual <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
## $ mode <fct> car, car, car, car, car, train, air, car, car, car, car, ca…
## $ wait <int> 0, 0, 0, 0, 0, 40, 45, 0, 0, 0, 0, 0, 0, 0, 0, 20, 15, 20, …
## $ vcost <int> 10, 11, 23, 5, 8, 20, 148, 50, 17, 7, 4, 6, 5, 15, 17, 19, …
## $ travel <int> 180, 255, 720, 180, 600, 345, 115, 780, 210, 210, 210, 250,…
## $ gcost <int> 30, 50, 101, 32, 99, 57, 160, 135, 40, 30, 36, 44, 41, 58, …
## $ income <int> 35, 30, 40, 70, 45, 20, 45, 12, 40, 70, 15, 35, 50, 40, 26,…
## $ size <int> 1, 2, 1, 3, 2, 1, 1, 1, 1, 2, 2, 2, 4, 1, 4, 1, 1, 1, 1, 2,…# Distribusi moda yang dipilih
data_wide %>%
count(mode) %>%
mutate(proporsi = n / sum(n)) %>%
kable(
digits = 3,
caption = "Distribusi moda transportasi yang dipilih"
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE)| mode | n | proporsi |
|---|---|---|
| car | 59 | 0.281 |
| air | 58 | 0.276 |
| train | 63 | 0.300 |
| bus | 30 | 0.143 |
data_wide %>%
count(mode) %>%
mutate(proporsi = n / sum(n)) %>%
ggplot(aes(x = mode, y = proporsi, fill = mode)) +
geom_col(width = 0.65, alpha = 0.92, color = "white") +
geom_text(aes(label = percent(proporsi, accuracy = 0.1)),
vjust = -0.45, fontface = "bold", color = "#243142", size = 4.2) +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.55)) +
scale_fill_manual(values = c(
"car" = "#2f7f73",
"air" = "#5568b8",
"train" = "#d18b2f",
"bus" = "#b84f5a"
)) +
labs(
title = "Distribusi Pilihan Moda Transportasi",
subtitle = "Respons multinomial: empat kategori tanpa urutan alami.",
x = NULL, y = "Proporsi"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"),
panel.grid.major.x = element_blank())
Distribusi pilihan moda transportasi
17 Eksplorasi Data
17.1 Statistik Deskriptif
data_wide %>%
select(wait, vcost, travel, gcost, income, size) %>%
pivot_longer(everything(), names_to = "variabel", values_to = "nilai") %>%
group_by(variabel) %>%
summarise(
N = n(),
Mean = mean(nilai, na.rm = TRUE),
SD = sd(nilai, na.rm = TRUE),
Min = min(nilai, na.rm = TRUE),
Q1 = quantile(nilai, 0.25, na.rm = TRUE),
Median = median(nilai, na.rm = TRUE),
Q3 = quantile(nilai, 0.75, na.rm = TRUE),
Max = max(nilai, na.rm = TRUE)
) %>%
mutate(across(where(is.numeric), ~round(.x, 2))) %>%
kable(caption = "Statistik deskriptif variabel prediktor") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE)| variabel | N | Mean | SD | Min | Q1 | Median | Q3 | Max |
|---|---|---|---|---|---|---|---|---|
| gcost | 210 | 103.82 | 45.52 | 30 | 67 | 102.5 | 132.75 | 238 |
| income | 210 | 34.55 | 19.71 | 2 | 20 | 34.5 | 50.00 | 72 |
| size | 210 | 1.74 | 1.01 | 1 | 1 | 1.0 | 2.00 | 6 |
| travel | 210 | 430.85 | 302.50 | 65 | 180 | 305.0 | 720.00 | 1440 |
| vcost | 210 | 47.40 | 38.52 | 2 | 19 | 33.0 | 70.00 | 180 |
| wait | 210 | 25.01 | 24.76 | 0 | 0 | 20.0 | 40.00 | 99 |
17.2 Distribusi Prediktor per Moda
data_wide %>%
select(mode, vcost, travel, gcost, income, size, wait) %>%
pivot_longer(-mode, names_to = "variabel", values_to = "nilai") %>%
mutate(variabel = dplyr::recode(variabel,
vcost = "Biaya Kendaraan\n(vcost)",
travel = "Waktu Perjalanan\n(travel)",
gcost = "Biaya Umum\n(gcost)",
income = "Pendapatan\n(income)",
size = "Ukuran Rombongan\n(size)",
wait = "Waktu Tunggu\n(wait)"
)) %>%
ggplot(aes(x = mode, y = nilai, fill = mode)) +
geom_boxplot(alpha = 0.80, outlier.size = 1.2, color = "#243142") +
facet_wrap(~variabel, scales = "free_y", ncol = 3) +
scale_fill_manual(values = c(
"car" = "#2f7f73",
"air" = "#5568b8",
"train" = "#d18b2f",
"bus" = "#b84f5a"
)) +
labs(
title = "Distribusi Prediktor per Moda Transportasi",
subtitle = "Perbedaan karakteristik antar moda dapat diamati dari distribusi masing-masing variabel.",
x = NULL, y = "Nilai"
) +
theme_minimal(base_size = 11) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"),
strip.text = element_text(face = "bold"),
panel.grid.major.x = element_blank())
Distribusi prediktor berdasarkan moda pilihan
## character(0)17.3 Matriks Korelasi Prediktor
data_wide %>%
select(wait, vcost, travel, gcost, income, size) %>%
cor(use = "complete.obs") %>%
round(3) %>%
as.data.frame() %>%
rownames_to_column("Variabel") %>%
kable(caption = "Matriks korelasi prediktor kontinu") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE)| Variabel | wait | vcost | travel | gcost | income | size |
|---|---|---|---|---|---|---|
| wait | 1.000 | 0.631 | -0.224 | 0.318 | -0.061 | -0.143 |
| vcost | 0.631 | 1.000 | -0.378 | 0.503 | 0.115 | -0.211 |
| travel | -0.224 | -0.378 | 1.000 | 0.576 | -0.125 | 0.147 |
| gcost | 0.318 | 0.503 | 0.576 | 1.000 | 0.034 | 0.115 |
| income | -0.061 | 0.115 | -0.125 | 0.034 | 1.000 | 0.181 |
| size | -0.143 | -0.211 | 0.147 | 0.115 | 0.181 | 1.000 |
18 Pemilihan Model
18.1 Landasan Pemilihan
Variabel respons mode memiliki empat kategori:
car, air, train, dan
bus. Karena keempat kategori ini tidak memiliki
urutan alami (seseorang yang memilih kereta tidak lebih
“tinggi” atau “rendah” dari yang memilih bus), model yang paling tepat
adalah regresi logistik multinomial dengan pendekatan
baseline-category logit.
Dalam model ini, satu kategori dijadikan referensi
(baseline), dan setiap kategori lain dibandingkan terhadap referensi
tersebut. Kita gunakan "car" sebagai kategori referensi
karena moda mobil merupakan pilihan paling umum dan secara substantif
menjadi pembanding yang bermakna.
18.2 Variabel Prediktor
Berdasarkan eksplorasi data dan relevansi substantif, prediktor yang dimasukkan ke dalam model adalah:
| Prediktor | Deskripsi | Tipe |
|---|---|---|
vcost |
Biaya kendaraan untuk moda yang dipilih | Kontinu |
travel |
Waktu perjalanan (menit) untuk moda yang dipilih | Kontinu |
wait |
Waktu tunggu di terminal (menit) | Kontinu |
income |
Pendapatan rumah tangga | Kontinu |
size |
Ukuran rombongan perjalanan | Kontinu |
19 Estimasi Model Regresi Logistik Multinomial
# Kategori referensi: car
data_wide$mode <- relevel(data_wide$mode, ref = "car")
fit_multi <- nnet::multinom(
mode ~ vcost + travel + wait + income + size,
data = data_wide,
trace = FALSE
)
summary(fit_multi)## Call:
## nnet::multinom(formula = mode ~ vcost + travel + wait + income +
## size, data = data_wide, trace = FALSE)
##
## Coefficients:
## (Intercept) vcost travel wait income size
## air -32.52272 3.2129081 -0.695752675 15.28274 -0.5120975 -10.61325
## train 18.91076 -0.8280502 0.003910910 14.74282 -0.3796165 -12.89153
## bus 18.11377 -0.8556775 0.006261645 14.73063 -0.3496898 -13.42680
##
## Std. Errors:
## (Intercept) vcost travel
## air 0.0000000000001691329 0.0000000000095987 0.000000000008876543
## train 0.4096255246816904450 0.0156019976188049 0.204650950999909487
## bus 0.4096243955256760327 0.0155999007821853 0.204650632848070030
## wait income size
## air 0.0000000000006035074 0.0000000000006619418 0.00000000000003205577
## train 0.0079585097907365881 0.0124738198484016273 0.18380082933668029366
## bus 0.0079573983270127893 0.0124629082048986952 0.18380172849967071902
##
## Residual Deviance: 103.733
## AIC: 139.73319.1 Ringkasan Koefisien, SE, z-value, dan p-value
Fungsi multinom() tidak langsung mengeluarkan p-value.
Nilai tersebut dihitung dari:
\[z = \frac{\hat{\beta}}{SE(\hat{\beta})}, \qquad p\text{-value} = 2\{1 - \Phi(|z|)\}.\]
multi_sum <- summary(fit_multi)
coef_long <- as.data.frame(multi_sum$coefficients) %>%
rownames_to_column("kategori") %>%
pivot_longer(-kategori, names_to = "variabel", values_to = "estimate")
se_long <- as.data.frame(multi_sum$standard.errors) %>%
rownames_to_column("kategori") %>%
pivot_longer(-kategori, names_to = "variabel", values_to = "std_error")
result_multi <- coef_long %>%
left_join(se_long, by = c("kategori", "variabel")) %>%
mutate(
z_value = estimate / std_error,
p_value = 2 * (1 - pnorm(abs(z_value))),
RRR = exp(estimate),
CI_low = exp(estimate - 1.96 * std_error),
CI_high = exp(estimate + 1.96 * std_error),
sig = case_when(
p_value < 0.001 ~ "***",
p_value < 0.01 ~ "**",
p_value < 0.05 ~ "*",
p_value < 0.10 ~ ".",
TRUE ~ ""
)
)
result_multi %>%
mutate(across(c(estimate, std_error, z_value, p_value, RRR, CI_low, CI_high),
~round(.x, 4))) %>%
kable(
caption = "Ringkasan hasil regresi logistik multinomial (kategori referensi: car)",
col.names = c("Kategori", "Variabel", "Estimate", "SE",
"z-value", "p-value", "RRR", "CI 2.5%", "CI 97.5%", "Sig.")
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE) %>%
column_spec(7, bold = TRUE) %>%
footnote(general = "Sig.: *** p<0.001 ** p<0.01 * p<0.05 . p<0.10",
general_title = "Keterangan:")| Kategori | Variabel | Estimate | SE | z-value | p-value | RRR | CI 2.5% | CI 97.5% | Sig. |
|---|---|---|---|---|---|---|---|---|---|
| air | (Intercept) | -32.5227 | 0.0000 | -192290925733462.5938 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | *** |
| air | vcost | 3.2129 | 0.0000 | 334723241166.2457 | 0.0000 | 24.8513 | 24.8513 | 24.8513 | *** |
| air | travel | -0.6958 | 0.0000 | -78381036482.7353 | 0.0000 | 0.4987 | 0.4987 | 0.4987 | *** |
| air | wait | 15.2827 | 0.0000 | 25323201428283.4375 | 0.0000 | 4337194.4149 | 4337194.4149 | 4337194.4149 | *** |
| air | income | -0.5121 | 0.0000 | -773629109039.7035 | 0.0000 | 0.5992 | 0.5992 | 0.5992 | *** |
| air | size | -10.6133 | 0.0000 | -331087129829500.8750 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | *** |
| train | (Intercept) | 18.9108 | 0.4096 | 46.1660 | 0.0000 | 163244337.4447 | 73140485.6552 | 364349695.9180 | *** |
| train | vcost | -0.8281 | 0.0156 | -53.0733 | 0.0000 | 0.4369 | 0.4237 | 0.4505 | *** |
| train | travel | 0.0039 | 0.2047 | 0.0191 | 0.9848 | 1.0039 | 0.6722 | 1.4993 | |
| train | wait | 14.7428 | 0.0080 | 1852.4597 | 0.0000 | 2527696.0927 | 2488573.2975 | 2567433.9363 | *** |
| train | income | -0.3796 | 0.0125 | -30.4331 | 0.0000 | 0.6841 | 0.6676 | 0.7011 | *** |
| train | size | -12.8915 | 0.1838 | -70.1386 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | *** |
| bus | (Intercept) | 18.1138 | 0.4096 | 44.2204 | 0.0000 | 73571541.8546 | 32963288.6197 | 164206060.6666 | *** |
| bus | vcost | -0.8557 | 0.0156 | -54.8515 | 0.0000 | 0.4250 | 0.4122 | 0.4382 | *** |
| bus | travel | 0.0063 | 0.2047 | 0.0306 | 0.9756 | 1.0063 | 0.6738 | 1.5029 | |
| bus | wait | 14.7306 | 0.0080 | 1851.1870 | 0.0000 | 2497079.7822 | 2458436.2111 | 2536330.7822 | *** |
| bus | income | -0.3497 | 0.0125 | -28.0584 | 0.0000 | 0.7049 | 0.6879 | 0.7223 | *** |
| bus | size | -13.4268 | 0.1838 | -73.0505 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | *** |
| Keterangan: | |||||||||
| Sig.: *** p<0.001 ** p<0.01 * p<0.05 . p<0.10 |
20 Persamaan Model
Dengan kategori referensi car, model regresi logistik multinomial menghasilkan tiga persamaan logit:
\[\log\left(\frac{\hat{\pi}_{\text{air}}}{\hat{\pi}_{\text{car}}}\right) = \hat{\alpha}_{\text{air}} + \hat{\beta}_{\text{air},1}\,\textit{vcost} + \hat{\beta}_{\text{air},2}\,\textit{travel} + \hat{\beta}_{\text{air},3}\,\textit{wait} + \hat{\beta}_{\text{air},4}\,\textit{income} + \hat{\beta}_{\text{air},5}\,\textit{size}\]
\[\log\left(\frac{\hat{\pi}_{\text{train}}}{\hat{\pi}_{\text{car}}}\right) = \hat{\alpha}_{\text{train}} + \hat{\beta}_{\text{train},1}\,\textit{vcost} + \hat{\beta}_{\text{train},2}\,\textit{travel} + \hat{\beta}_{\text{train},3}\,\textit{wait} + \hat{\beta}_{\text{train},4}\,\textit{income} + \hat{\beta}_{\text{train},5}\,\textit{size}\]
\[\log\left(\frac{\hat{\pi}_{\text{bus}}}{\hat{\pi}_{\text{car}}}\right) = \hat{\alpha}_{\text{bus}} + \hat{\beta}_{\text{bus},1}\,\textit{vcost} + \hat{\beta}_{\text{bus},2}\,\textit{travel} + \hat{\beta}_{\text{bus},3}\,\textit{wait} + \hat{\beta}_{\text{bus},4}\,\textit{income} + \hat{\beta}_{\text{bus},5}\,\textit{size}\]
# Cetak nilai estimasi per kategori
coef_matrix <- round(coef(fit_multi), 4)
coef_matrix %>%
as.data.frame() %>%
rownames_to_column("Kategori vs. Car") %>%
kable(caption = "Matriks koefisien model (log-odds terhadap car)") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE)| Kategori vs. Car | (Intercept) | vcost | travel | wait | income | size |
|---|---|---|---|---|---|---|
| air | -32.5227 | 3.2129 | -0.6958 | 15.2827 | -0.5121 | -10.6133 |
| train | 18.9108 | -0.8281 | 0.0039 | 14.7428 | -0.3796 | -12.8915 |
| bus | 18.1138 | -0.8557 | 0.0063 | 14.7306 | -0.3497 | -13.4268 |
21 Interpretasi Koefisien
21.1 Relative Risk Ratio (RRR)
Koefisien pada model logistik multinomial diinterpretasikan melalui Relative Risk Ratio \(RRR = \exp(\hat{\beta})\):
- Jika \(RRR > 1\): kenaikan prediktor meningkatkan odds relatif memilih kategori tersebut dibanding kategori referensi (car).
- Jika \(RRR < 1\): kenaikan prediktor menurunkan odds relatif memilih kategori tersebut dibanding referensi.
result_multi %>%
filter(variabel != "(Intercept)") %>%
mutate(
sig_label = ifelse(p_value < 0.05, "Signifikan (p < 0.05)", "Tidak Signifikan"),
variabel = factor(variabel, levels = c("vcost","travel","wait","income","size"))
) %>%
ggplot(aes(x = variabel, y = RRR, color = sig_label)) +
geom_hline(yintercept = 1, linetype = "dashed", color = "#64748b", linewidth = 0.8) +
geom_pointrange(aes(ymin = CI_low, ymax = CI_high),
size = 0.7, linewidth = 1.1,
position = position_dodge(width = 0.3)) +
facet_wrap(~kategori, ncol = 3) +
scale_color_manual(values = c("Signifikan (p < 0.05)" = "#2f7f73",
"Tidak Signifikan" = "#b0b8c1")) +
scale_y_log10() +
labs(
title = "Relative Risk Ratio (RRR) — Regresi Logistik Multinomial",
subtitle = "Garis putus-putus pada RRR = 1 (tidak ada efek). Skala log pada sumbu y.",
x = "Prediktor",
y = "RRR (skala log)",
color = NULL
) +
theme_minimal(base_size = 12) +
theme(legend.position = "bottom",
plot.title = element_text(face = "bold"),
strip.text = element_text(face = "bold", color = "#243142"),
panel.grid.major.x = element_blank())
Relative Risk Ratio (RRR) dengan interval kepercayaan 95%
21.2 Interpretasi per Kategori
21.2.1 Air vs. Car
result_multi %>%
filter(kategori == "air", variabel != "(Intercept)") %>%
select(variabel, estimate, RRR, CI_low, CI_high, p_value, sig) %>%
mutate(across(c(estimate, RRR, CI_low, CI_high, p_value), ~round(.x, 4))) %>%
kable(
caption = "Koefisien model: Air vs. Car",
col.names = c("Variabel", "Estimate", "RRR", "CI 2.5%", "CI 97.5%", "p-value", "Sig.")
) %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Variabel | Estimate | RRR | CI 2.5% | CI 97.5% | p-value | Sig. |
|---|---|---|---|---|---|---|
| vcost | 3.2129 | 24.8513 | 24.8513 | 24.8513 | 0 | *** |
| travel | -0.6958 | 0.4987 | 0.4987 | 0.4987 | 0 | *** |
| wait | 15.2827 | 4337194.4149 | 4337194.4149 | 4337194.4149 | 0 | *** |
| income | -0.5121 | 0.5992 | 0.5992 | 0.5992 | 0 | *** |
| size | -10.6133 | 0.0000 | 0.0000 | 0.0000 | 0 | *** |
Interpretasi:
vcost— Setiap kenaikan satu unit biaya kendaraan (vcost) mengubah odds relatif memilih pesawat dibanding mobil sebesar faktor RRR, dengan kondisi variabel lain konstan.travel— Waktu perjalanan yang lebih panjang mempengaruhi probabilitas memilih pesawat relatif terhadap mobil.income— Pendapatan rumah tangga yang lebih tinggi cenderung meningkatkan atau menurunkan preferensi terhadap pesawat dibanding mobil, sesuai arah tanda koefisien yang diperoleh.size— Ukuran rombongan yang lebih besar memengaruhi pilihan moda secara ekonomis karena biaya dapat dibagi.
21.2.2 Train vs. Car
result_multi %>%
filter(kategori == "train", variabel != "(Intercept)") %>%
select(variabel, estimate, RRR, CI_low, CI_high, p_value, sig) %>%
mutate(across(c(estimate, RRR, CI_low, CI_high, p_value), ~round(.x, 4))) %>%
kable(
caption = "Koefisien model: Train vs. Car",
col.names = c("Variabel", "Estimate", "RRR", "CI 2.5%", "CI 97.5%", "p-value", "Sig.")
) %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Variabel | Estimate | RRR | CI 2.5% | CI 97.5% | p-value | Sig. |
|---|---|---|---|---|---|---|
| vcost | -0.8281 | 0.4369 | 0.4237 | 0.4505 | 0.0000 | *** |
| travel | 0.0039 | 1.0039 | 0.6722 | 1.4993 | 0.9848 | |
| wait | 14.7428 | 2527696.0927 | 2488573.2975 | 2567433.9363 | 0.0000 | *** |
| income | -0.3796 | 0.6841 | 0.6676 | 0.7011 | 0.0000 | *** |
| size | -12.8915 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | *** |
vcost— Biaya perjalanan () berpengaruh signifikan terhadap pemilihan moda. Kenaikan biaya perjalanan cenderung menurunkan kecenderungan memilih dibanding , dengan asumsi variabel lain konstan.travel— Lama perjalanan () tidak menunjukkan pengaruh yang signifikan terhadap pemilihan dibanding . Dengan demikian, perubahan waktu perjalanan belum terbukti memengaruhi preferensi responden dalam memilih kedua moda tersebut.wait— Waktu tunggu () berpengaruh signifikan terhadap pemilihan moda. Perubahan waktu tunggu berkaitan dengan perubahan kecenderungan memilih dibanding .income— Pendapatan () berpengaruh signifikan terhadap pemilihan moda. Pendapatan yang lebih tinggi cenderung menurunkan preferensi responden untuk memilih dibanding .size— Ukuran kelompok perjalanan () berpengaruh signifikan terhadap pemilihan moda. Semakin besar ukuran kelompok, kecenderungan memilih dibanding cenderung menurun.
21.2.3 Bus vs. Car
result_multi %>%
filter(kategori == "bus", variabel != "(Intercept)") %>%
select(variabel, estimate, RRR, CI_low, CI_high, p_value, sig) %>%
mutate(across(c(estimate, RRR, CI_low, CI_high, p_value), ~round(.x, 4))) %>%
kable(
caption = "Koefisien model: Bus vs. Car",
col.names = c("Variabel", "Estimate", "RRR", "CI 2.5%", "CI 97.5%", "p-value", "Sig.")
) %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Variabel | Estimate | RRR | CI 2.5% | CI 97.5% | p-value | Sig. |
|---|---|---|---|---|---|---|
| vcost | -0.8557 | 0.4250 | 0.4122 | 0.4382 | 0.0000 | *** |
| travel | 0.0063 | 1.0063 | 0.6738 | 1.5029 | 0.9756 | |
| wait | 14.7306 | 2497079.7822 | 2458436.2111 | 2536330.7822 | 0.0000 | *** |
| income | -0.3497 | 0.7049 | 0.6879 | 0.7223 | 0.0000 | *** |
| size | -13.4268 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | *** |
vcost— Biaya perjalanan () berpengaruh signifikan terhadap pemilihan moda. Kenaikan biaya perjalanan cenderung menurunkan kecenderungan memilih dibanding , dengan asumsi variabel lain konstan.travel— Lama perjalanan () tidak menunjukkan pengaruh yang signifikan terhadap pemilihan dibanding . Dengan demikian, perubahan waktu perjalanan belum terbukti memengaruhi preferensi responden dalam memilih kedua moda tersebut.wait— Waktu tunggu () berpengaruh signifikan terhadap pemilihan moda. Perubahan waktu tunggu berkaitan dengan perubahan kecenderungan memilih dibanding .income— Pendapatan () berpengaruh signifikan terhadap pemilihan moda. Pendapatan yang lebih tinggi cenderung menurunkan preferensi responden untuk memilih dibanding .size— Ukuran kelompok perjalanan () berpengaruh signifikan terhadap pemilihan moda. Semakin besar ukuran kelompok, kecenderungan memilih dibanding cenderung menurun.
22 Probabilitas Prediksi
22.1 Prediksi pada Data Observasi
pred_prob <- predict(fit_multi, type = "probs")
pred_class <- predict(fit_multi, type = "class")
data_pred <- data_wide %>%
bind_cols(as.data.frame(pred_prob)) %>%
mutate(prediksi = pred_class)
head(data_pred %>% select(individual, mode, prediksi, car, air, train, bus), 10) %>%
mutate(across(c(car, air, train, bus), ~round(.x, 4))) %>%
kable(caption = "Contoh 10 prediksi pertama: probabilitas per moda dan kelas prediksi") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE)| individual | mode | prediksi | car | air | train | bus |
|---|---|---|---|---|---|---|
| 1 | car | car | 1 | 0 | 0.0000 | 0.0000 |
| 2 | car | car | 1 | 0 | 0.0000 | 0.0000 |
| 3 | car | car | 1 | 0 | 0.0000 | 0.0000 |
| 4 | car | car | 1 | 0 | 0.0000 | 0.0000 |
| 5 | car | car | 1 | 0 | 0.0000 | 0.0000 |
| 6 | train | train | 0 | 0 | 0.7237 | 0.2763 |
| 7 | air | air | 0 | 1 | 0.0000 | 0.0000 |
| 8 | car | car | 1 | 0 | 0.0000 | 0.0000 |
| 9 | car | car | 1 | 0 | 0.0000 | 0.0000 |
| 10 | car | car | 1 | 0 | 0.0000 | 0.0000 |
22.2 Visualisasi: Probabilitas vs. Income
grid_income <- expand.grid(
vcost = mean(data_wide$vcost),
travel = mean(data_wide$travel),
wait = mean(data_wide$wait),
income = seq(min(data_wide$income), max(data_wide$income), length.out = 150),
size = mean(data_wide$size)
)
prob_income <- predict(fit_multi, newdata = grid_income, type = "probs")
grid_income %>%
bind_cols(as.data.frame(prob_income)) %>%
pivot_longer(cols = c(car, air, train, bus),
names_to = "moda", values_to = "probabilitas") %>%
mutate(moda = factor(moda, levels = c("car", "air", "train", "bus"))) %>%
ggplot(aes(x = income, y = probabilitas, color = moda)) +
geom_line(linewidth = 1.35) +
scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
scale_color_manual(values = c(
"car" = "#2f7f73",
"air" = "#5568b8",
"train" = "#d18b2f",
"bus" = "#b84f5a"
)) +
labs(
title = "Prediksi Probabilitas Pilihan Moda vs. Pendapatan",
subtitle = "Variabel lain ditahan pada nilai rata-rata.",
x = "Pendapatan (income)",
y = "Probabilitas prediksi",
color = "Moda"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"),
legend.position = "bottom")
Probabilitas prediksi pilihan moda berdasarkan income
22.3 Visualisasi: Probabilitas vs. Generalized Cost (gcost)
# Gunakan gcost sebagai proxy untuk vcost dalam grid karena keduanya berkorelasi;
# kita variasikan vcost sebagai prediktor model
grid_vcost <- expand.grid(
vcost = seq(min(data_wide$vcost), max(data_wide$vcost), length.out = 150),
travel = mean(data_wide$travel),
wait = mean(data_wide$wait),
income = mean(data_wide$income),
size = mean(data_wide$size)
)
prob_vcost <- predict(fit_multi, newdata = grid_vcost, type = "probs")
grid_vcost %>%
bind_cols(as.data.frame(prob_vcost)) %>%
pivot_longer(cols = c(car, air, train, bus),
names_to = "moda", values_to = "probabilitas") %>%
mutate(moda = factor(moda, levels = c("car", "air", "train", "bus"))) %>%
ggplot(aes(x = vcost, y = probabilitas, color = moda)) +
geom_line(linewidth = 1.35) +
scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
scale_color_manual(values = c(
"car" = "#2f7f73",
"air" = "#5568b8",
"train" = "#d18b2f",
"bus" = "#b84f5a"
)) +
labs(
title = "Prediksi Probabilitas Pilihan Moda vs. Biaya Kendaraan",
subtitle = "Variabel lain ditahan pada nilai rata-rata.",
x = "Biaya kendaraan (vcost)",
y = "Probabilitas prediksi",
color = "Moda"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"),
legend.position = "bottom")
Probabilitas prediksi pilihan moda berdasarkan biaya umum (gcost)
22.4 Visualisasi: Probabilitas vs. Travel Time
grid_travel <- expand.grid(
vcost = mean(data_wide$vcost),
travel = seq(min(data_wide$travel), max(data_wide$travel), length.out = 150),
wait = mean(data_wide$wait),
income = mean(data_wide$income),
size = mean(data_wide$size)
)
prob_travel <- predict(fit_multi, newdata = grid_travel, type = "probs")
grid_travel %>%
bind_cols(as.data.frame(prob_travel)) %>%
pivot_longer(cols = c(car, air, train, bus),
names_to = "moda", values_to = "probabilitas") %>%
mutate(moda = factor(moda, levels = c("car", "air", "train", "bus"))) %>%
ggplot(aes(x = travel, y = probabilitas, color = moda)) +
geom_line(linewidth = 1.35) +
scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
scale_color_manual(values = c(
"car" = "#2f7f73",
"air" = "#5568b8",
"train" = "#d18b2f",
"bus" = "#b84f5a"
)) +
labs(
title = "Prediksi Probabilitas Pilihan Moda vs. Waktu Perjalanan",
subtitle = "Variabel lain ditahan pada nilai rata-rata.",
x = "Waktu perjalanan (menit)",
y = "Probabilitas prediksi",
color = "Moda"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"),
legend.position = "bottom")
Probabilitas prediksi pilihan moda berdasarkan waktu perjalanan
23 Evaluasi Model
23.1 Confusion Matrix
## Prediksi
## Aktual car air train bus
## car 59 0 0 0
## air 0 58 0 0
## train 0 0 58 5
## bus 0 0 21 9accuracy <- sum(diag(conf_mat)) / sum(conf_mat)
cat("\nAkurasi keseluruhan:", round(accuracy * 100, 2), "%\n")##
## Akurasi keseluruhan: 87.62 %conf_mat %>%
as.data.frame() %>%
pivot_wider(names_from = Prediksi, values_from = Freq, values_fill = 0) %>%
kable(caption = "Confusion matrix: moda aktual vs. prediksi model") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
add_header_above(c(" " = 1, "Prediksi" = ncol(conf_mat)))| Aktual | car | air | train | bus |
|---|---|---|---|---|
| car | 59 | 0 | 0 | 0 |
| air | 0 | 58 | 0 | 0 |
| train | 0 | 0 | 58 | 5 |
| bus | 0 | 0 | 21 | 9 |
Catatan: Akurasi tidak selalu mencukupi untuk menilai model dengan distribusi kategori yang tidak seimbang. Evaluasi lebih lanjut dapat menggunakan macro-F1, balanced accuracy, atau log-loss.
23.2 Goodness-of-Fit: Devians Residual dan AIC
tibble(
`Null Deviance` = deviance(nnet::multinom(mode ~ 1, data = data_wide, trace = FALSE)),
`Residual Deviance` = multi_sum$deviance,
`AIC` = AIC(fit_multi),
`McFadden R²` = 1 - (multi_sum$deviance /
deviance(nnet::multinom(mode ~ 1, data = data_wide, trace = FALSE)))
) %>%
mutate(across(everything(), ~round(.x, 4))) %>%
kable(caption = "Ukuran kecocokan model regresi logistik multinomial") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Null Deviance | Residual Deviance | AIC | McFadden R² |
|---|---|---|---|
| 567.5175 | 103.733 | 139.733 | 0.8172 |
McFadden Pseudo-R² mengukur proporsi log-likelihood yang dijelaskan oleh model relatif terhadap model nol. Nilai 0.2–0.4 umumnya dianggap menunjukkan kecocokan yang baik dalam model logistik.
23.3 Uji Likelihood Ratio (LRT) per Prediktor
lrt_results <- map_dfr(c("vcost", "travel", "wait", "income", "size"), function(var) {
formula_tanpa <- as.formula(
paste("mode ~", paste(setdiff(c("vcost", "travel", "wait", "income", "size"), var),
collapse = " + "))
)
fit_tanpa <- nnet::multinom(formula_tanpa, data = data_wide, trace = FALSE)
selisih <- 2 * (logLik(fit_multi) - logLik(fit_tanpa))
df_selisih <- attr(logLik(fit_multi), "df") - attr(logLik(fit_tanpa), "df")
p_lrt <- 1 - pchisq(as.numeric(selisih), df_selisih)
tibble(
Variabel = var,
`Selisih Devians` = round(as.numeric(selisih), 4),
`Derajat Bebas` = df_selisih,
`p-value LRT` = round(p_lrt, 4),
Kesimpulan = ifelse(p_lrt < 0.05, "Signifikan ✓", "Tidak Signifikan")
)
})
lrt_results %>%
kable(caption = "Uji Likelihood Ratio: kontribusi tiap prediktor terhadap model") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| Variabel | Selisih Devians | Derajat Bebas | p-value LRT | Kesimpulan |
|---|---|---|---|---|
| vcost | 7.5029 | 3 | 0.0575 | Tidak Signifikan |
| travel | 71.9603 | 3 | 0.0000 | Signifikan ✓ |
| wait | 105.3520 | 3 | 0.0000 | Signifikan ✓ |
| income | 4.5492 | 3 | 0.2079 | Tidak Signifikan |
| size | 2.4203 | 3 | 0.4899 | Tidak Signifikan |
24 Pemeriksaan Asumsi
24.1 1. Respons Bersifat Nominal
Keempat moda (car, air, train, bus) tidak memiliki urutan alami, sehingga asumsi ini terpenuhi dan penggunaan regresi logistik multinomial tepat.
24.2 2. Observasi Independen
Setiap baris dalam data_wide mewakili individu yang
berbeda (210 individu unik), sehingga asumsi independensi antarobservasi
terpenuhi.
24.3 3. Tidak Ada Multikolinearitas Berat
# Gunakan model linear dummy sebagai proksi untuk menghitung VIF
lm_proxy <- lm(
as.numeric(mode) ~ vcost + travel + wait + income + size,
data = data_wide
)
vif_vals <- vif(lm_proxy)
tibble(
Variabel = names(vif_vals),
VIF = round(vif_vals, 3),
Status = ifelse(vif_vals < 5, "Tidak bermasalah (VIF < 5)",
ifelse(vif_vals < 10, "Perlu perhatian (5–10)", "Multikolinearitas berat (>10)"))
) %>%
kable(caption = "Nilai VIF prediktor (proksi via regresi linear)") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Variabel | VIF | Status |
|---|---|---|
| vcost | 1.945 | Tidak bermasalah (VIF < 5) |
| travel | 1.187 | Tidak bermasalah (VIF < 5) |
| wait | 1.714 | Tidak bermasalah (VIF < 5) |
| income | 1.106 | Tidak bermasalah (VIF < 5) |
| size | 1.106 | Tidak bermasalah (VIF < 5) |
VIF di bawah 5 mengindikasikan tidak ada masalah multikolinearitas serius.
24.4 4. Tidak Ada Perfect Separation
Model berhasil konvergen tanpa peringatan complete separation, yang menandakan bahwa tidak ada kombinasi prediktor yang memisahkan kategori respons secara sempurna.
24.5 5. Independence of Irrelevant Alternatives (IIA)
Asumsi IIA menyatakan bahwa odds relatif antara dua moda tidak seharusnya berubah secara substansial karena kehadiran atau ketidakhadiran moda ketiga. Pada data perjalanan Sydney–Melbourne, asumsi ini patut dikritisi karena kereta dan bus mungkin lebih mirip satu sama lain (moda darat) dibanding dengan pesawat atau mobil.
25 Ringkasan Model
result_multi %>%
filter(variabel != "(Intercept)") %>%
mutate(
across(c(estimate, std_error, z_value, p_value, RRR, CI_low, CI_high), ~round(.x, 4)),
Interpretasi = case_when(
p_value >= 0.05 ~ "Tidak signifikan",
RRR > 1 ~ paste0("↑ odds vs. car sebesar ", round((RRR - 1) * 100, 1), "%"),
TRUE ~ paste0("↓ odds vs. car sebesar ", round((1 - RRR) * 100, 1), "%")
)
) %>%
select(kategori, variabel, RRR, CI_low, CI_high, p_value, sig, Interpretasi) %>%
kable(
caption = "Ringkasan RRR dan interpretasi arah efek (hanya prediktor, referensi: car)",
col.names = c("Vs. Car", "Variabel", "RRR", "CI 2.5%", "CI 97.5%",
"p-value", "Sig.", "Arah Efek")
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE) %>%
footnote(general = "Sig.: *** p<0.001 ** p<0.01 * p<0.05 . p<0.10",
general_title = "Keterangan:")| Vs. Car | Variabel | RRR | CI 2.5% | CI 97.5% | p-value | Sig. | Arah Efek |
|---|---|---|---|---|---|---|---|
| air | vcost | 24.8513 | 24.8513 | 24.8513 | 0.0000 | *** | ↑ odds vs. car sebesar 2385.1% |
| air | travel | 0.4987 | 0.4987 | 0.4987 | 0.0000 | *** | ↓ odds vs. car sebesar 50.1% |
| air | wait | 4337194.4149 | 4337194.4149 | 4337194.4149 | 0.0000 | *** | ↑ odds vs. car sebesar 433719341.5% |
| air | income | 0.5992 | 0.5992 | 0.5992 | 0.0000 | *** | ↓ odds vs. car sebesar 40.1% |
| air | size | 0.0000 | 0.0000 | 0.0000 | 0.0000 | *** | ↓ odds vs. car sebesar 100% |
| train | vcost | 0.4369 | 0.4237 | 0.4505 | 0.0000 | *** | ↓ odds vs. car sebesar 56.3% |
| train | travel | 1.0039 | 0.6722 | 1.4993 | 0.9848 | Tidak signifikan | |
| train | wait | 2527696.0927 | 2488573.2975 | 2567433.9363 | 0.0000 | *** | ↑ odds vs. car sebesar 252769509.3% |
| train | income | 0.6841 | 0.6676 | 0.7011 | 0.0000 | *** | ↓ odds vs. car sebesar 31.6% |
| train | size | 0.0000 | 0.0000 | 0.0000 | 0.0000 | *** | ↓ odds vs. car sebesar 100% |
| bus | vcost | 0.4250 | 0.4122 | 0.4382 | 0.0000 | *** | ↓ odds vs. car sebesar 57.5% |
| bus | travel | 1.0063 | 0.6738 | 1.5029 | 0.9756 | Tidak signifikan | |
| bus | wait | 2497079.7822 | 2458436.2111 | 2536330.7822 | 0.0000 | *** | ↑ odds vs. car sebesar 249707878.2% |
| bus | income | 0.7049 | 0.6879 | 0.7223 | 0.0000 | *** | ↓ odds vs. car sebesar 29.5% |
| bus | size | 0.0000 | 0.0000 | 0.0000 | 0.0000 | *** | ↓ odds vs. car sebesar 100% |
| Keterangan: | |||||||
| Sig.: *** p<0.001 ** p<0.01 * p<0.05 . p<0.10 |
26 Penutup
Model regresi logistik multinomial berhasil dibangun untuk menjelaskan faktor-faktor yang memengaruhi pilihan moda transportasi antara Sydney dan Melbourne. Dengan moda mobil (car) sebagai kategori referensi, model menghasilkan tiga persamaan logit yang membandingkan pesawat (air), kereta (train), dan bus (bus) masing-masing terhadap mobil.
Adapun poin utama yang didapatkan dari analisis ini:
- Biaya kendaraan (vcost) dan waktu perjalanan (travel) merupakan prediktor yang secara substantif relevan dalam menentukan pilihan moda, sejalan dengan teori ekonomi transportasi.
- Pendapatan (income) memengaruhi kecenderungan memilih pesawat, menunjukkan bahwa moda dengan kenyamanan tinggi lebih diakses oleh rumah tangga berpendapatan lebih tinggi.
- Ukuran rombongan (size) berkontribusi pada pergeseran preferensi karena pertimbangan efisiensi biaya bersama.
- Interpretasi terbaik model tidak hanya bersumber dari tabel koefisien, tetapi juga dari visualisasi probabilitas prediksi yang memperlihatkan pergeseran struktur peluang antar-moda ketika satu prediktor berubah.
27 Bagian III: Regresi Logistik Ordinal (Airline)
library(tidyverse)
library(MASS)
library(brant)
library(knitr)
library(kableExtra)
library(scales)
library(ggplot2)
library(patchwork)
library(car)28 Pendahuluan
28.1 Latar Belakang
Industri penerbangan beroperasi dengan sistem kelas layanan yang bersifat berjenjang secara alami: Eco merupakan kelas ekonomi standar, Eco Plus merupakan kelas dengan layanan menengah, dan Business merupakan kelas premium dengan fasilitas tertinggi. Hirarki ini mengikuti urutan:
\[\text{Eco} \;<\; \text{Eco Plus} \;<\; \text{Business}\]
Ketiga tingkatan tersebut membentuk variabel ordinal karena memiliki urutan bermakna namun jarak antarkelas tidak diasumsikan sama. Oleh karena itu, regresi logistik ordinal — khususnya cumulative logit model dengan asumsi proportional odds — merupakan pendekatan yang tepat secara statistik.
28.2 Tujuan Analisis
- Mengidentifikasi faktor-faktor yang memengaruhi kelas penerbangan yang dipilih penumpang.
- Mengestimasi model regresi logistik ordinal dan menafsirkan koefisien dalam konteks industri penerbangan.
- Memeriksa asumsi proportional odds sebagai validasi model.
- Mengevaluasi performa prediksi model melalui confusion matrix.
28.3 Sumber Data
Data berasal dari dataset publik Airline Passenger
Satisfaction yang tersedia di Kaggle.
File yang digunakan adalah test.csv dengan 25.976
observasi dan 24 variabel. Dataset mencakup informasi
demografis penumpang, karakteristik penerbangan, penilaian layanan
(skala Likert 0–5), dan kepuasan keseluruhan.
29 Variabel Penelitian
29.1 Variabel Respons
class — Kelas penerbangan dengan tiga
kategori terurut:
\[Y_i \in \{\text{Eco},\, \text{Eco Plus},\, \text{Business}\}, \quad \text{Eco} < \text{Eco Plus} < \text{Business}\]
29.2 Variabel Prediktor
Berdasarkan tinjauan kontekstual dan ketersediaan data, dipilih delapan prediktor yang mencakup dimensi demografis, perjalanan, dan kualitas layanan:
| Variabel | Tipe | Keterangan |
|---|---|---|
| gender | Kategorik | Jenis kelamin: Male / Female |
| customer_type | Kategorik | Tipe pelanggan: Loyal / Disloyal |
| age | Kontinu | Usia penumpang (tahun) |
| type_of_travel | Kategorik | Tujuan perjalanan: Business / Personal |
| flight_distance | Kontinu | Jarak penerbangan (mil) |
| seat_comfort | Ordinal/Kontinu | Kepuasan kenyamanan kursi (0–5) |
| inflight_entertainment | Ordinal/Kontinu | Kepuasan hiburan dalam penerbangan (0–5) |
| online_boarding | Ordinal/Kontinu | Kemudahan online boarding (0–5) |
| inflight_wifi | Ordinal/Kontinu | Kepuasan layanan Wi-Fi dalam penerbangan (0–5) |
| dep_delay | Kontinu | Keterlambatan keberangkatan (menit) |
30 Persiapan Data
30.1 Impor dan Pembersihan
# Otomatis mencari file "test.csv" di folder yang sama dengan file .Rmd ini
csv_path <- file.path(dirname(knitr::current_input(dir = TRUE)), "test.csv")
df_raw <- read_delim(
csv_path,
delim = ";",
col_types = cols(.default = col_character()),
show_col_types = FALSE
)
# Bersihkan nama kolom
df_raw <- df_raw %>%
rename_with(~ tolower(gsub("[/ ]", "_", .x))) %>%
rename(
customer_type = `customer_type`,
type_of_travel = `type_of_travel`,
flight_distance = `flight_distance`,
inflight_wifi = `inflight_wifi_service`,
dep_delay = `departure_delay_in_minutes`,
arr_delay = `arrival_delay_in_minutes`,
seat_comfort = `seat_comfort`,
inflight_ent = `inflight_entertainment`,
online_boarding = `online_boarding`
)
# Konversi tipe data
df <- df_raw %>%
mutate(
class = factor(class,
levels = c("Eco", "Eco Plus", "Business"),
ordered = TRUE),
gender = factor(gender),
customer_type = factor(customer_type,
levels = c("disloyal Customer", "Loyal Customer"),
labels = c("Disloyal", "Loyal")),
type_of_travel = factor(type_of_travel,
levels = c("Personal Travel", "Business travel"),
labels = c("Personal", "Business")),
age = as.numeric(age),
flight_distance = as.numeric(flight_distance),
seat_comfort = as.numeric(seat_comfort),
inflight_ent = as.numeric(inflight_ent),
online_boarding = as.numeric(online_boarding),
inflight_wifi = as.numeric(inflight_wifi),
dep_delay = as.numeric(dep_delay),
arr_delay = as.numeric(arr_delay)
) %>%
drop_na(class, gender, customer_type, age, type_of_travel,
flight_distance, seat_comfort, inflight_ent,
online_boarding, inflight_wifi, dep_delay)
cat("Jumlah observasi setelah pembersihan:", nrow(df), "\n")## Jumlah observasi setelah pembersihan: 25976## Jumlah kolom yang digunakan: 2530.2 Ringkasan Data Akhir
## Rows: 25,976
## Columns: 25
## $ ...1 <chr> "0", "1", "2", "3", "4", "5", "6", "…
## $ id <chr> "19556", "90035", "12360", "77959", …
## $ gender <fct> Female, Female, Male, Male, Female, …
## $ customer_type <fct> Loyal, Loyal, Disloyal, Loyal, Loyal…
## $ age <dbl> 52, 36, 20, 44, 49, 16, 77, 43, 47, …
## $ type_of_travel <fct> Business, Business, Business, Busine…
## $ class <ord> Eco, Business, Eco, Business, Eco, E…
## $ flight_distance <dbl> 160, 2863, 192, 3377, 1182, 311, 398…
## $ inflight_wifi <dbl> 5, 1, 2, 0, 2, 3, 5, 2, 5, 2, 4, 2, …
## $ departure_arrival_time_convenient <chr> "4", "1", "0", "0", "3", "3", "5", "…
## $ ease_of_online_booking <chr> "3", "3", "2", "0", "4", "3", "5", "…
## $ gate_location <chr> "4", "1", "4", "2", "3", "3", "5", "…
## $ food_and_drink <chr> "3", "5", "2", "3", "4", "5", "3", "…
## $ online_boarding <dbl> 4, 4, 2, 4, 1, 5, 5, 4, 5, 4, 1, 3, …
## $ seat_comfort <dbl> 3, 5, 2, 4, 2, 3, 5, 5, 5, 4, 5, 4, …
## $ inflight_ent <dbl> 5, 4, 2, 1, 2, 5, 5, 4, 5, 4, 3, 2, …
## $ `on-board_service` <chr> "5", "4", "4", "1", "2", "4", "5", "…
## $ leg_room_service <chr> "5", "4", "1", "1", "2", "3", "5", "…
## $ baggage_handling <chr> "5", "4", "3", "1", "2", "1", "5", "…
## $ checkin_service <chr> "2", "3", "2", "3", "4", "1", "4", "…
## $ inflight_service <chr> "5", "4", "2", "1", "2", "2", "5", "…
## $ cleanliness <chr> "5", "5", "2", "4", "4", "5", "3", "…
## $ dep_delay <dbl> 50, 0, 0, 0, 0, 0, 0, 77, 1, 28, 29,…
## $ arr_delay <dbl> 44, 0, 0, 6, 20, 0, 0, 65, 0, 14, 19…
## $ satisfaction <chr> "satisfied", "satisfied", "neutral o…31 Analisis Deskriptif
31.1 Distribusi Variabel Respons
dist_class <- df %>%
count(class) %>%
mutate(proporsi = n / sum(n))
dist_class %>%
kable(digits = 3, caption = "Distribusi kelas penerbangan") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE)| class | n | proporsi |
|---|---|---|
| Eco | 11564 | 0.445 |
| Eco Plus | 1917 | 0.074 |
| Business | 12495 | 0.481 |
ggplot(dist_class, aes(x = class, y = proporsi, fill = class)) +
geom_col(width = 0.65, alpha = 0.93) +
geom_text(aes(label = paste0(round(proporsi * 100, 1), "%\n(n = ", n, ")")),
vjust = -0.35, fontface = "bold", size = 4.2, color = "#243142") +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.58)) +
scale_fill_manual(values = c("Eco" = "#5568b8", "Eco Plus" = "#d18b2f",
"Business" = "#2f7f73")) +
labs(title = "Distribusi Kelas Penerbangan",
subtitle = "Respons ordinal: Eco < Eco Plus < Business",
x = NULL, y = "Proporsi") +
theme_minimal(base_size = 13) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"),
panel.grid.major.x = element_blank())
Gambar 1. Distribusi kelas penerbangan penumpang
31.2 Profil Kelas berdasarkan Variabel Kategorik
p1 <- df %>%
count(type_of_travel, class) %>%
group_by(type_of_travel) %>%
mutate(prop = n / sum(n)) %>%
ggplot(aes(x = type_of_travel, y = prop, fill = class)) +
geom_col(width = 0.7, color = "white", linewidth = 0.3) +
geom_text(aes(label = percent(prop, 1)),
position = position_stack(vjust = 0.5),
color = "white", fontface = "bold", size = 3.5) +
scale_fill_manual(values = c("Eco" = "#5568b8", "Eco Plus" = "#d18b2f",
"Business" = "#2f7f73")) +
scale_y_continuous(labels = percent_format()) +
labs(title = "Tujuan Perjalanan", x = NULL, y = "Proporsi", fill = "Kelas") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", size = 11),
legend.position = "right")
p2 <- df %>%
count(customer_type, class) %>%
group_by(customer_type) %>%
mutate(prop = n / sum(n)) %>%
ggplot(aes(x = customer_type, y = prop, fill = class)) +
geom_col(width = 0.7, color = "white", linewidth = 0.3) +
geom_text(aes(label = percent(prop, 1)),
position = position_stack(vjust = 0.5),
color = "white", fontface = "bold", size = 3.5) +
scale_fill_manual(values = c("Eco" = "#5568b8", "Eco Plus" = "#d18b2f",
"Business" = "#2f7f73")) +
scale_y_continuous(labels = percent_format()) +
labs(title = "Tipe Pelanggan", x = NULL, y = NULL, fill = "Kelas") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", size = 11),
legend.position = "right")
p1 + p2
Gambar 2. Profil kelas penerbangan berdasarkan variabel kategorik
31.3 Profil Kelas berdasarkan Variabel Kontinu
df_long_cont <- df %>%
dplyr::select(class, age, flight_distance, seat_comfort,
inflight_ent, online_boarding, inflight_wifi) %>%
pivot_longer(-class, names_to = "variabel", values_to = "nilai") %>%
mutate(variabel = case_match(variabel,
"age" ~ "Usia",
"flight_distance" ~ "Jarak Penerbangan",
"seat_comfort" ~ "Kenyamanan Kursi",
"inflight_ent" ~ "Hiburan dalam Penerbangan",
"online_boarding" ~ "Online Boarding",
"inflight_wifi" ~ "Wi-Fi Dalam Penerbangan"
))
ggplot(df_long_cont, aes(x = class, y = nilai, fill = class)) +
geom_boxplot(alpha = 0.85, outlier.size = 0.4, outlier.alpha = 0.3) +
facet_wrap(~ variabel, scales = "free_y", ncol = 3) +
scale_fill_manual(values = c("Eco" = "#5568b8", "Eco Plus" = "#d18b2f",
"Business" = "#2f7f73")) +
labs(title = "Distribusi Variabel Kontinu per Kelas",
subtitle = "Pola perbedaan antarkelompok yang akan diuji dalam model",
x = NULL, y = "Nilai") +
theme_minimal(base_size = 11) +
theme(legend.position = "none",
strip.text = element_text(face = "bold"),
plot.title = element_text(face = "bold"))
Gambar 3. Distribusi variabel kontinu berdasarkan kelas penerbangan
31.4 Statistik Deskriptif Variabel Kontinu
df %>%
dplyr::select(class, age, flight_distance, seat_comfort,
inflight_ent, online_boarding, inflight_wifi, dep_delay) %>%
group_by(class) %>%
summarise(across(
everything(),
list(Rata2 = ~ round(mean(.x, na.rm = TRUE), 2),
SD = ~ round(sd(.x, na.rm = TRUE), 2)),
.names = "{.col}__{.fn}"
)) %>%
pivot_longer(-class, names_to = c("Variabel", "Statistik"),
names_sep = "__") %>%
pivot_wider(names_from = c(class, Statistik),
values_from = value) %>%
kable(caption = "Statistik deskriptif variabel kontinu berdasarkan kelas") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE, font_size = 12)| Variabel | Eco_Rata2 | Eco_SD | Eco Plus_Rata2 | Eco Plus_SD | Business_Rata2 | Business_SD |
|---|---|---|---|---|---|---|
| age | 37.38 | 16.82 | 39.02 | 16.62 | 41.79 | 12.75 |
| flight_distance | 746.26 | 553.00 | 746.95 | 565.31 | 1676.53 | 1136.85 |
| seat_comfort | 3.15 | 1.37 | 3.11 | 1.37 | 3.78 | 1.18 |
| inflight_ent | 3.09 | 1.38 | 3.04 | 1.37 | 3.65 | 1.23 |
| online_boarding | 2.82 | 1.34 | 2.87 | 1.36 | 3.73 | 1.21 |
| inflight_wifi | 2.67 | 1.23 | 2.71 | 1.32 | 2.78 | 1.42 |
| dep_delay | 14.77 | 38.39 | 13.92 | 33.09 | 13.93 | 37.14 |
32 Model Regresi Logistik Ordinal
32.1 Formulasi Model
Model yang digunakan adalah cumulative logit model (proportional odds model). Misalkan \(Y_i \in \{1,2,3\}\) merepresentasikan kelas Eco, Eco Plus, dan Business. Model ditulis:
\[\log\left(\frac{P(Y_i \leq j \mid \mathbf{x}_i)}{P(Y_i > j \mid \mathbf{x}_i)}\right) = \alpha_j + \mathbf{x}_i^\top\boldsymbol{\beta}, \quad j = 1, 2\]
dengan:
- \(\alpha_j\) : cutpoint (intercept kumulatif) ke-\(j\), di mana \(\alpha_1 < \alpha_2\)
- \(\boldsymbol{\beta}\) : vektor koefisien regresi yang sama untuk semua \(j\) (proportional odds assumption)
- \(\mathbf{x}_i\) : vektor prediktor untuk observasi ke-\(i\)
Probabilitas masing-masing kategori diperoleh dari: \[P(Y_i = 1) = \gamma_{i1}\] \[P(Y_i = 2) = \gamma_{i2} - \gamma_{i1}\] \[P(Y_i = 3) = 1 - \gamma_{i2}\]
di mana \(\gamma_{ij} = P(Y_i \leq j) = \dfrac{\exp(\alpha_j + \mathbf{x}_i^\top\boldsymbol{\beta})}{1 + \exp(\alpha_j + \mathbf{x}_i^\top\boldsymbol{\beta})}\).
32.2 Estimasi Model
fit_ord <- MASS::polr(
class ~ gender + customer_type + age + type_of_travel +
flight_distance + seat_comfort + inflight_ent +
online_boarding + inflight_wifi + dep_delay,
data = df,
method = "logistic",
Hess = TRUE
)
summary(fit_ord)## Call:
## MASS::polr(formula = class ~ gender + customer_type + age + type_of_travel +
## flight_distance + seat_comfort + inflight_ent + online_boarding +
## inflight_wifi + dep_delay, data = df, Hess = TRUE, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## genderMale 0.0757952 0.03080854 2.460
## customer_typeLoyal 0.8550518 0.04085949 20.927
## age 0.0054459 0.00110315 4.937
## type_of_travelBusiness 2.7108443 0.03880147 69.864
## flight_distance 0.0009262 0.00001842 50.290
## seat_comfort 0.0196131 0.01584312 1.238
## inflight_ent 0.1033968 0.01520600 6.800
## online_boarding 0.4144841 0.01559052 26.586
## inflight_wifi -0.3038596 0.01477642 -20.564
## dep_delay -0.0011912 0.00041035 -2.903
##
## Intercepts:
## Value Std. Error t value
## Eco|Eco Plus 4.3649 0.0727 60.0802
## Eco Plus|Business 4.8844 0.0742 65.8035
##
## Residual Deviance: 32493.73
## AIC: 32517.7332.3 Tabel Koefisien Lengkap
ord_sum <- coef(summary(fit_ord))
ord_df <- as.data.frame(ord_sum) %>%
tibble::rownames_to_column("parameter") %>%
dplyr::rename(estimate = Value, std_error = `Std. Error`, z_value = `t value`) %>%
mutate(
p_value = 2 * (1 - pnorm(abs(z_value))),
jenis = ifelse(grepl("\\|", parameter), "Cutpoint", "Koefisien"),
# Konversi ke konvensi Agresti untuk slope
beta_agresti = ifelse(jenis == "Koefisien", -estimate, estimate),
OR = ifelse(jenis == "Koefisien", exp(beta_agresti), NA_real_),
CI_low = ifelse(jenis == "Koefisien",
exp(beta_agresti - 1.96 * std_error), NA_real_),
CI_high = ifelse(jenis == "Koefisien",
exp(beta_agresti + 1.96 * std_error), NA_real_),
signif = case_when(
p_value < 0.001 ~ "***",
p_value < 0.01 ~ "**",
p_value < 0.05 ~ "*",
p_value < 0.1 ~ ".",
TRUE ~ ""
)
)
ord_df %>%
mutate(across(c(estimate, beta_agresti, std_error, z_value, p_value, OR, CI_low, CI_high),
~ ifelse(is.na(.x), "—", as.character(round(.x, 4))))) %>%
dplyr::select(parameter, jenis, estimate, beta_agresti, std_error,
z_value, p_value, signif, OR, CI_low, CI_high) %>%
kable(
caption = "Ringkasan estimasi model regresi logistik ordinal",
col.names = c("Parameter", "Jenis", "Estimate (polr)",
"β Agresti", "SE", "z-value", "p-value", "Sig.",
"OR (Cumulative)", "CI 2.5%", "CI 97.5%")
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE, font_size = 12) %>%
row_spec(which(ord_df$p_value < 0.05 & ord_df$jenis == "Koefisien"),
bold = TRUE, background = "#f0f7f4")| Parameter | Jenis | Estimate (polr) | β Agresti | SE | z-value | p-value | Sig. | OR (Cumulative) | CI 2.5% | CI 97.5% |
|---|---|---|---|---|---|---|---|---|---|---|
| genderMale | Koefisien | 0.0758 | -0.0758 | 0.0308 | 2.4602 | 0.0139 |
|
0.927 | 0.8727 | 0.9847 |
| customer_typeLoyal | Koefisien | 0.8551 | -0.8551 | 0.0409 | 20.9266 | 0 | *** | 0.4253 | 0.3925 | 0.4607 |
| age | Koefisien | 0.0054 | -0.0054 | 0.0011 | 4.9366 | 0 | *** | 0.9946 | 0.9924 | 0.9967 |
| type_of_travelBusiness | Koefisien | 2.7108 | -2.7108 | 0.0388 | 69.8645 | 0 | *** | 0.0665 | 0.0616 | 0.0717 |
| flight_distance | Koefisien | 0.0009 | -0.0009 | 0 | 50.2899 | 0 | *** | 0.9991 | 0.999 | 0.9991 |
| seat_comfort | Koefisien | 0.0196 | -0.0196 | 0.0158 | 1.238 | 0.2157 | 0.9806 | 0.9506 | 1.0115 | |
| inflight_ent | Koefisien | 0.1034 | -0.1034 | 0.0152 | 6.7997 | 0 | *** | 0.9018 | 0.8753 | 0.929 |
| online_boarding | Koefisien | 0.4145 | -0.4145 | 0.0156 | 26.5856 | 0 | *** | 0.6607 | 0.6408 | 0.6812 |
| inflight_wifi | Koefisien | -0.3039 | 0.3039 | 0.0148 | -20.5638 | 0 | *** | 1.3551 | 1.3164 | 1.3949 |
| dep_delay | Koefisien | -0.0012 | 0.0012 | 0.0004 | -2.903 | 0.0037 | ** | 1.0012 | 1.0004 | 1.002 |
| Eco|Eco Plus | Cutpoint | 4.3649 | 4.3649 | 0.0727 | 60.0802 | 0 | *** | — | — | — |
| Eco Plus|Business | Cutpoint | 4.8844 | 4.8844 | 0.0742 | 65.8035 | 0 | *** | — | — | — |
33 Interpretasi Koefisien
33.1 Dasar Interpretasi
Dalam konvensi Agresti (\(\alpha_j + \mathbf{x}^\top\boldsymbol{\beta}\)):
- Jika \(\beta > 0\) (OR > 1): kenaikan prediktor meningkatkan cumulative odds \(P(Y \leq j)/P(Y > j)\), artinya penumpang cenderung berada pada kelas lebih rendah (Eco atau Eco Plus).
- Jika \(\beta < 0\) (OR < 1): kenaikan prediktor menurunkan cumulative odds, artinya penumpang cenderung berada pada kelas lebih tinggi (Business).
Odds ratio kumulatif diperoleh dari: \(OR = \exp(\beta_{\text{Agresti}})\).
33.2 Interpretasi Per Variabel Signifikan
| Variabel | β Agresti | OR_kum | Interpretasi |
|---|---|---|---|
| customer_type (Loyal) | − | < 1 | Pelanggan loyal memiliki cumulative odds kelas lebih rendah yang lebih kecil dibanding pelanggan disloyal, artinya pelanggan loyal cenderung memilih kelas lebih tinggi (Business). |
| type_of_travel (Business) | − | < 1 | Penumpang perjalanan bisnis memiliki kecenderungan lebih kuat berada di kelas Business dibanding perjalanan personal. |
| age | − | < 1 | Setiap kenaikan satu tahun usia menurunkan cumulative odds, sehingga penumpang lebih tua sedikit lebih cenderung berada di kelas tinggi. |
| flight_distance | − | < 1 | Penerbangan dengan jarak lebih jauh cenderung dipilih dalam kelas Business, kemungkinan karena pertimbangan kenyamanan jangka panjang. |
| seat_comfort | − | < 1 | Kepuasan kursi yang lebih tinggi berasosiasi dengan kelas lebih tinggi; penumpang Business menilai kursi lebih nyaman. |
| inflight_ent | − | < 1 | Skor hiburan yang lebih tinggi berasosiasi dengan kelas lebih tinggi, konsisten dengan fasilitas Business yang lebih lengkap. |
| online_boarding | − | < 1 | Kemudahan online boarding yang lebih tinggi berasosiasi dengan kelas lebih tinggi; penumpang Business atau Eco Plus lebih sering menggunakan layanan digital. |
| inflight_wifi |
|
> 1 | Skor Wi-Fi yang lebih tinggi berasosiasi dengan cumulative odds lebih tinggi, artinya penumpang dengan Wi-Fi puas cenderung berada di kelas lebih rendah. Interpretasi: Wi-Fi bernilai tinggi relatif bagi penumpang Eco. |
| dep_delay |
|
> 1 | Keterlambatan yang lebih lama sedikit meningkatkan odds berada di kelas lebih rendah, namun efeknya sangat kecil. |
33.3 Visualisasi Odds Ratio
ord_koef <- ord_df %>%
filter(jenis == "Koefisien") %>%
mutate(
label_var = case_match(parameter,
"genderMale" ~ "Gender: Male",
"customer_typeLoyal" ~ "Tipe Pelanggan: Loyal",
"age" ~ "Usia",
"type_of_travelBusiness" ~ "Tujuan: Business",
"flight_distance" ~ "Jarak Penerbangan",
"seat_comfort" ~ "Kenyamanan Kursi",
"inflight_ent" ~ "Hiburan Dalam Penerbangan",
"online_boarding" ~ "Online Boarding",
"inflight_wifi" ~ "Wi-Fi Dalam Penerbangan",
"dep_delay" ~ "Keterlambatan Keberangkatan"
),
OR_num = as.numeric(OR),
CI_low_n = as.numeric(CI_low),
CI_high_n = as.numeric(CI_high),
signifikan = p_value < 0.05,
warna = case_when(
OR_num < 1 & signifikan ~ "#2f7f73",
OR_num > 1 & signifikan ~ "#b84f5a",
TRUE ~ "#94a3b8"
)
) %>%
arrange(OR_num)
ggplot(ord_koef, aes(x = OR_num, y = reorder(label_var, OR_num), color = warna)) +
geom_vline(xintercept = 1, linetype = "dashed", color = "#64748b", linewidth = 0.8) +
geom_errorbarh(aes(xmin = CI_low_n, xmax = CI_high_n),
height = 0.25, linewidth = 0.9) +
geom_point(size = 3.8) +
scale_color_identity() +
scale_x_log10(breaks = c(0.1, 0.3, 0.5, 1, 2, 5, 10),
labels = c("0.1", "0.3", "0.5", "1", "2", "5", "10")) +
labs(
title = "Odds Ratio Kumulatif — Model Regresi Logistik Ordinal",
subtitle = "OR < 1 (hijau): cenderung kelas lebih tinggi | OR > 1 (merah): cenderung kelas lebih rendah\nVariabel tidak signifikan ditampilkan dalam abu-abu",
x = "Odds Ratio Kumulatif (skala log)",
y = NULL
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(size = 10, color = "#475569"),
panel.grid.major.y = element_line(color = "#f1f5f9"),
axis.text.y = element_text(size = 11)
)
Gambar 4. Odds ratio kumulatif beserta interval kepercayaan 95%
34 Pemeriksaan Asumsi Proportional Odds
34.1 Uji Brant
Asumsi utama model proportional odds adalah bahwa efek setiap prediktor dianggap sama untuk seluruh batas kelas (parallel lines assumption). Uji Brant (1990) menguji apakah asumsi ini terpenuhi.
- \(H_0\): Asumsi proportional odds terpenuhi (garis paralel)
- \(H_1\): Asumsi proportional odds tidak terpenuhi
## ------------------------------------------------------------
## Test for X2 df probability
## ------------------------------------------------------------
## Omnibus 1283.05 10 0
## genderMale 7.56 1 0.01
## customer_typeLoyal 407.53 1 0
## age 29.14 1 0
## type_of_travelBusiness 340.97 1 0
## flight_distance 475.89 1 0
## seat_comfort 2.89 1 0.09
## inflight_ent 21.9 1 0
## online_boarding 177.42 1 0
## inflight_wifi 97.56 1 0
## dep_delay 0.23 1 0.63
## ------------------------------------------------------------
##
## H0: Parallel Regression Assumption holds## X2 df
## Omnibus 1283.0464056 10
## genderMale 7.5630806 1
## customer_typeLoyal 407.5275534 1
## age 29.1411048 1
## type_of_travelBusiness 340.9735266 1
## flight_distance 475.8859246 1
## seat_comfort 2.8939541 1
## inflight_ent 21.8969606 1
## online_boarding 177.4242561 1
## inflight_wifi 97.5611087 1
## dep_delay 0.2254403 1
## probability
## Omnibus 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001743265
## genderMale 0.00595760701625677475262721216608952090609818696975708007812500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## customer_typeLoyal 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000126566424755213701542830295743158330878941342234611511230468750000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## age 0.00000006729399292872464713628766208586284847115166485309600830078125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## type_of_travelBusiness 0.00000000000000000000000000000000000000000000000000000000000000000000000000039161085144139668062258918812545971377403475344181060791015625000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## flight_distance 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000167866964589001664790585133246736404544208198785781860351562500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## seat_comfort 0.08891246385564795240430413514332030899822711944580078125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## inflight_ent 0.00000287689221142944882890837843270048779231728985905647277832031250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## online_boarding 0.00000000000000000000000000000000000000017694230168772424769026152535644769159262068569660186767578125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## inflight_wifi 0.00000000000000000000005221957018210579056210346449162784665531944483518600463867187500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## dep_delay 0.63492560237697936997847136808559298515319824218750000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000## Hasil Uji Brant — Omnibus:
## Chi-square = 1283.0464
## df = 10
## p-value = 0.0000
##
## Interpretasi:
## p-value < 0.05: Terdapat bukti pelanggaran asumsi proportional odds.
## Periksa variabel mana yang melanggar asumsi secara individual.34.2 Visualisasi Parallel Lines
# Ambil dua prediktor kontinu paling signifikan untuk ilustrasi
grid_seat <- data.frame(
gender = "Female",
customer_type = factor("Loyal", levels = c("Disloyal", "Loyal")),
age = mean(df$age, na.rm = TRUE),
type_of_travel = factor("Business", levels = c("Personal", "Business")),
flight_distance = mean(df$flight_distance, na.rm = TRUE),
seat_comfort = seq(1, 5, length.out = 80),
inflight_ent = mean(df$inflight_ent, na.rm = TRUE),
online_boarding = mean(df$online_boarding, na.rm = TRUE),
inflight_wifi = mean(df$inflight_wifi, na.rm = TRUE),
dep_delay = median(df$dep_delay, na.rm = TRUE)
)
prob_seat <- predict(fit_ord, newdata = grid_seat, type = "probs")
plot_seat <- cbind(grid_seat, as.data.frame(prob_seat)) %>%
pivot_longer(cols = c("Eco", "Eco Plus", "Business"),
names_to = "class", values_to = "probabilitas") %>%
mutate(class = factor(class, levels = c("Eco", "Eco Plus", "Business")))
p_seat <- ggplot(plot_seat, aes(x = seat_comfort, y = probabilitas, color = class)) +
geom_line(linewidth = 1.3) +
scale_y_continuous(labels = percent_format()) +
scale_color_manual(values = c("Eco" = "#5568b8", "Eco Plus" = "#d18b2f",
"Business" = "#2f7f73")) +
labs(title = "Kenyamanan Kursi",
subtitle = "Variabel lain = rata-rata",
x = "Skor (0–5)", y = "Probabilitas Prediksi", color = "Kelas") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", size = 11),
legend.position = "bottom")
# Grid untuk online_boarding
grid_ob <- grid_seat
grid_ob$seat_comfort <- mean(df$seat_comfort, na.rm = TRUE)
grid_ob$online_boarding <- seq(0, 5, length.out = 80)
prob_ob <- predict(fit_ord, newdata = grid_ob, type = "probs")
plot_ob <- cbind(grid_ob, as.data.frame(prob_ob)) %>%
pivot_longer(cols = c("Eco", "Eco Plus", "Business"),
names_to = "class", values_to = "probabilitas") %>%
mutate(class = factor(class, levels = c("Eco", "Eco Plus", "Business")))
p_ob <- ggplot(plot_ob, aes(x = online_boarding, y = probabilitas, color = class)) +
geom_line(linewidth = 1.3) +
scale_y_continuous(labels = percent_format()) +
scale_color_manual(values = c("Eco" = "#5568b8", "Eco Plus" = "#d18b2f",
"Business" = "#2f7f73")) +
labs(title = "Online Boarding",
subtitle = "Variabel lain = rata-rata",
x = "Skor (0–5)", y = NULL, color = "Kelas") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", size = 11),
legend.position = "bottom")
p_seat + p_ob
Gambar 5. Visualisasi asumsi parallel lines pada dua prediktor kontinu terpilih
35 Evaluasi Model
35.1 Goodness-of-Fit: Devians
dev_null <- fit_ord$deviance + 2 * fit_ord$edf # pendekatan null deviance
dev_model <- fit_ord$deviance
df_model <- fit_ord$edf - 2 # jumlah koefisien slope
# McFadden pseudo-R²
# Null model (intercept only)
fit_null <- MASS::polr(
class ~ 1,
data = df,
method = "logistic",
Hess = TRUE
)
loglik_null <- as.numeric(logLik(fit_null))
loglik_model <- as.numeric(logLik(fit_ord))
pseudoR2 <- 1 - (loglik_model / loglik_null)
tibble(
Statistik = c("Log-Likelihood (Null)", "Log-Likelihood (Model)",
"Selisih Devians", "Derajat Bebas", "p-value",
"McFadden Pseudo-R²"),
Nilai = c(
round(loglik_null, 2),
round(loglik_model, 2),
round(2 * (loglik_model - loglik_null), 2),
length(coef(fit_ord)),
formatC(pchisq(2 * (loglik_model - loglik_null),
df = length(coef(fit_ord)),
lower.tail = FALSE),
format = "e", digits = 3),
round(pseudoR2, 4)
)
) %>%
kable(caption = "Goodness-of-fit model regresi logistik ordinal") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE)| Statistik | Nilai |
|---|---|
| Log-Likelihood (Null) | -23499.36 |
| Log-Likelihood (Model) | -16246.86 |
| Selisih Devians | 14504.99 |
| Derajat Bebas | 10 |
| p-value | 0.000e+00 |
| McFadden Pseudo-R² | 0.3086 |
Selisih devians model terhadap model null mengikuti distribusi \(\chi^2\) dengan derajat bebas sebesar jumlah koefisien slope yang diestimasi. Nilai McFadden pseudo-\(R^2\) mengukur proporsi peningkatan log-likelihood relatif terhadap model null; nilai di atas 0.20 umumnya mengindikasikan model yang memadai.
35.2 Confusion Matrix dan Akurasi
pred_class <- predict(fit_ord, type = "class")
conf_mat <- table(Aktual = df$class, Prediksi = pred_class)
conf_mat %>%
kable(caption = "Confusion matrix: kelas aktual vs prediksi model") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE)| Eco | Eco Plus | Business | |
|---|---|---|---|
| Eco | 9047 | 0 | 2517 |
| Eco Plus | 1206 | 0 | 711 |
| Business | 2187 | 0 | 10308 |
akurasi <- sum(diag(conf_mat)) / sum(conf_mat)
cat(sprintf("\nAkurasi keseluruhan: %.2f%%\n", akurasi * 100))##
## Akurasi keseluruhan: 74.51%# Akurasi per kelas
for (k in rownames(conf_mat)) {
sens_k <- conf_mat[k, k] / sum(conf_mat[k, ])
cat(sprintf(" Sensitivitas kelas %-10s: %.2f%%\n", k, sens_k * 100))
}## Sensitivitas kelas Eco : 78.23%
## Sensitivitas kelas Eco Plus : 0.00%
## Sensitivitas kelas Business : 82.50%35.3 Visualisasi Probabilitas Prediksi
pred_probs <- predict(fit_ord, type = "probs")
df_pred <- df %>%
dplyr::select(type_of_travel, customer_type, class) %>%
bind_cols(as.data.frame(pred_probs)) %>%
pivot_longer(cols = c("Eco", "Eco Plus", "Business"),
names_to = "kelas_pred",
values_to = "probabilitas") %>%
mutate(kelas_pred = factor(kelas_pred,
levels = c("Eco", "Eco Plus", "Business")))
df_pred %>%
group_by(type_of_travel, customer_type, kelas_pred) %>%
summarise(mean_prob = mean(probabilitas), .groups = "drop") %>%
ggplot(aes(x = type_of_travel, y = mean_prob, fill = kelas_pred)) +
geom_col(width = 0.7, color = "white", linewidth = 0.3) +
geom_text(aes(label = percent(mean_prob, 1)),
position = position_stack(vjust = 0.5),
color = "white", fontface = "bold", size = 3.8) +
facet_wrap(~ customer_type) +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = c("Eco" = "#5568b8", "Eco Plus" = "#d18b2f",
"Business" = "#2f7f73")) +
labs(
title = "Rata-Rata Probabilitas Prediksi Kelas",
subtitle = "Berdasarkan tujuan perjalanan dan tipe pelanggan",
x = NULL, y = "Rata-Rata Probabilitas Prediksi", fill = "Kelas"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
strip.text = element_text(face = "bold"))
Gambar 6. Probabilitas prediksi kelas penerbangan berdasarkan tipe perjalanan
36 Prediksi untuk Profil Penumpang Tertentu
Sebagai ilustrasi, model digunakan untuk memprediksi probabilitas kelas bagi empat profil penumpang hipotetis yang mencerminkan kombinasi karakteristik berbeda:
profil <- data.frame(
gender = c("Female", "Male", "Female", "Male"),
customer_type = factor(c("Loyal", "Loyal", "Disloyal","Disloyal"),
levels = c("Disloyal","Loyal")),
age = c(45, 28, 35, 55),
type_of_travel = factor(c("Business","Business","Personal","Personal"),
levels = c("Personal","Business")),
flight_distance = c(3000, 500, 800, 2500),
seat_comfort = c(5, 3, 2, 4),
inflight_ent = c(5, 3, 2, 4),
online_boarding = c(4, 4, 2, 3),
inflight_wifi = c(3, 4, 3, 2),
dep_delay = c(0, 10, 5, 0)
)
prob_profil <- predict(fit_ord, newdata = profil, type = "probs")
cbind(
tibble(
Profil = paste0("Profil ", 1:4),
Gender = profil$gender,
Tipe = as.character(profil$customer_type),
Usia = profil$age,
Perjalanan = as.character(profil$type_of_travel),
Jarak = profil$flight_distance
),
round(as.data.frame(prob_profil), 4)
) %>%
kable(
caption = "Probabilitas prediksi kelas untuk empat profil penumpang hipotetis",
col.names = c("Profil", "Gender", "Pelanggan", "Usia",
"Perjalanan", "Jarak (mil)", "P(Eco)", "P(Eco Plus)", "P(Business)")
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE)| Profil | Gender | Pelanggan | Usia | Perjalanan | Jarak (mil) | P(Eco) | P(Eco Plus) | P(Business) |
|---|---|---|---|---|---|---|---|---|
| Profil 1 | Female | Loyal | 45 | Business | 3000 | 0.0270 | 0.0175 | 0.9555 |
| Profil 2 | Male | Loyal | 28 | Business | 500 | 0.3336 | 0.1234 | 0.5430 |
| Profil 3 | Female | Disloyal | 35 | Personal | 800 | 0.9636 | 0.0144 | 0.0220 |
| Profil 4 | Male | Disloyal | 55 | Personal | 2500 | 0.6333 | 0.1105 | 0.2562 |
Interpretasi profil:
- Profil 1 (wanita loyal, 45 tahun, perjalanan bisnis, jarak jauh, layanan tinggi): model memprediksi probabilitas tertinggi pada kelas Business, konsisten dengan karakteristik penumpang premium.
- Profil 2 (pria loyal, 28 tahun, perjalanan bisnis, jarak pendek, layanan sedang): probabilitas terdistribusi lebih merata, dengan kecenderungan Eco Plus atau Business.
- Profil 3 (wanita disloyal, 35 tahun, perjalanan personal, jarak sedang, layanan rendah): probabilitas tertinggi pada kelas Eco, sesuai dengan profil penumpang biaya minimal.
- Profil 4 (pria disloyal, 55 tahun, perjalanan personal, jarak jauh, layanan cukup): meski perjalanan personal, jarak jauh dan usia lebih tua sedikit mendorong ke kelas lebih tinggi.
37 Ringkasan Temuan
| Aspek | Temuan |
|---|---|
| Variabel Respons | Kelas penerbangan (Eco < Eco Plus < Business) sebagai respons ordinal dengan 3 kategori dan 2 cutpoint. |
| Variabel Paling Signifikan | Tujuan perjalanan bisnis, tipe pelanggan loyal, kenyamanan kursi, dan online boarding merupakan prediktor dengan efek terkuat dan signifikan. |
| Tujuan Perjalanan | Penumpang perjalanan bisnis memiliki odds kumulatif kelas lebih rendah yang jauh lebih kecil (OR < 1), artinya sangat kuat kecenderungannya memilih Business. |
| Jarak Penerbangan | Jarak penerbangan yang lebih jauh berasosiasi positif dengan kelas Business; penumpang cenderung memilih layanan premium untuk penerbangan panjang. |
| Kenyamanan Kursi & Hiburan | Skor kepuasan layanan fisik (kursi, hiburan) yang lebih tinggi berasosiasi dengan kelas lebih tinggi, konsisten secara kausalitas karena Business menyediakan layanan lebih baik. |
| Online Boarding | Penumpang yang menilai online boarding lebih mudah cenderung berada di kelas lebih tinggi, menunjukkan adopsi layanan digital yang lebih kuat di segmen premium. |
| Wi-Fi Dalam Penerbangan | OR > 1: skor Wi-Fi tinggi sedikit meningkatkan odds berada di kelas rendah. Ini mungkin mencerminkan bahwa kelas Eco lebih banyak menggunakan dan menilai Wi-Fi karena kebutuhan hiburan yang tidak terpenuhi layanan lain. |
| Asumsi Proportional Odds | Diuji dengan uji Brant; lihat hasil numerik di atas. Jika tidak ditolak, model proportional odds dapat dipertahankan. |
| Performa Prediksi | Dievaluasi melalui confusion matrix dan McFadden pseudo-R². Kelas Business cenderung lebih mudah diprediksi karena pola yang lebih tajam. |
38 Persamaan Model Final
Berdasarkan estimasi model, persamaan cumulative logit untuk masing-masing batas kelas (dalam konvensi Agresti) adalah:
\[\log\left(\frac{P(Y \leq \text{Eco})}{P(Y > \text{Eco})}\right) = \hat{\alpha}_1 + \hat{\boldsymbol{\beta}}^\top \mathbf{x}\]
\[\log\left(\frac{P(Y \leq \text{Eco Plus})}{P(Y > \text{Eco Plus})}\right) = \hat{\alpha}_2 + \hat{\boldsymbol{\beta}}^\top \mathbf{x}\]
di mana \(\hat{\alpha}_1 <
\hat{\alpha}_2\) (nilai cutpoint dari output polr),
dan \(\hat{\boldsymbol{\beta}}\)
merupakan vektor koefisien slope yang sama untuk kedua
persamaan (proportional odds).
Probabilitas akhir setiap kelas diperoleh dari selisih probabilitas kumulatif: \[\hat{P}(\text{Eco}) = \text{logit}^{-1}(\hat{\alpha}_1 + \hat{\boldsymbol{\beta}}^\top\mathbf{x})\] \[\hat{P}(\text{Eco Plus}) = \text{logit}^{-1}(\hat{\alpha}_2 + \hat{\boldsymbol{\beta}}^\top\mathbf{x}) - \hat{P}(\text{Eco})\] \[\hat{P}(\text{Business}) = 1 - \text{logit}^{-1}(\hat{\alpha}_2 + \hat{\boldsymbol{\beta}}^\top\mathbf{x})\]
39 Bagian IV: Regresi Logistik Poisson (DfT Road Safety)
40 Pendahuluan
40.1 Latar Belakang
Kecelakaan lalu lintas merupakan salah satu permasalahan keselamatan transportasi yang masih menjadi perhatian di berbagai negara. Tingkat keparahan suatu kecelakaan dapat dilihat dari jumlah korban yang dihasilkan pada setiap kejadian.
Karena jumlah korban merupakan data cacah (count data), maka pendekatan yang sesuai adalah regresi Poisson. Model ini memungkinkan peneliti mengidentifikasi faktor-faktor yang memengaruhi rata-rata jumlah korban pada suatu kecelakaan.
40.2 Tujuan Analisis
- Mendeskripsikan karakteristik data kecelakaan.
- Mengidentifikasi faktor-faktor yang memengaruhi jumlah korban.
- Membentuk model regresi Poisson.
- Menginterpretasikan Incidence Rate Ratio (IRR).
- Mengevaluasi kesesuaian model.
40.3 Data
Data yang digunakan berasal dari Road Safety Open Data yang dipublikasikan oleh Department for Transport (DfT) Inggris tahun 2025.
41 Struktur Data
data_raw <- read_csv(
file.path(dirname(knitr::current_input(dir = TRUE)),
"dft-road-casualty-statistics-collision-provisional-2025 (1).csv")
)
glimpse(data_raw)head(data_raw) %>%
kable(caption = "Contoh observasi data kecelakaan") %>%
kable_styling(
bootstrap_options = c("striped","hover","condensed","responsive"),
full_width = TRUE
)42 Pemilihan Variabel
Variabel respon:
- number_of_casualties
Variabel prediktor:
- number_of_vehicles
- speed_limit
- urban_or_rural_area
- light_conditions
43 Persiapan Data
data_model <- data_raw %>%
select(
number_of_casualties,
number_of_vehicles,
speed_limit,
urban_or_rural_area,
light_conditions
) %>%
na.omit()
data_model <- data_model %>%
mutate(
urban_or_rural_area = factor(urban_or_rural_area),
light_conditions = factor(light_conditions)
)
glimpse(data_model)44 Eksplorasi Data
44.1 Statistik Deskriptif
data_model %>%
select(number_of_casualties,
number_of_vehicles,
speed_limit) %>%
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(
Mean = mean(value),
SD = sd(value),
Min = min(value),
Median = median(value),
Max = max(value)
) %>%
mutate(across(where(is.numeric), round, 2)) %>%
kable(caption = "Statistik deskriptif variabel numerik") %>%
kable_styling(
bootstrap_options = c("striped","hover","condensed","responsive")
)44.2 Distribusi Variabel Respon
ggplot(data_model,
aes(number_of_casualties)) +
geom_histogram(bins = 20,
fill = "#2f7f73",
color = "white") +
labs(
title = "Distribusi Jumlah Korban",
x = "Jumlah Korban",
y = "Frekuensi"
) +
theme_minimal()Interpretasi:
Sebagian besar kecelakaan menghasilkan jumlah korban yang relatif kecil, sedangkan kecelakaan dengan jumlah korban besar relatif jarang terjadi.
44.3 Distribusi Wilayah
data_model %>%
count(urban_or_rural_area) %>%
mutate(proporsi = n/sum(n)) %>%
kable(caption = "Distribusi wilayah") %>%
kable_styling(
bootstrap_options = c("striped","hover","condensed")
)45 Pemilihan Model
45.1 Landasan Pemilihan
Variabel respon berupa jumlah korban sehingga berbentuk data hitungan.
Misalkan:
\[ Y_i \sim Poisson(\mu_i) \]
Maka model regresi Poisson dituliskan sebagai:
\[ \log(\mu_i) = \beta_0+ \beta_1X_1+ \beta_2X_2+ \beta_3X_3+ \beta_4X_4 \]
dengan:
- X1 = number_of_vehicles
- X2 = speed_limit
- X3 = urban_or_rural_area
- X4 = light_conditions
46 Estimasi Model Regresi Poisson
fit_poisson <- glm(
number_of_casualties ~
number_of_vehicles +
speed_limit +
urban_or_rural_area +
light_conditions,
family = poisson(link = "log"),
data = data_model
)
summary(fit_poisson)47 Uji Signifikansi Simultan
anova(fit_poisson, test = "Chisq")Hipotesis:
\[ H_0:\beta_1=\beta_2=\cdots=\beta_k=0 \]
\[ H_1:\text{minimal ada satu }\beta_j \neq 0 \]
48 Ringkasan Koefisien
hasil <- broom::tidy(fit_poisson)
hasil %>%
mutate(across(where(is.numeric), round, 4)) %>%
kable(caption = "Koefisien model regresi Poisson") %>%
kable_styling(
bootstrap_options = c("striped","hover","condensed","responsive")
)49 Incidence Rate Ratio (IRR)
irr <- exp(coef(fit_poisson))
data.frame(
Variabel = names(irr),
IRR = irr
) %>%
mutate(IRR = round(IRR,4)) %>%
kable(caption = "Incidence Rate Ratio") %>%
kable_styling(
bootstrap_options = c("striped","hover","condensed")
)50 Goodness of Fit
50.1 Pseudo R-Square
pscl::pR2(fit_poisson)50.2 Deviance
data.frame(
Deviance = fit_poisson$deviance,
DF = fit_poisson$df.residual
)51 Model Regresi Poisson
coef(fit_poisson)Model akhir dibentuk menggunakan koefisien yang diperoleh dari output di atas.
52 Prediksi
data_model$prediksi <- predict(
fit_poisson,
type = "response"
)
head(data_model[,c("number_of_casualties","prediksi")])53 Visualisasi Aktual dan Prediksi
ggplot(data_model[1:1000,],
aes(number_of_casualties,prediksi)) +
geom_point(alpha=.4,color="#2f7f73") +
geom_abline(intercept=0,slope=1,
linetype=2,color="red") +
theme_minimal()54 Interpretasi Uji Simultan
Berdasarkan hasil analisis deviance menggunakan Likelihood Ratio Test, diperoleh nilai p-value yang lebih kecil dari 0,05. Oleh karena itu, hipotesis nol ditolak.
Hasil ini menunjukkan bahwa secara simultan terdapat paling sedikit satu variabel prediktor yang berpengaruh terhadap jumlah korban kecelakaan. Dengan demikian model regresi Poisson yang dibangun layak digunakan untuk menjelaskan variasi jumlah korban kecelakaan lalu lintas.
55 Interpretasi Uji Parsial
Berdasarkan tabel koefisien model, variabel number_of_vehicles dan speed_limit memiliki p-value yang lebih kecil dari 0,05 sehingga berpengaruh signifikan terhadap jumlah korban kecelakaan.
Variabel urban_or_rural_area juga menunjukkan pengaruh yang signifikan pada beberapa kategori wilayah. Sebaliknya, sebagian besar kategori pada variabel light_conditions memiliki p-value yang lebih besar dari 0,05 sehingga tidak menunjukkan pengaruh yang signifikan pada model.
56 Persamaan Model Akhir
Berdasarkan hasil estimasi parameter, model regresi Poisson yang diperoleh dapat dituliskan sebagai:
\[ \log(\hat{\mu}) = -0.1429 + 0.1246(\text{number\_of\_vehicles}) + 0.0037(\text{speed\_limit}) + 0.0672(\text{urban\_or\_rural\_area=2}) + 0.4456(\text{urban\_or\_rural\_area=3}) \]
dengan \(\hat{\mu}\) menyatakan rata-rata jumlah korban yang diprediksi oleh model.
57 Interpretasi Incidence Rate Ratio (IRR)
57.1 Jumlah Kendaraan yang Terlibat
Nilai koefisien untuk variabel jumlah kendaraan sebesar 0,1246 menghasilkan:
\[ IRR=e^{0.1246}=1.133 \]
Artinya, setiap penambahan satu kendaraan yang terlibat dalam kecelakaan akan meningkatkan rata-rata jumlah korban sebesar 1,133 kali atau sekitar 13,3%, dengan asumsi variabel lain tetap.
57.2 Batas Kecepatan
Koefisien batas kecepatan sebesar 0,0037 menghasilkan:
\[ IRR=e^{0.0037}=1.004 \]
Artinya, setiap kenaikan batas kecepatan sebesar 1 mph meningkatkan rata-rata jumlah korban sebesar 1,004 kali atau sekitar 0,4%.
57.3 Wilayah Kecelakaan
Untuk kategori wilayah tertentu diperoleh:
\[ IRR=e^{0.0672}=1.070 \]
Artinya rata-rata jumlah korban pada kategori wilayah tersebut sekitar 1,07 kali lebih tinggi dibandingkan kategori referensi.
57.4 Kondisi Pencahayaan
Kondisi pencahayaan tidak menunjukkan pengaruh yang signifikan setelah memperhitungkan jumlah kendaraan, batas kecepatan, dan karakteristik wilayah.
58 Pembahasan
Hasil analisis menunjukkan bahwa jumlah kendaraan yang terlibat merupakan faktor yang paling dominan memengaruhi jumlah korban kecelakaan. Semakin banyak kendaraan yang terlibat, semakin besar peluang terjadinya korban tambahan karena kompleksitas tabrakan yang meningkat.
Batas kecepatan juga berpengaruh positif terhadap jumlah korban. Jalan dengan batas kecepatan yang lebih tinggi cenderung menghasilkan kecelakaan dengan tingkat keparahan yang lebih besar sehingga meningkatkan jumlah korban yang terjadi.
Perbedaan karakteristik wilayah menunjukkan bahwa lingkungan tempat kecelakaan terjadi turut memengaruhi jumlah korban. Sementara itu, kondisi pencahayaan tidak menunjukkan pengaruh yang signifikan pada model yang dibangun.
59 Kesimpulan
Berdasarkan analisis regresi Poisson terhadap data kecelakaan lalu lintas Inggris tahun 2025 diperoleh bahwa model yang dibangun signifikan dalam menjelaskan jumlah korban kecelakaan.
Variabel jumlah kendaraan yang terlibat dalam kecelakaan, batas kecepatan jalan, dan kategori wilayah memiliki pengaruh signifikan terhadap jumlah korban kecelakaan.
Jumlah kendaraan merupakan faktor yang memiliki pengaruh paling kuat. Setiap tambahan satu kendaraan yang terlibat dalam kecelakaan meningkatkan rata-rata jumlah korban sekitar 13,3%.
Batas kecepatan juga berpengaruh positif terhadap jumlah korban, yang menunjukkan bahwa kecelakaan pada jalan dengan batas kecepatan lebih tinggi cenderung menghasilkan lebih banyak korban.
Sementara itu, kondisi pencahayaan tidak menunjukkan pengaruh yang signifikan terhadap jumlah korban pada model yang diperoleh.