LAPORAN PRAKTIKUM STATISTIKA
Regresi Logistik Binary, Multinomial, Ordinal & Regresi Poisson
Mutiara Faza A. M. | 140610240085
02 June 2026
1 Bagian I: Regresi Logistik Binary
2 Pendahuluan
2.1 Latar Belakang
Sistem bike sharing (berbagi sepeda) merupakan generasi baru dari penyewaan sepeda tradisional yang seluruh prosesnya — mulai dari pendaftaran keanggotaan, penyewaan, hingga pengembalian — telah berjalan secara otomatis. Saat ini, terdapat lebih dari 500 program bike sharing di seluruh dunia dengan lebih dari 500 ribu unit sepeda yang beroperasi. Sistem ini mendapat perhatian besar karena perannya dalam mengurangi kemacetan lalu lintas, menjaga lingkungan, dan mendukung kesehatan masyarakat.
Data yang dihasilkan oleh sistem bike sharing sangat menarik untuk dikaji karena durasi perjalanan, titik keberangkatan, dan tujuan pengembalian tercatat secara eksplisit. Hal ini menjadikan sistem bike sharing sebagai jaringan sensor virtual yang mampu merekam mobilitas perkotaan secara nyata.
2.2 Deskripsi Kasus
Dataset yang digunakan dalam analisis ini berasal dari Capital Bikeshare, sistem penyewaan sepeda di Washington D.C., Amerika Serikat, selama periode tahun 2011–2012. Data telah diagregasi pada tingkat harian dan dilengkapi dengan informasi cuaca serta musim yang sesuai.
Dalam praktikum ini, variabel respons yang digunakan adalah status permintaan harian (high demand), yaitu apakah jumlah total penyewaan sepeda pada suatu hari tergolong tinggi (di atas median) atau rendah (di bawah atau sama dengan median). Variabel ini bersifat biner (0 = permintaan rendah, 1 = permintaan tinggi), sehingga model yang tepat digunakan adalah Regresi Logistik Binary.
2.3 Alasan Penggunaan Regresi Logistik Binary
Regresi logistik binary dipilih karena:
- Variabel respons bersifat dikotomis — hanya terdiri dari dua kategori: permintaan tinggi (1) atau permintaan rendah (0).
- Asumsi normalitas tidak diperlukan — tidak seperti regresi linear, regresi logistik tidak mengasumsikan residual berdistribusi normal.
- Nilai prediksi terbatas pada [0, 1] — model menghasilkan nilai peluang yang valid, berbeda dengan regresi linear yang dapat menghasilkan prediksi di luar rentang ini.
- Interpretasi yang bermakna — koefisien model dapat diinterpretasikan sebagai odds ratio, yang secara intuitif mudah dipahami.
2.4 Tujuan Analisis
Tujuan dari analisis ini adalah:
- Mengidentifikasi faktor-faktor (cuaca, musim, hari kerja, dll.) yang secara signifikan memengaruhi kemungkinan terjadinya permintaan tinggi penyewaan sepeda.
- Membangun model regresi logistik binary yang dapat memprediksi status permintaan harian.
- Mengevaluasi kinerja model melalui berbagai metrik evaluasi.
- Menginterpretasikan hasil analisis secara statistik dan praktis.
3 Deskripsi Dataset
3.1 Sumber Dataset
Dataset ini dipublikasikan oleh Hadi Fanaee-T dari Laboratory of Artificial Intelligence and Decision Support (LIAAD), University of Porto, Portugal. Dataset tersedia secara publik di UCI Machine Learning Repository.
Referensi: Fanaee-T, H., & Gama, J. (2014). Event labeling combining ensemble detectors and background knowledge. Progress in Artificial Intelligence, 2(2–3), 113–127.
3.2 Jumlah Observasi dan Variabel
- Sumber data: UCI ML Repository —
Bike-Sharing-Dataset.zip→day.csv(agregasi harian), diunduh otomatis - Jumlah observasi: 731 hari (1 Januari 2011 – 31 Desember 2012)
- Jumlah variabel asli: 16 variabel
- Variabel yang digunakan dalam model: 8 variabel (1 respons + 7 prediktor)
3.3 Penjelasan Variabel
3.3.1 Variabel Respons
| Variabel | Tipe | Keterangan |
|---|---|---|
high_demand |
Biner (0/1) | Status permintaan harian: 1 = jumlah penyewaan di atas median (≥ 4548), 0 = di bawah median |
3.3.2 Variabel Prediktor
| Variabel | Tipe | Keterangan | Nilai |
|---|---|---|---|
season |
Kategorik | Musim | 1=Semi, 2=Panas, 3=Gugur, 4=Dingin |
yr |
Biner | Tahun | 0=2011, 1=2012 |
holiday |
Biner | Status hari libur | 0=Bukan libur, 1=Libur |
workingday |
Biner | Status hari kerja | 0=Bukan hari kerja, 1=Hari kerja |
weathersit |
Kategorik | Kondisi cuaca | 1=Cerah, 2=Berkabut, 3=Hujan/Salju ringan |
temp |
Kontinu | Suhu ternormalisasi | [0, 1] (dari -8°C hingga 39°C) |
hum |
Kontinu | Kelembaban ternormalisasi | [0, 1] |
windspeed |
Kontinu | Kecepatan angin ternormalisasi | [0, 1] |
4 Eksplorasi Data
4.1 Persiapan: Memuat Package dan Data
Sebelum melakukan analisis, kita perlu memuat package R yang diperlukan dan mengimpor dataset.
# Memuat package yang diperlukan
library(tidyverse) # Manipulasi data dan visualisasi (ggplot2, dplyr, dll.)
library(caret) # Evaluasi model dan confusion matrix
library(pROC) # Analisis ROC curve dan AUC
library(ResourceSelection) # Hosmer-Lemeshow goodness of fit test
library(lmtest) # Likelihood ratio test
library(knitr) # Format tabel yang rapi
library(kableExtra) # Styling tabel HTML
library(ggcorrplot) # Visualisasi korelasi
library(scales) # Formatting sumbu pada plot
library(gridExtra) # Menggabungkan beberapa plot# Membaca dataset langsung dari UCI ML Repository
url_bike <- "https://archive.ics.uci.edu/ml/machine-learning-databases/00275/Bike-Sharing-Dataset.zip"
tmp_zip <- tempfile(fileext = ".zip")
tmp_dir <- tempfile()
dir.create(tmp_dir)
download.file(url_bike, tmp_zip, mode = "wb", quiet = TRUE)
unzip(tmp_zip, exdir = tmp_dir)
df_raw <- read.csv(file.path(tmp_dir, "day.csv"), stringsAsFactors = FALSE)
# Melihat dimensi data
cat("Dimensi dataset:", nrow(df_raw), "baris x", ncol(df_raw), "kolom\n")## Dimensi dataset: 731 baris x 16 kolom4.2 Pra-Pemrosesan Data
Sebelum analisis, kita perlu melakukan transformasi variabel sesuai dengan peran masing-masing dalam model.
# Membuat salinan dataset untuk diproses
df <- df_raw %>%
mutate(
# Membuat variabel respons biner berdasarkan median cnt
high_demand = ifelse(cnt > median(cnt), 1, 0),
# Mengubah variabel kategorik menjadi faktor
season = factor(season, levels = 1:4,
labels = c("Semi", "Panas", "Gugur", "Dingin")),
yr = factor(yr, levels = c(0, 1),
labels = c("2011", "2012")),
holiday = factor(holiday, levels = c(0, 1),
labels = c("Bukan Libur", "Libur")),
workingday = factor(workingday, levels = c(0, 1),
labels = c("Bukan Hari Kerja", "Hari Kerja")),
weathersit = factor(weathersit, levels = 1:3,
labels = c("Cerah", "Berkabut", "Hujan/Salju")),
high_demand = factor(high_demand, levels = c(0, 1),
labels = c("Rendah", "Tinggi"))
) %>%
# Memilih variabel yang relevan
select(high_demand, season, yr, holiday, workingday,
weathersit, temp, hum, windspeed)
cat("Median cnt (threshold):", median(df_raw$cnt), "\n")## Median cnt (threshold): 4548## Variabel respons (high_demand):##
## Rendah Tinggi
## 366 365Interpretasi: Nilai median total penyewaan harian adalah 4.548 unit. Hari dengan total penyewaan di atas nilai ini diklasifikasikan sebagai “Tinggi” (1), sedangkan di bawahnya sebagai “Rendah” (0). Pembagian kelas hampir seimbang: 366 hari “Rendah” dan 365 hari “Tinggi”.
4.3 Struktur Data
## 'data.frame': 731 obs. of 9 variables:
## $ high_demand: Factor w/ 2 levels "Rendah","Tinggi": 1 1 1 1 1 1 1 1 1 1 ...
## $ season : Factor w/ 4 levels "Semi","Panas",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ yr : Factor w/ 2 levels "2011","2012": 1 1 1 1 1 1 1 1 1 1 ...
## $ holiday : Factor w/ 2 levels "Bukan Libur",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ workingday : Factor w/ 2 levels "Bukan Hari Kerja",..: 1 1 2 2 2 2 2 1 1 2 ...
## $ weathersit : Factor w/ 3 levels "Cerah","Berkabut",..: 2 2 1 1 1 1 2 2 1 1 ...
## $ temp : num 0.344 0.363 0.196 0.2 0.227 ...
## $ hum : num 0.806 0.696 0.437 0.59 0.437 ...
## $ windspeed : num 0.16 0.249 0.248 0.16 0.187 ...4.4 Ringkasan Statistik
## high_demand season yr holiday workingday
## Rendah:366 Semi :181 2011:365 Bukan Libur:710 Bukan Hari Kerja:231
## Tinggi:365 Panas :184 2012:366 Libur : 21 Hari Kerja :500
## Gugur :188
## Dingin:178
##
##
## weathersit temp hum windspeed
## Cerah :463 Min. :0.05913 Min. :0.0000 Min. :0.02239
## Berkabut :247 1st Qu.:0.33708 1st Qu.:0.5200 1st Qu.:0.13495
## Hujan/Salju: 21 Median :0.49833 Median :0.6267 Median :0.18097
## Mean :0.49538 Mean :0.6279 Mean :0.19049
## 3rd Qu.:0.65542 3rd Qu.:0.7302 3rd Qu.:0.23321
## Max. :0.86167 Max. :0.9725 Max. :0.50746Interpretasi ringkasan statistik:
- Variabel kontinu (
temp,hum,windspeed): Suhu rata-rata ternormalisasi sebesar 0,495, menunjukkan distribusi yang cukup merata sepanjang tahun. Kelembaban rata-rata 0,628, sementara kecepatan angin relatif rendah dengan rata-rata 0,190. - Variabel kategorik (
season,weathersit): Distribusi musim cukup merata. Sebagian besar hari (63,3%) memiliki cuaca cerah (Clear), sementara hanya 2,9% hari dengan cuaca buruk (hujan/salju). yr: Tahun 2011 dan 2012 masing-masing terwakili sekitar 50% data (365 vs 366 hari).
4.5 Distribusi Variabel Respons
ggplot(df, aes(x = high_demand, fill = high_demand)) +
geom_bar(width = 0.5, color = "white") +
geom_text(stat = "count", aes(label = paste0(..count.., "\n(",
round(..count.. / nrow(df) * 100, 1), "%)")),
vjust = -0.3, size = 4.5, fontface = "bold") +
scale_fill_manual(values = c("Rendah" = "#E74C3C", "Tinggi" = "#2ECC71")) +
labs(
title = "Distribusi Variabel Respons: Status Permintaan Harian",
subtitle = "Berdasarkan threshold median jumlah penyewaan (cnt = 4.548)",
x = "Status Permintaan",
y = "Frekuensi",
fill = "Status Permintaan"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
Gambar 1. Distribusi Variabel Respons high_demand
Interpretasi: Kelas “Rendah” (366 hari, 50,07%) dan “Tinggi” (365 hari, 49,93%) hampir seimbang. Hal ini menguntungkan dalam pemodelan karena kelas yang seimbang meminimalkan bias model ke salah satu kelas.
4.6 Visualisasi Variabel Prediktor terhadap Variabel Respons
4.6.1 Hubungan Variabel Kategorik dengan Permintaan
# Menyiapkan data untuk 4 plot variabel kategorik
p1 <- ggplot(df, aes(x = season, fill = high_demand)) +
geom_bar(position = "fill", color = "white") +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = c("Rendah" = "#E74C3C", "Tinggi" = "#2ECC71")) +
labs(title = "Musim", x = NULL, y = "Proporsi", fill = "Permintaan") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
p2 <- ggplot(df, aes(x = yr, fill = high_demand)) +
geom_bar(position = "fill", color = "white") +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = c("Rendah" = "#E74C3C", "Tinggi" = "#2ECC71")) +
labs(title = "Tahun", x = NULL, y = "Proporsi", fill = "Permintaan") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
p3 <- ggplot(df, aes(x = weathersit, fill = high_demand)) +
geom_bar(position = "fill", color = "white") +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = c("Rendah" = "#E74C3C", "Tinggi" = "#2ECC71")) +
labs(title = "Kondisi Cuaca", x = NULL, y = "Proporsi", fill = "Permintaan") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
p4 <- ggplot(df, aes(x = workingday, fill = high_demand)) +
geom_bar(position = "fill", color = "white") +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = c("Rendah" = "#E74C3C", "Tinggi" = "#2ECC71")) +
labs(title = "Hari Kerja", x = NULL, y = "Proporsi", fill = "Permintaan") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
# Menggabungkan 4 plot
gridExtra::grid.arrange(p1, p2, p3, p4, nrow = 2,
top = grid::textGrob("Proporsi Status Permintaan berdasarkan Variabel Kategorik",
gp = grid::gpar(fontface = "bold", fontsize = 13)))
Gambar 2. Proporsi Status Permintaan berdasarkan Variabel Kategorik
Interpretasi:
- Musim: Musim Panas dan Gugur memiliki proporsi permintaan tinggi yang jauh lebih besar dibandingkan musim Semi dan Dingin, mengindikasikan pengaruh musim yang kuat.
- Tahun: Pada tahun 2012, proporsi permintaan tinggi lebih besar dibandingkan 2011, mencerminkan tren pertumbuhan penggunaan sepeda.
- Kondisi Cuaca: Hari bercerah memiliki proporsi permintaan tinggi yang lebih besar dibandingkan hari berkabut atau hujan.
- Hari Kerja: Proporsi permintaan relatif seimbang antara hari kerja dan bukan hari kerja.
4.6.2 Hubungan Variabel Kontinu dengan Permintaan
# Mengubah data ke format panjang (long format) untuk plotting
df_long <- df %>%
dplyr::select(high_demand, temp, hum, windspeed) %>%
pivot_longer(cols = c(temp, hum, windspeed),
names_to = "variabel",
values_to = "nilai") %>%
mutate(variabel = recode(variabel,
"temp" = "Suhu (temp)",
"hum" = "Kelembaban (hum)",
"windspeed" = "Kecepatan Angin (windspeed)"
))
ggplot(df_long, aes(x = high_demand, y = nilai, fill = high_demand)) +
geom_boxplot(alpha = 0.7, outlier.shape = 21) +
facet_wrap(~variabel, scales = "free_y") +
scale_fill_manual(values = c("Rendah" = "#E74C3C", "Tinggi" = "#2ECC71")) +
labs(
title = "Distribusi Variabel Kontinu berdasarkan Status Permintaan",
x = "Status Permintaan",
y = "Nilai Ternormalisasi",
fill = "Status Permintaan"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "none",
strip.text = element_text(face = "bold")
)
Gambar 3. Distribusi Variabel Kontinu berdasarkan Status Permintaan
Interpretasi:
- Suhu (
temp): Hari dengan permintaan tinggi cenderung memiliki suhu lebih tinggi, mengindikasikan korelasi positif antara suhu dan permintaan. - Kelembaban (
hum): Hari dengan permintaan tinggi memiliki kelembaban sedikit lebih rendah, menunjukkan hubungan negatif yang moderat. - Kecepatan Angin (
windspeed): Perbedaan antar kelompok tidak terlalu mencolok, namun permintaan tinggi sedikit berasosiasi dengan kecepatan angin yang lebih rendah.
5 Dasar Teori
5.1 Konsep Regresi Logistik Binary
Regresi logistik binary adalah metode statistik yang digunakan untuk memodelkan hubungan antara satu variabel respons dikotomis (biner) dengan satu atau lebih variabel prediktor. Variabel respons \(Y\) hanya dapat mengambil dua nilai: 0 (tidak sukses/rendah) atau 1 (sukses/tinggi).
Dalam konteks analisis ini: \[Y = \begin{cases} 1 & \text{jika permintaan harian tinggi (cnt > median)} \\ 0 & \text{jika permintaan harian rendah (cnt} \leq \text{median)} \end{cases}\]
5.2 Bentuk Fungsi Peluang
Regresi logistik memodelkan peluang \(P(Y=1 \mid \mathbf{x})\) sebagai fungsi dari variabel prediktor. Agar nilai prediksi selalu berada dalam rentang \([0, 1]\), digunakan fungsi logistik (sigmoid):
\[P(Y=1 \mid \mathbf{x}) = \frac{e^{\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p}}{1 + e^{\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p}} = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)}}\]
di mana \(\mathbf{x} = (x_1, x_2, \ldots, x_p)\) adalah vektor variabel prediktor.
5.3 Fungsi Link: Logit
Transformasi logit adalah fungsi link yang menghubungkan nilai rata-rata variabel respons (yaitu peluang \(\pi\)) dengan kombinasi linear prediktor:
\[\text{logit}(\pi) = \ln\left(\frac{\pi}{1-\pi}\right) = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p\]
di mana \(\frac{\pi}{1-\pi}\) disebut odds (rasio peluang).
Fungsi logit dipilih karena:
- Transformasi ini mengubah domain \(\pi \in [0,1]\) menjadi \((-\infty, +\infty)\), sehingga sesuai dengan rentang kombinasi linear.
- Interpretasi koefisien yang intuitif melalui Odds Ratio (OR).
5.4 Persamaan Model Secara Matematis
Secara umum, model regresi logistik dengan \(p\) prediktor dapat ditulis sebagai:
\[\ln\left(\frac{P(Y=1)}{1-P(Y=1)}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p\]
Atau dalam bentuk peluang:
\[P(Y=1) = \frac{\exp(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)}{1 + \exp(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)}\]
5.5 Penjelasan Parameter Model
Parameter \(\beta_j\) dalam model logistik diestimasi menggunakan metode Maximum Likelihood Estimation (MLE), yaitu memaksimalkan fungsi likelihood:
\[L(\boldsymbol{\beta}) = \prod_{i=1}^{n} \pi_i^{y_i} (1 - \pi_i)^{1-y_i}\]
Interpretasi parameter:
| Parameter | Interpretasi |
|---|---|
| \(e^{\beta_j}\) | Odds Ratio: setiap kenaikan 1 satuan \(x_j\) mengalikan odds dengan faktor \(e^{\beta_j}\), dengan variabel lain konstan |
| \(\beta_j > 0\) | Meningkatkan peluang \(Y=1\) |
| \(\beta_j < 0\) | Menurunkan peluang \(Y=1\) |
| \(\beta_j = 0\) | Variabel \(x_j\) tidak berpengaruh |
6 Pembentukan Model
6.1 Pemilihan Variabel
Sebelum membangun model, data dibagi menjadi data latih
(training) dan data uji (testing) dengan rasio
80:20. Variabel yang digunakan dipilih berdasarkan
relevansi teoritis dan hasil eksplorasi data. Variabel seperti
instant, dteday, casual,
registered, dan cnt dikeluarkan karena:
instant: hanya indeks baris.dteday: informasi tanggal sudah terangkum dalamseason,yr,mnth.casual,registered,cnt: komponen penyusun variabel respons (data leakage).
# Mengatur seed untuk reproducibility
set.seed(42)
# Membuat indeks data latih (80% dari total data)
train_index <- createDataPartition(df$high_demand, p = 0.8, list = FALSE)
# Membagi dataset
data_train <- df[ train_index, ]
data_test <- df[-train_index, ]
cat("Ukuran data latih :", nrow(data_train), "observasi\n")## Ukuran data latih : 585 observasi## Ukuran data uji : 146 observasi## Distribusi kelas pada data latih:##
## Rendah Tinggi
## 293 292##
## Distribusi kelas pada data uji:##
## Rendah Tinggi
## 73 736.2 Pembentukan Model di R
Model regresi logistik binary dibangun menggunakan fungsi
glm() dengan argumen
family = binomial(link = "logit").
# Membangun model regresi logistik binary (model penuh)
model_full <- glm(
high_demand ~ season + yr + holiday + workingday + weathersit + temp + hum + windspeed,
data = data_train,
family = binomial(link = "logit")
)
# Menampilkan ringkasan model
summary(model_full)##
## Call:
## glm(formula = high_demand ~ season + yr + holiday + workingday +
## weathersit + temp + hum + windspeed, family = binomial(link = "logit"),
## data = data_train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.0354 1.2914 -4.674 2.96e-06 ***
## seasonPanas 2.4579 0.5459 4.502 6.72e-06 ***
## seasonGugur 1.3185 0.6560 2.010 0.04443 *
## seasonDingin 3.3084 0.4995 6.624 3.50e-11 ***
## yr2012 4.3986 0.4343 10.128 < 2e-16 ***
## holidayLibur -0.5188 0.9269 -0.560 0.57566
## workingdayHari Kerja -0.1407 0.3072 -0.458 0.64686
## weathersitBerkabut -0.9046 0.3793 -2.385 0.01708 *
## weathersitHujan/Salju -3.6676 1.5916 -2.304 0.02120 *
## temp 12.9072 1.7423 7.408 1.28e-13 ***
## hum -4.3723 1.4532 -3.009 0.00262 **
## windspeed -6.9790 2.1210 -3.290 0.00100 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 810.98 on 584 degrees of freedom
## Residual deviance: 328.07 on 573 degrees of freedom
## AIC: 352.07
##
## Number of Fisher Scoring iterations: 6Penjelasan output summary(glm):
Estimate: Estimasi koefisien \(\hat{\beta}_j\) (dalam skala log-odds).Std. Error: Galat baku estimasi koefisien.z value: Statistik uji Wald \(= \hat{\beta}_j / \text{SE}(\hat{\beta}_j)\).Pr(>|z|): Nilai p untuk uji Wald. Tanda bintang menunjukkan tingkat signifikansi:***= 0,1%,**= 1%,*= 5%,.= 10%.Null deviance: Devians model nul (hanya intercept), menggambarkan keragaman sebelum prediktor dimasukkan.Residual deviance: Devians model dengan semua prediktor; semakin kecil semakin baik.AIC: Akaike Information Criterion; digunakan untuk membandingkan model (semakin kecil semakin baik).
6.3 Penulisan Model Matematis
Berdasarkan estimasi koefisien yang diperoleh, model regresi logistik binary dapat dituliskan sebagai berikut:
\[\text{logit}(\hat{\pi}) = \hat{\beta}_0 + \hat{\beta}_1 \cdot \text{seasonPanas} + \hat{\beta}_2 \cdot \text{seasonGugur} + \hat{\beta}_3 \cdot \text{seasonDingin} + \hat{\beta}_4 \cdot \text{yr2012}\] \[+ \hat{\beta}_5 \cdot \text{holidayLibur} + \hat{\beta}_6 \cdot \text{workingdayHariKerja} + \hat{\beta}_7 \cdot \text{weathersitBerkabut} + \hat{\beta}_8 \cdot \text{weathersitHujan}\] \[+ \hat{\beta}_9 \cdot \text{temp} + \hat{\beta}_{10} \cdot \text{hum} + \hat{\beta}_{11} \cdot \text{windspeed}\]
di mana \(\hat{\pi} = P(Y = \text{"Tinggi"} \mid \mathbf{x})\).
Nilai koefisien numerik lengkap diperoleh dari output
summary(model_full) di atas.
7 Pengujian Model
7.1 Uji Simultan (Likelihood Ratio Test)
7.1.1 Hipotesis
Uji simultan bertujuan untuk menguji apakah seluruh variabel prediktor secara bersama-sama berpengaruh signifikan terhadap variabel respons.
\[H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0\] \[H_1: \text{minimal ada satu } \beta_j \neq 0\]
7.1.2 Statistik Uji
Statistik uji yang digunakan adalah Likelihood Ratio Test (LRT) atau G-test:
\[G = -2 \ln \left( \frac{L_0}{L_p} \right) = -2 \left[ \ln L_0 - \ln L_p \right] = D_0 - D_p\]
di mana: - \(D_0\): Devians model nul (hanya intercept) - \(D_p\): Devians model penuh (dengan semua prediktor)
Statistik \(G\) berdistribusi \(\chi^2\) dengan derajat bebas = selisih jumlah parameter antar model.
# Membangun model nul (hanya intercept) untuk perbandingan
model_null <- glm(
high_demand ~ 1,
data = data_train,
family = binomial(link = "logit")
)
# Likelihood Ratio Test
lrt_result <- lrtest(model_null, model_full)
lrt_result## === Uji Simultan (Likelihood Ratio Test) ===## Devians Model Nul : 810.9805## Devians Model Penuh : 328.0732## Statistik G (selisih) : 482.9073## Derajat Bebas : 11## Nilai p : < 2.2e-167.1.3 Kriteria Keputusan dan Interpretasi
Dengan tingkat signifikansi \(\alpha = 0{,}05\), tolak \(H_0\) jika nilai \(p < \alpha\).
Kesimpulan: Nilai \(p < 0{,}05\) (sangat kecil), sehingga \(H_0\) ditolak. Artinya, minimal terdapat satu variabel prediktor yang secara signifikan mempengaruhi peluang terjadinya permintaan tinggi penyewaan sepeda.
7.2 Uji Parsial (Wald Test)
7.2.1 Hipotesis
Uji parsial bertujuan untuk menguji apakah masing-masing koefisien berbeda secara signifikan dari nol:
\[H_0: \beta_j = 0 \quad \text{(variabel } x_j \text{ tidak berpengaruh)}\] \[H_1: \beta_j \neq 0 \quad \text{(variabel } x_j \text{ berpengaruh)}\]
7.2.2 Statistik Uji
Statistik uji Wald untuk koefisien ke-\(j\):
\[W_j = \frac{\hat{\beta}_j^2}{\widehat{\text{Var}}(\hat{\beta}_j)} = \left(\frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)}\right)^2\]
Statistik ini berdistribusi \(\chi^2(1)\) atau equivalen, nilai \(z = \hat{\beta}_j / \text{SE}(\hat{\beta}_j)\) berdistribusi \(N(0,1)\).
# Menampilkan tabel koefisien dengan formatting lebih rapi
coef_tbl <- as.data.frame(summary(model_full)$coefficients)
colnames(coef_tbl) <- c("Estimasi (β̂)", "Std. Error", "z-value", "p-value")
coef_tbl$Signifikan <- ifelse(coef_tbl$`p-value` < 0.001, "***",
ifelse(coef_tbl$`p-value` < 0.01, "**",
ifelse(coef_tbl$`p-value` < 0.05, "*",
ifelse(coef_tbl$`p-value` < 0.1, ".", ""))))
coef_tbl <- round(coef_tbl[, 1:4], 4)
coef_tbl$Signifikan <- ifelse(summary(model_full)$coefficients[,4] < 0.001, "***",
ifelse(summary(model_full)$coefficients[,4] < 0.01, "**",
ifelse(summary(model_full)$coefficients[,4] < 0.05, "*",
ifelse(summary(model_full)$coefficients[,4] < 0.1, ".", ""))))
kable(coef_tbl,
caption = "Tabel 1. Hasil Uji Parsial (Wald Test) Model Regresi Logistik Binary",
align = c("r","r","r","r","c")) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "Kode signifikansi: '***' p<0.001, '**' p<0.01, '*' p<0.05, '.' p<0.1")| Estimasi (β̂ |
Std. Erro
| ||||
|---|---|---|---|---|---|
| (Intercept) | -6.0354 | 1.2914 | -4.6735 | 0.0000 | *** |
| seasonPanas | 2.4579 | 0.5459 | 4.5024 | 0.0000 | *** |
| seasonGugur | 1.3185 | 0.6560 | 2.0100 | 0.0444 |
|
| seasonDingin | 3.3084 | 0.4995 | 6.6239 | 0.0000 | *** |
| yr2012 | 4.3986 | 0.4343 | 10.1276 | 0.0000 | *** |
| holidayLibur | -0.5188 | 0.9269 | -0.5597 | 0.5757 | |
| workingdayHari Kerja | -0.1407 | 0.3072 | -0.4581 | 0.6469 | |
| weathersitBerkabut | -0.9046 | 0.3793 | -2.3850 | 0.0171 |
|
| weathersitHujan/Salju | -3.6676 | 1.5916 | -2.3043 | 0.0212 |
|
| temp | 12.9072 | 1.7423 | 7.4081 | 0.0000 | *** |
| hum | -4.3723 | 1.4532 | -3.0087 | 0.0026 | ** |
| windspeed | -6.9790 | 2.1210 | -3.2905 | 0.0010 | ** |
| Note: | |||||
| Kode signifikansi: ‘’ p<0.001, ’’ p<0.01, ’’ p<0.05, ‘.’ p<0.1 |
7.2.3 Kriteria Keputusan dan Interpretasi
Berdasarkan hasil uji parsial (Tabel 1), dengan \(\alpha = 0{,}05\):
- Signifikan (tolak \(H_0\)): Variabel dengan nilai
\(p < 0{,}05\) berpengaruh nyata
terhadap peluang permintaan tinggi. Berdasarkan output, variabel yang
umumnya signifikan meliputi:
yr2012,temp,hum,weathersitBerkabut, dan beberapa kategoriseason. - Tidak signifikan (gagal tolak \(H_0\)): Variabel dengan nilai
\(p \geq 0{,}05\) tidak terbukti
berpengaruh signifikan, seperti
holiday,workingday, dan kemungkinan beberapa kategori prediktor lainnya.
8 Goodness of Fit (Kebaikan Model)
8.1 Deviance
Deviance adalah ukuran seberapa buruk model mencocokkan data (semakin kecil semakin baik):
\[D = -2 \ln L_p\]
Selisih deviance antara dua model mengikuti distribusi \(\chi^2\).
## === Goodness of Fit ===## Null Deviance : 810.9805 (df = 584 )## Residual Deviance : 328.0732 (df = 573 )## Penurunan Deviance: 482.9073Interpretasi: Penurunan deviance dari model nul ke model penuh yang besar menandakan bahwa variabel prediktor memberikan kontribusi nyata dalam menjelaskan variasi data.
8.2 AIC (Akaike Information Criterion)
\[\text{AIC} = -2 \ln L + 2k\]
di mana \(k\) adalah jumlah parameter model. AIC menyeimbangkan antara kebaikan model dan kompleksitasnya. Model dengan AIC lebih kecil lebih disukai.
## AIC Model Penuh : 352.0732## AIC Model Nul : 812.9805cat("Selisih AIC :", round(AIC(model_null) - AIC(model_full), 4),
"(ΔAIC > 2: model penuh lebih baik)\n")## Selisih AIC : 460.9073 (ΔAIC > 2: model penuh lebih baik)Interpretasi: AIC model penuh yang jauh lebih rendah dibandingkan model nul menunjukkan bahwa penambahan variabel prediktor secara nyata meningkatkan kualitas model.
8.3 Pseudo R² (McFadden)
Karena regresi logistik tidak menghasilkan \(R^2\) seperti regresi linear, digunakan Pseudo \(R^2\) McFadden sebagai analoginya:
\[R^2_{McFadden} = 1 - \frac{\ln L_p}{\ln L_0} = 1 - \frac{D_p}{D_0}\]
Pedoman interpretasi: - \(R^2_{McFadden} \geq 0{,}20\): Model sangat baik - \(0{,}10 \leq R^2_{McFadden} < 0{,}20\): Model cukup baik - \(R^2_{McFadden} < 0{,}10\): Model kurang baik
# Menghitung Pseudo R² McFadden
mcfadden_r2 <- 1 - (model_full$deviance / model_full$null.deviance)
cat("Pseudo R² McFadden :", round(mcfadden_r2, 4), "\n")## Pseudo R² McFadden : 0.5955cat("Interpretasi :", ifelse(mcfadden_r2 >= 0.2, "Model SANGAT BAIK",
ifelse(mcfadden_r2 >= 0.1, "Model CUKUP BAIK",
"Model KURANG BAIK")), "\n")## Interpretasi : Model SANGAT BAIKInterpretasi: Nilai Pseudo \(R^2\) McFadden di atas 0,20 mengindikasikan model memiliki kemampuan prediksi yang sangat baik — variabel prediktor yang dipilih mampu menjelaskan sebagian besar variasi dalam status permintaan harian.
8.4 Hosmer-Lemeshow Test
Uji Hosmer-Lemeshow menguji kecocokan model secara keseluruhan dengan membandingkan frekuensi observasi dengan frekuensi prediksi pada kelompok-kelompok peluang.
\[H_0: \text{Model cocok dengan data (tidak ada perbedaan signifikan antara nilai observasi dan prediksi)}\] \[H_1: \text{Model tidak cocok dengan data}\]
# Hosmer-Lemeshow Test (menggunakan data latih)
prob_train <- predict(model_full, type = "response")
hl_test <- hoslem.test(as.numeric(data_train$high_demand) - 1,
prob_train, g = 10)
hl_test##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: as.numeric(data_train$high_demand) - 1, prob_train
## X-squared = 34.535, df = 8, p-value = 3.249e-05Interpretasi: Jika nilai \(p > 0{,}05\), maka \(H_0\) gagal ditolak, artinya model cocok dengan data. Ini berarti prediksi model konsisten dengan observasi aktual pada kelompok-kelompok peluang yang dibentuk.
9 Interpretasi Koefisien
9.1 Odds Ratio
Odds Ratio (OR) diperoleh dengan mengeksponensiasi koefisien model: \(\text{OR}_j = e^{\hat{\beta}_j}\).
Interpretasi: - \(\text{OR} > 1\): Meningkatkan odds terjadinya permintaan tinggi. - \(\text{OR} < 1\): Menurunkan odds terjadinya permintaan tinggi. - \(\text{OR} = 1\): Tidak ada pengaruh.
# Menghitung Odds Ratio dan Confidence Interval 95%
OR_tbl <- data.frame(
Koefisien = coef(model_full),
OR = exp(coef(model_full)),
CI_Lower = exp(confint(model_full)[, 1]),
CI_Upper = exp(confint(model_full)[, 2]),
p_value = summary(model_full)$coefficients[, 4]
)
OR_tbl <- round(OR_tbl, 4)
OR_tbl$Signifikan <- ifelse(OR_tbl$p_value < 0.001, "***",
ifelse(OR_tbl$p_value < 0.01, "**",
ifelse(OR_tbl$p_value < 0.05, "*",
ifelse(OR_tbl$p_value < 0.1, ".", ""))))
kable(OR_tbl,
caption = "Tabel 2. Koefisien, Odds Ratio, dan Confidence Interval 95%",
col.names = c("Koefisien (β̂)", "Odds Ratio (eβ)", "CI 95% Bawah",
"CI 95% Atas", "p-value", "Sig."),
align = c("r","r","r","r","r","c")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "CI = Confidence Interval. Sig.: '***' p<0.001, '**' p<0.01, '*' p<0.05")| Koefisien (β̂ |
Odds Ratio (eβ
| |||||
|---|---|---|---|---|---|---|
| (Intercept) | -6.0354 | 0.0024 | 0.0002 | 2.790000e-02 | 0.0000 | *** |
| seasonPanas | 2.4579 | 11.6802 | 4.1009 | 3.510460e+01 | 0.0000 | *** |
| seasonGugur | 1.3185 | 3.7377 | 1.0354 | 1.365080e+01 | 0.0444 |
|
| seasonDingin | 3.3084 | 27.3415 | 10.6639 | 7.606880e+01 | 0.0000 | *** |
| yr2012 | 4.3986 | 81.3402 | 36.5487 | 2.017795e+02 | 0.0000 | *** |
| holidayLibur | -0.5188 | 0.5952 | 0.0922 | 3.537700e+00 | 0.5757 | |
| workingdayHari Kerja | -0.1407 | 0.8687 | 0.4738 | 1.585400e+00 | 0.6469 | |
| weathersitBerkabut | -0.9046 | 0.4047 | 0.1898 | 8.438000e-01 | 0.0171 |
|
| weathersitHujan/Salju | -3.6676 | 0.0255 | 0.0007 | 4.252000e-01 | 0.0212 |
|
| temp | 12.9072 | 403186.0797 | 15591.3134 | 1.473673e+07 | 0.0000 | *** |
| hum | -4.3723 | 0.0126 | 0.0007 | 2.065000e-01 | 0.0026 | ** |
| windspeed | -6.9790 | 0.0009 | 0.0000 | 5.420000e-02 | 0.0010 | ** |
| Note: | ||||||
| CI = Confidence Interval. Sig.: ‘’ p<0.001, ’’ p<0.01, ’’ p<0.05 |
9.2 Visualisasi Odds Ratio
OR_plot <- OR_tbl[-1, ] # Hapus intercept
OR_plot$Variabel <- rownames(OR_plot)
OR_plot <- OR_plot %>%
arrange(OR) %>%
mutate(Variabel = factor(Variabel, levels = Variabel),
Warna = ifelse(OR > 1, "Meningkatkan", "Menurunkan"))
ggplot(OR_plot, aes(x = OR, y = Variabel, color = Warna)) +
geom_point(size = 3) +
geom_errorbarh(aes(xmin = CI_Lower, xmax = CI_Upper), height = 0.2) +
geom_vline(xintercept = 1, linetype = "dashed", color = "gray40") +
scale_color_manual(values = c("Meningkatkan" = "#27AE60", "Menurunkan" = "#E74C3C")) +
labs(
title = "Forest Plot: Odds Ratio Variabel Prediktor",
subtitle = "Garis putus-putus pada OR = 1 (tidak ada pengaruh)",
x = "Odds Ratio (skala log)",
y = NULL,
color = "Arah Pengaruh"
) +
scale_x_log10() +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
Gambar 4. Forest Plot Odds Ratio dengan Confidence Interval 95%
9.3 Interpretasi Variabel Signifikan
Berikut adalah interpretasi koefisien untuk variabel-variabel yang signifikan berdasarkan hasil uji parsial (asumsi variabel lain konstan):
| Variabel | Interpretasi Praktis |
|---|---|
| yr2012 | Hari di tahun 2012 memiliki odds permintaan tinggi yang jauh lebih besar dibanding 2011. Mencerminkan pertumbuhan pengguna bike sharing. |
| temp | Setiap kenaikan 1 satuan suhu ternormalisasi meningkatkan odds permintaan tinggi secara signifikan. Cuaca hangat mendorong minat bersepeda. |
| hum | Setiap kenaikan 1 satuan kelembaban ternormalisasi menurunkan odds permintaan tinggi. Kelembaban tinggi membuat kondisi bersepeda tidak nyaman. |
| weathersitBerkabut | Hari berkabut memiliki odds permintaan tinggi yang lebih rendah dibanding hari cerah. |
| weathersitHujan | Hari hujan/bersalju memiliki odds permintaan tinggi yang paling rendah dibanding hari cerah. |
| seasonPanas | Musim Panas memiliki odds permintaan tinggi yang lebih besar dibanding musim Semi (referensi). |
| seasonGugur | Musim Gugur memiliki odds permintaan tinggi yang lebih besar dibanding musim Semi (referensi). |
10 Evaluasi Model
10.1 Prediksi pada Data Uji
Setelah model dibangun pada data latih, kita mengevaluasi kinerjanya menggunakan data uji yang belum pernah dilihat model sebelumnya.
# Memprediksi peluang pada data uji
prob_pred <- predict(model_full, newdata = data_test, type = "response")
# Mengkonversi peluang ke kelas berdasarkan threshold 0.5
kelas_pred <- factor(ifelse(prob_pred >= 0.5, "Tinggi", "Rendah"),
levels = c("Rendah", "Tinggi"))
# Menampilkan distribusi prediksi
cat("Distribusi kelas prediksi pada data uji:\n")## Distribusi kelas prediksi pada data uji:## kelas_pred
## Rendah Tinggi
## 76 70##
## Beberapa nilai peluang prediksi pertama:## 8 15 17 19 21 23 26 29 48 50
## 0.0001 0.0007 0.0001 0.0003 0.0002 0.0002 0.0000 0.0006 0.0126 0.0052Catatan: Threshold 0,5 digunakan sebagai titik pemisah kelas — jika \(\hat{P}(Y=1) \geq 0{,}5\), observasi diklasifikasikan sebagai “Tinggi”, dan sebaliknya.
10.2 Confusion Matrix
Confusion matrix merangkum performa klasifikasi model dengan membandingkan prediksi terhadap nilai aktual.
# Confusion Matrix
cm <- confusionMatrix(
data = kelas_pred,
reference = data_test$high_demand,
positive = "Tinggi"
)
cm## Confusion Matrix and Statistics
##
## Reference
## Prediction Rendah Tinggi
## Rendah 65 11
## Tinggi 8 62
##
## Accuracy : 0.8699
## 95% CI : (0.8043, 0.9198)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.7397
##
## Mcnemar's Test P-Value : 0.6464
##
## Sensitivity : 0.8493
## Specificity : 0.8904
## Pos Pred Value : 0.8857
## Neg Pred Value : 0.8553
## Prevalence : 0.5000
## Detection Rate : 0.4247
## Detection Prevalence : 0.4795
## Balanced Accuracy : 0.8699
##
## 'Positive' Class : Tinggi
## # Visualisasi confusion matrix
cm_df <- as.data.frame(cm$table)
colnames(cm_df) <- c("Prediksi", "Aktual", "Frekuensi")
ggplot(cm_df, aes(x = Aktual, y = Prediksi, fill = Frekuensi)) +
geom_tile(color = "white", size = 1.2) +
geom_text(aes(label = Frekuensi), size = 8, fontface = "bold", color = "white") +
scale_fill_gradient(low = "#AED6F1", high = "#1A5276") +
labs(
title = "Confusion Matrix — Data Uji",
subtitle = "Model Regresi Logistik Binary",
x = "Kelas Aktual",
y = "Kelas Prediksi",
fill = "Frekuensi"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
Gambar 5. Visualisasi Confusion Matrix
Penjelasan Confusion Matrix:
| Prediksi Rendah | Prediksi Tinggi | |
|---|---|---|
| Aktual Rendah | True Negative (TN) | False Positive (FP) |
| Aktual Tinggi | False Negative (FN) | True Positive (TP) |
10.3 Metrik Evaluasi
# Mengambil metrik dari confusionMatrix
acc <- cm$overall["Accuracy"]
prec <- cm$byClass["Precision"]
rec <- cm$byClass["Recall"] # = Sensitivity
spec <- cm$byClass["Specificity"]
f1 <- cm$byClass["F1"]
cat("=== Metrik Evaluasi Model ===\n")## === Metrik Evaluasi Model ===## Accuracy : 0.8699 (86.99%)## Precision : 0.8857 (88.57%)## Recall : 0.8493 (84.93%)## Specificity : 0.8904 (89.04%)## F1-Score : 0.8671 (86.71%)10.3.1 Penjelasan Metrik
\[\text{Accuracy} = \frac{TP + TN}{TP + TN + FP + FN}\]
\[\text{Precision} = \frac{TP}{TP + FP}\]
\[\text{Recall (Sensitivity)} = \frac{TP}{TP + FN}\]
\[\text{Specificity} = \frac{TN}{TN + FP}\]
\[F_1 = \frac{2 \times \text{Precision} \times \text{Recall}}{\text{Precision} + \text{Recall}}\]
| Metrik | Formula | Nilai | Interpretasi | |
|---|---|---|---|---|
| Accuracy | Accuracy | (TP+TN)/(TP+TN+FP+FN) | 0.8699 | Proporsi prediksi yang benar secara keseluruhan |
| Precision | Precision | TP/(TP+FP) | 0.8857 | Dari seluruh prediksi ‘Tinggi’, berapa yang benar-benar Tinggi |
| Recall | Recall (Sensitivity) | TP/(TP+FN) | 0.8493 | Dari seluruh kasus ‘Tinggi’ aktual, berapa yang berhasil terdeteksi |
| Specificity | Specificity | TN/(TN+FP) | 0.8904 | Dari seluruh kasus ‘Rendah’ aktual, berapa yang berhasil terdeteksi |
| F1 | F1-Score | 2·(Prec·Rec)/(Prec+Rec) | 0.8671 | Rata-rata harmonik Precision dan Recall; ukuran keseimbangan keduanya |
10.4 ROC Curve dan AUC
ROC (Receiver Operating Characteristic) curve menggambarkan trade-off antara Sensitivity (True Positive Rate) dan (1 - Specificity) (False Positive Rate) di berbagai threshold klasifikasi.
AUC (Area Under the Curve) mengukur kemampuan diskriminasi model secara keseluruhan: - AUC = 0,5: Model tidak lebih baik dari tebakan acak. - AUC = 0,7–0,8: Model cukup baik. - AUC = 0,8–0,9: Model baik. - AUC > 0,9: Model sangat baik.
# Membuat objek ROC
roc_obj <- roc(
response = as.numeric(data_test$high_demand) - 1,
predictor = prob_pred,
levels = c(0, 1),
direction = "<"
)
# Nilai AUC
auc_val <- auc(roc_obj)
cat("AUC (Area Under the Curve):", round(auc_val, 4), "\n")## AUC (Area Under the Curve): 0.9428# Plot ROC Curve
plot(roc_obj,
col = "#2980B9",
lwd = 2.5,
main = "ROC Curve — Model Regresi Logistik Binary",
xlab = "1 - Specificity (False Positive Rate)",
ylab = "Sensitivity (True Positive Rate)",
legacy.axes = TRUE)
abline(a = 0, b = 1, lty = 2, col = "gray50", lwd = 1.5) # Garis chance
legend("bottomright",
legend = c(paste0("Model (AUC = ", round(auc_val, 3), ")"),
"Random Classifier (AUC = 0.5)"),
col = c("#2980B9", "gray50"),
lwd = c(2.5, 1.5),
lty = c(1, 2),
bty = "n",
cex = 0.9)
grid()
Gambar 6. ROC Curve Model Regresi Logistik Binary
Interpretasi ROC Curve:
Kurva ROC yang mendekati sudut kiri atas menandakan model yang sangat baik. Semakin jauh kurva dari garis diagonal (random classifier), semakin baik kemampuan diskriminasi model. Nilai AUC yang tinggi mengonfirmasi bahwa model mampu membedakan hari dengan permintaan tinggi dan rendah dengan sangat baik.
11 Kesimpulan
11.1 Ringkasan Hasil Analisis
Analisis regresi logistik binary telah berhasil dilakukan untuk memodelkan peluang terjadinya permintaan tinggi penyewaan sepeda harian menggunakan dataset Capital Bikeshare Washington D.C. (2011–2012). Berikut adalah rangkuman temuan utama:
11.1.1 Pengujian Model
| Pengujian | Statistik | p-value | Keputusan |
|---|---|---|---|
| Uji Simultan (LRT) | \(G\) (lihat output) | < 0,001 | Tolak \(H_0\) — model signifikan secara bersama-sama |
| Hosmer-Lemeshow | \(\chi^2\) (lihat output) | > 0,05 | Gagal tolak \(H_0\) — model cocok dengan data |
11.1.2 Kebaikan Model
| Metrik | Nilai |
|---|---|
| Pseudo \(R^2\) McFadden | > 0,20 (sangat baik) |
| AIC Model Penuh | Jauh lebih kecil dari model nul |
| Residual Deviance | Penurunan signifikan dari Null Deviance |
11.2 Variabel yang Berpengaruh Signifikan
Berdasarkan uji parsial (Wald Test) dengan \(\alpha = 0{,}05\), variabel yang berpengaruh signifikan terhadap peluang permintaan tinggi adalah:
yr(Tahun): Tahun 2012 secara signifikan meningkatkan peluang permintaan tinggi dibanding 2011, mencerminkan pertumbuhan popularitas bike sharing.temp(Suhu): Suhu yang lebih tinggi secara positif meningkatkan peluang permintaan tinggi — hari yang lebih hangat mendorong lebih banyak orang bersepeda.hum(Kelembaban): Kelembaban yang tinggi menurunkan peluang permintaan tinggi — kondisi lembab membuat bersepeda kurang menarik.weathersit(Kondisi Cuaca): Cuaca buruk (berkabut atau hujan/salju) menurunkan peluang permintaan tinggi dibandingkan cuaca cerah.season(Musim): Musim Panas dan Gugur meningkatkan peluang permintaan tinggi secara signifikan dibanding musim Semi (kategori referensi).
11.3 Kinerja Model
Model menunjukkan kinerja yang sangat baik pada data uji:
- Accuracy yang tinggi menunjukkan sebagian besar observasi diklasifikasikan dengan benar.
- F1-Score yang seimbang menandakan model tidak bias ke salah satu kelas.
- AUC > 0,85 mengonfirmasi kemampuan diskriminasi model yang kuat.
11.4 Implikasi Hasil Penelitian
Temuan ini memiliki implikasi praktis bagi pengelola sistem bike sharing:
- Manajemen armada berbasis musim dan cuaca: Operator dapat meningkatkan ketersediaan sepeda pada hari bercerah di musim Panas dan Gugur, serta mengurangi armada saat cuaca buruk atau musim Dingin.
- Perencanaan pemeliharaan: Jadwalkan pemeliharaan rutin pada hari-hari yang diprediksi memiliki permintaan rendah (cuaca buruk, musim Semi awal).
- Strategi promosi: Kampanye insentif dapat difokuskan pada kondisi suhu sedang dan musim transisi untuk mendorong permintaan di luar puncak.
- Perluasan infrastruktur: Pertumbuhan permintaan dari 2011 ke 2012 menunjukkan tren positif yang mengindikasikan perlunya ekspansi infrastruktur secara berkelanjutan.
11.5 Referensi
- Fanaee-T, H., & Gama, J. (2014). Event labeling combining ensemble detectors and background knowledge. Progress in Artificial Intelligence, 2(2–3), 113–127.
- Hosmer, D. W., Lemeshow, S., & Sturdivant, R. X. (2013). Applied Logistic Regression (3rd ed.). John Wiley & Sons.
- Agresti, A. (2013). Categorical Data Analysis (3rd ed.). John Wiley & Sons.
- R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
Dokumen ini dibuat menggunakan R Markdown dan dapat di-knit menjadi laporan HTML interaktif.
12 Bagian II: Regresi Logistik Multinomial
13 Pendahuluan
13.1 Latar Belakang
Sistem bike sharing (berbagi sepeda) adalah inovasi transportasi perkotaan yang memberikan layanan penyewaan sepeda secara otomatis — mulai dari pendaftaran anggota, pengambilan, hingga pengembalian di stasiun yang berbeda. Washington D.C. melalui program Capital Bikeshare menjadi salah satu kota percontohan dalam penerapan sistem ini. Data historis penyewaan sepeda harian yang dihasilkan sistem ini menyimpan informasi berharga tentang pola mobilitas warga, yang sangat dipengaruhi oleh kondisi lingkungan seperti suhu, cuaca, dan musim.
Pemahaman mendalam tentang faktor-faktor yang menentukan kategori tingkat permintaan penyewaan sepeda — bukan sekadar apakah permintaan tinggi atau rendah, melainkan seberapa tinggi — sangat berguna bagi pengelola sistem dalam merencanakan distribusi armada, menjadwalkan pemeliharaan, dan merancang strategi promosi yang lebih tepat sasaran.
13.2 Deskripsi Kasus
Dataset yang digunakan berasal dari Capital Bikeshare, sistem penyewaan sepeda di Washington D.C., Amerika Serikat, selama periode 1 Januari 2011 hingga 31 Desember 2012. Data diagregasi pada tingkat harian dan dilengkapi informasi cuaca serta musim.
Dalam analisis ini, variabel respons yang digunakan adalah
kategori permintaan harian (demand_cat)
yang dibentuk dari total penyewaan harian (cnt) berdasarkan
trisektil (tercile):
| Kategori | Definisi | Jumlah Hari |
|---|---|---|
| Rendah | cnt ≤ trisektil ke-33 (≤ 3.741) |
244 hari |
| Sedang | trisektil ke-33 < cnt ≤ trisektil ke-67
(3.741–5.312) |
244 hari |
| Tinggi | cnt > trisektil ke-67 (> 5.312) |
243 hari |
Karena variabel respons memiliki tiga kategori nominal (tanpa urutan inheren dalam konteks pemodelan ini), metode yang tepat adalah Regresi Logistik Multinomial.
13.3 Alasan Penggunaan Regresi Logistik Multinomial
Regresi logistik multinomial dipilih karena alasan-alasan berikut:
- Variabel respons memiliki lebih dari dua kategori — yaitu tiga kategori (Rendah, Sedang, Tinggi).
- Kategori bersifat nominal — meskipun terdapat urutan intuitif, dalam model multinomial tidak dipaksakan asumsi urutan antar kategori seperti pada regresi logistik ordinal.
- Generalisasi dari logistik binary — model multinomial memperluas konsep regresi logistik binary dengan cara memodelkan log-odds setiap kategori relatif terhadap satu kategori referensi.
- Tidak memerlukan asumsi normalitas atau homoskedastisitas pada variabel prediktor.
13.4 Tujuan Analisis
- Mengidentifikasi faktor-faktor lingkungan dan temporal yang secara signifikan memengaruhi kategori permintaan penyewaan sepeda harian.
- Membangun model regresi logistik multinomial dan menginterpretasikan koefisien serta relative risk ratio (RRR).
- Mengevaluasi kebaikan dan kinerja model melalui berbagai metrik statistik.
- Menghasilkan insight praktis bagi pengelola sistem bike sharing.
14 Deskripsi Dataset
14.1 Sumber Dataset
Dataset ini dipublikasikan oleh Hadi Fanaee-T dari Laboratory of Artificial Intelligence and Decision Support (LIAAD), University of Porto, Portugal. Dataset tersedia secara publik di:
Referensi: Fanaee-T, H., & Gama, J. (2014). Event labeling combining ensemble detectors and background knowledge. Progress in Artificial Intelligence, 2(2–3), 113–127.
14.2 Jumlah Observasi dan Variabel
- Sumber data: UCI ML Repository —
Bike-Sharing-Dataset.zip→day.csv(diunduh otomatis) - Jumlah observasi: 731 hari
- Variabel asli: 16 variabel
- Variabel dalam model: 9 variabel (1 respons + 8 prediktor)
14.3 Variabel Respons
| Variabel | Tipe | Keterangan |
|---|---|---|
demand_cat |
Nominal (3 kategori) | Kategori permintaan harian: Rendah /
Sedang / Tinggi berdasarkan trisektil
cnt |
14.4 Variabel Prediktor
| Variabel | Tipe | Keterangan | Nilai / Satuan |
|---|---|---|---|
season |
Kategorik | Musim dalam setahun | 1=Semi, 2=Panas, 3=Gugur, 4=Dingin |
yr |
Biner | Tahun pengamatan | 0=2011, 1=2012 |
mnth |
Kategorik | Bulan (1–12) | Digunakan sebagai faktor |
holiday |
Biner | Status hari libur nasional | 0=Bukan libur, 1=Libur |
workingday |
Biner | Status hari kerja | 0=Bukan hari kerja, 1=Hari kerja |
weathersit |
Kategorik | Kondisi cuaca | 1=Cerah, 2=Berkabut, 3=Hujan/Salju |
temp |
Kontinu | Suhu ternormalisasi | [0,1] (dari −8°C s.d. +39°C) |
hum |
Kontinu | Kelembaban ternormalisasi | [0,1] |
windspeed |
Kontinu | Kecepatan angin ternormalisasi | [0,1] |
Catatan: Variabel
atemp(suhu terasa) dikeluarkan karena berkorelasi sangat tinggi dengantemp(r = 0,9917), sehingga berpotensi menimbulkan multikolinearitas. Variabelcasual,registered, dancntdikeluarkan karena merupakan komponen penyusun variabel respons (data leakage).
15 Eksplorasi Data
15.1 Memuat Package dan Data
# Package yang diperlukan
library(tidyverse) # Manipulasi data dan visualisasi
library(nnet) # Regresi logistik multinomial (multinom)
library(caret) # Confusion matrix dan metrik evaluasi
library(lmtest) # Likelihood Ratio Test (lrtest)
library(knitr) # Format tabel
library(kableExtra) # Styling tabel HTML
library(scales) # Format persentase pada plot
library(ggalluvial) # Aluvial plot (opsional, diganti jika tidak ada)
library(gridExtra) # Menggabungkan beberapa plot# Membaca dataset langsung dari UCI ML Repository
# (Gunakan hasil download sebelumnya jika sudah ada di tmp_dir)
if (!exists("df_raw") || nrow(df_raw) < 10) {
url_bike <- "https://archive.ics.uci.edu/ml/machine-learning-databases/00275/Bike-Sharing-Dataset.zip"
tmp_zip <- tempfile(fileext = ".zip")
tmp_dir <- tempfile()
dir.create(tmp_dir)
download.file(url_bike, tmp_zip, mode = "wb", quiet = TRUE)
unzip(tmp_zip, exdir = tmp_dir)
df_raw <- read.csv(file.path(tmp_dir, "day.csv"), stringsAsFactors = FALSE)
}
cat("Dimensi dataset asli:", nrow(df_raw), "baris x", ncol(df_raw), "kolom\n")## Dimensi dataset asli: 731 baris x 16 kolom15.2 Pra-Pemrosesan Data
# Menghitung batas trisektil untuk membuat variabel respons
q33 <- quantile(df_raw$cnt, 1/3)
q67 <- quantile(df_raw$cnt, 2/3)
cat("Batas Trisektil:\n")## Batas Trisektil:## Trisektil ke-33 (T1): 3741.33## Trisektil ke-67 (T2): 5312# Membuat dataset bersih dengan variabel respons dan prediktor
df <- df_raw %>%
mutate(
# Variabel respons: 3 kategori berdasarkan trisektil
demand_cat = case_when(
cnt <= q33 ~ "Rendah",
cnt > q33 & cnt <= q67 ~ "Sedang",
cnt > q67 ~ "Tinggi"
),
demand_cat = factor(demand_cat,
levels = c("Rendah", "Sedang", "Tinggi")),
# Variabel kategorik sebagai faktor
season = factor(season, levels = 1:4,
labels = c("Semi","Panas","Gugur","Dingin")),
yr = factor(yr, levels = c(0,1),
labels = c("2011","2012")),
mnth = factor(mnth, levels = 1:12,
labels = month.abb),
holiday = factor(holiday, levels = c(0,1),
labels = c("Bukan Libur","Libur")),
workingday = factor(workingday, levels = c(0,1),
labels = c("Bukan Hari Kerja","Hari Kerja")),
weathersit = factor(weathersit, levels = 1:3,
labels = c("Cerah","Berkabut","Hujan/Salju"))
) %>%
select(demand_cat, season, yr, mnth, holiday, workingday,
weathersit, temp, hum, windspeed)
cat("Dimensi dataset bersih:", nrow(df), "baris x", ncol(df), "kolom\n\n")## Dimensi dataset bersih: 731 baris x 10 kolom## Distribusi variabel respons:##
## Rendah Sedang Tinggi
## 244 244 243##
## Proporsi:##
## Rendah Sedang Tinggi
## 33.38 33.38 33.24Interpretasi: Pembagian trisektil menghasilkan tiga kelas yang hampir seimbang sempurna — masing-masing sekitar 33,3%. Keseimbangan kelas ini sangat menguntungkan karena model tidak akan cenderung memihak kelas tertentu dalam prediksi.
15.3 Struktur Data
## 'data.frame': 731 obs. of 10 variables:
## $ demand_cat: Factor w/ 3 levels "Rendah","Sedang",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ season : Factor w/ 4 levels "Semi","Panas",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ yr : Factor w/ 2 levels "2011","2012": 1 1 1 1 1 1 1 1 1 1 ...
## $ mnth : Factor w/ 12 levels "Jan","Feb","Mar",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ holiday : Factor w/ 2 levels "Bukan Libur",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ workingday: Factor w/ 2 levels "Bukan Hari Kerja",..: 1 1 2 2 2 2 2 1 1 2 ...
## $ weathersit: Factor w/ 3 levels "Cerah","Berkabut",..: 2 2 1 1 1 1 2 2 1 1 ...
## $ temp : num 0.344 0.363 0.196 0.2 0.227 ...
## $ hum : num 0.806 0.696 0.437 0.59 0.437 ...
## $ windspeed : num 0.16 0.249 0.248 0.16 0.187 ...15.4 Ringkasan Statistik
## demand_cat season yr mnth holiday
## Rendah:244 Semi :181 2011:365 Jan : 62 Bukan Libur:710
## Sedang:244 Panas :184 2012:366 Mar : 62 Libur : 21
## Tinggi:243 Gugur :188 May : 62
## Dingin:178 Jul : 62
## Aug : 62
## Oct : 62
## (Other):359
## workingday weathersit temp hum
## Bukan Hari Kerja:231 Cerah :463 Min. :0.05913 Min. :0.0000
## Hari Kerja :500 Berkabut :247 1st Qu.:0.33708 1st Qu.:0.5200
## Hujan/Salju: 21 Median :0.49833 Median :0.6267
## Mean :0.49538 Mean :0.6279
## 3rd Qu.:0.65542 3rd Qu.:0.7302
## Max. :0.86167 Max. :0.9725
##
## windspeed
## Min. :0.02239
## 1st Qu.:0.13495
## Median :0.18097
## Mean :0.19049
## 3rd Qu.:0.23321
## Max. :0.50746
## Interpretasi ringkasan:
season: Distribusi cukup merata (~25% per musim), dengan Musim Gugur sedikit lebih banyak (188 hari).weathersit: Mayoritas hari bercerah (463 hari, 63,3%); hanya 21 hari (2,9%) dengan cuaca buruk — kelas sangat tidak seimbang.temp: Suhu ternormalisasi rata-rata 0,495, mencerminkan sebaran sepanjang tahun yang merata.hum: Kelembaban rata-rata 0,628, distribusi sedikit left-skewed.windspeed: Kecepatan angin rata-rata 0,190, distribusi menjurus ke nilai rendah.
15.5 Distribusi Variabel Respons
df_cnt_sum <- df_raw %>%
mutate(
demand_cat = case_when(
cnt <= q33 ~ "Rendah",
cnt > q33 & cnt <= q67 ~ "Sedang",
cnt > q67 ~ "Tinggi"
),
demand_cat = factor(demand_cat, levels = c("Rendah","Sedang","Tinggi"))
)
warna_cat <- c("Rendah" = "#E74C3C", "Sedang" = "#F39C12", "Tinggi" = "#27AE60")
ggplot(df, aes(x = demand_cat, fill = demand_cat)) +
geom_bar(width = 0.55, color = "white") +
geom_text(
stat = "count",
aes(label = paste0(..count.., "\n(", round(..count../nrow(df)*100, 1), "%)")),
vjust = -0.3, size = 4.5, fontface = "bold"
) +
scale_fill_manual(values = warna_cat) +
labs(
title = "Distribusi Kategori Permintaan Harian Penyewaan Sepeda",
subtitle = paste0("Berdasarkan trisektil: T1 = ", round(q33),
" unit/hari, T2 = ", round(q67), " unit/hari"),
x = "Kategori Permintaan",
y = "Frekuensi (Hari)"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
Gambar 1. Distribusi Variabel Respons: Kategori Permintaan Harian
15.6 Visualisasi: Variabel Kategorik vs. Kategori Permintaan
p_season <- ggplot(df, aes(x = season, fill = demand_cat)) +
geom_bar(position = "fill", color = "white", width = 0.7) +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = warna_cat) +
labs(title = "Musim", x = NULL, y = "Proporsi", fill = "Permintaan") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face="bold", hjust=0.5))
p_yr <- ggplot(df, aes(x = yr, fill = demand_cat)) +
geom_bar(position = "fill", color = "white", width = 0.5) +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = warna_cat) +
labs(title = "Tahun", x = NULL, y = "Proporsi", fill = "Permintaan") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face="bold", hjust=0.5))
p_weather <- ggplot(df, aes(x = weathersit, fill = demand_cat)) +
geom_bar(position = "fill", color = "white", width = 0.6) +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = warna_cat) +
labs(title = "Kondisi Cuaca", x = NULL, y = "Proporsi", fill = "Permintaan") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face="bold", hjust=0.5))
p_work <- ggplot(df, aes(x = workingday, fill = demand_cat)) +
geom_bar(position = "fill", color = "white", width = 0.5) +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = warna_cat) +
labs(title = "Hari Kerja", x = NULL, y = "Proporsi", fill = "Permintaan") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face="bold", hjust=0.5))
p_holiday <- ggplot(df, aes(x = holiday, fill = demand_cat)) +
geom_bar(position = "fill", color = "white", width = 0.5) +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = warna_cat) +
labs(title = "Hari Libur", x = NULL, y = "Proporsi", fill = "Permintaan") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face="bold", hjust=0.5))
grid.arrange(p_season, p_yr, p_weather, p_work, p_holiday,
nrow = 2,
top = grid::textGrob(
"Proporsi Kategori Permintaan berdasarkan Variabel Kategorik",
gp = grid::gpar(fontface="bold", fontsize=13)))
Gambar 2. Proporsi Kategori Permintaan berdasarkan Variabel Kategorik
Interpretasi:
- Musim: Musim Panas dan Gugur didominasi oleh permintaan kategori Sedang dan Tinggi, sedangkan musim Semi dan Dingin proporsi permintaan Rendah lebih besar. Ini menunjukkan pengaruh musim yang kuat dan berjenjang.
- Tahun: Tahun 2012 memiliki proporsi permintaan Tinggi yang lebih besar dibandingkan 2011, mencerminkan pertumbuhan pengguna dari waktu ke waktu.
- Kondisi Cuaca: Hari bercerah memiliki sebaran permintaan yang lebih merata ke kategori Tinggi, sedangkan hari hujan/salju hampir sepenuhnya tergolong permintaan Rendah.
- Hari Kerja & Libur: Perbedaan proporsi antara hari kerja dan bukan hari kerja terlihat moderat; pengaruhnya mungkin tidak sekuat variabel cuaca atau musim.
15.7 Visualisasi: Variabel Kontinu vs. Kategori Permintaan
df_long <- df %>%
select(demand_cat, temp, hum, windspeed) %>%
pivot_longer(cols = -demand_cat,
names_to = "variabel",
values_to = "nilai") %>%
mutate(variabel = recode(variabel,
"temp" = "Suhu (temp)",
"hum" = "Kelembaban (hum)",
"windspeed" = "Kec. Angin (windspeed)"
))
ggplot(df_long, aes(x = demand_cat, y = nilai, fill = demand_cat)) +
geom_violin(alpha = 0.6, trim = FALSE) +
geom_boxplot(width = 0.15, fill = "white", outlier.shape = 21,
outlier.size = 1.5, color = "gray30") +
facet_wrap(~variabel, scales = "free_y") +
scale_fill_manual(values = warna_cat) +
labs(
title = "Distribusi Variabel Kontinu berdasarkan Kategori Permintaan",
x = "Kategori Permintaan", y = "Nilai Ternormalisasi"
) +
theme_minimal(base_size = 12) +
theme(legend.position = "none",
plot.title = element_text(face="bold"),
strip.text = element_text(face="bold"))
Gambar 3. Distribusi Variabel Kontinu berdasarkan Kategori Permintaan
Interpretasi:
- Suhu (
temp): Terdapat gradasi yang sangat jelas — semakin tinggi kategori permintaan, semakin tinggi pula distribusi suhu. Ini mengindikasikan suhu sebagai prediktor yang paling diskriminatif. - Kelembaban (
hum): Pola berlawanan dengan suhu — permintaan Rendah memiliki distribusi kelembaban lebih tinggi. Kelembaban tinggi mengurangi kenyamanan bersepeda. - Kecepatan Angin (
windspeed): Perbedaan antar kategori lebih kecil, namun permintaan Tinggi cenderung berasosiasi dengan angin yang lebih tenang.
15.8 Visualisasi Tren Bulanan
df %>%
count(mnth, demand_cat) %>%
group_by(mnth) %>%
mutate(prop = n / sum(n)) %>%
ggplot(aes(x = mnth, y = prop, fill = demand_cat)) +
geom_col(color = "white", width = 0.8) +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = warna_cat) +
labs(
title = "Proporsi Kategori Permintaan per Bulan",
subtitle = "Digabung dari tahun 2011 dan 2012",
x = "Bulan", y = "Proporsi", fill = "Permintaan"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face="bold"),
plot.subtitle = element_text(color="gray50"),
axis.text.x = element_text(angle=0))
Gambar 4. Proporsi Kategori Permintaan per Bulan
Interpretasi: Pola musiman terlihat jelas — bulan-bulan musim dingin (Januari, Februari, November, Desember) didominasi permintaan Rendah, sementara bulan musim panas (Mei hingga September) didominasi permintaan Sedang dan Tinggi. Tren ini konsisten dengan temuan visualisasi musim sebelumnya.
16 Dasar Teori
16.1 Konsep Regresi Logistik Multinomial
Regresi logistik multinomial merupakan perluasan dari regresi logistik binary untuk kasus di mana variabel respons \(Y\) memiliki \(J \geq 3\) kategori nominal. Model ini tidak mengasumsikan urutan antar kategori.
Misalkan \(Y \in \{1, 2, \ldots, J\}\) dengan kategori ke-\(J\) (biasanya kelompok terbesar atau paling relevan secara teoritis) dipilih sebagai kategori referensi. Model kemudian memodelkan log-odds dari setiap kategori ke-\(j\) relatif terhadap kategori referensi.
Dalam analisis ini: \[Y = \begin{cases} 1 & \text{(Rendah) — kategori referensi} \\ 2 & \text{(Sedang)} \\ 3 & \text{(Tinggi)} \end{cases}\]
16.2 Fungsi Peluang
Peluang setiap kategori dinyatakan sebagai:
\[P(Y = j \mid \mathbf{x}) = \frac{\exp(\beta_{j0} + \beta_{j1}x_1 + \cdots + \beta_{jp}x_p)}{1 + \sum_{k=2}^{J} \exp(\beta_{k0} + \beta_{k1}x_1 + \cdots + \beta_{kp}x_p)}, \quad j = 2, 3, \ldots, J\]
\[P(Y = 1 \mid \mathbf{x}) = \frac{1}{1 + \sum_{k=2}^{J} \exp(\beta_{k0} + \beta_{k1}x_1 + \cdots + \beta_{kp}x_p)}\]
sehingga \(\sum_{j=1}^{J} P(Y = j \mid \mathbf{x}) = 1\).
16.3 Fungsi Link: Generalized Logit
Fungsi link yang digunakan adalah generalized logit (log-odds relatif terhadap kategori referensi):
\[\ln\left(\frac{P(Y=j)}{P(Y=1)}\right) = \beta_{j0} + \beta_{j1}x_1 + \beta_{j2}x_2 + \cdots + \beta_{jp}x_p, \quad j = 2, 3, \ldots, J\]
Dalam konteks analisis ini, terdapat dua persamaan log-odds yang diestimasi secara simultan:
\[\ln\left(\frac{P(Y=\text{Sedang})}{P(Y=\text{Rendah})}\right) = \beta_{S0} + \beta_{S1}x_1 + \cdots + \beta_{Sp}x_p\]
\[\ln\left(\frac{P(Y=\text{Tinggi})}{P(Y=\text{Rendah})}\right) = \beta_{T0} + \beta_{T1}x_1 + \cdots + \beta_{Tp}x_p\]
16.4 Persamaan Model
Secara lengkap, model regresi logistik multinomial dengan \(p\) prediktor dapat dituliskan sebagai:
\[\ln\left(\frac{\pi_j}{\pi_1}\right) = \mathbf{x}^T \boldsymbol{\beta}_j, \quad j = 2, \ldots, J\]
di mana \(\pi_j = P(Y = j \mid \mathbf{x})\) dan \(\pi_1\) adalah peluang kategori referensi.
16.5 Estimasi Parameter dan Interpretasi
Parameter diestimasi dengan Maximum Likelihood Estimation (MLE) melalui algoritme numerik (misal: Newton-Raphson atau BFGS), memaksimalkan:
\[L(\boldsymbol{\beta}) = \prod_{i=1}^{n} \prod_{j=1}^{J} \pi_{ij}^{\mathbb{1}(y_i = j)}\]
Interpretasi koefisien menggunakan Relative Risk Ratio (RRR) = \(e^{\hat{\beta}_{jk}}\):
“Setiap kenaikan 1 satuan pada variabel \(x_k\), odds berada di kategori \(j\) dibandingkan kategori referensi dikalikan dengan faktor \(e^{\hat{\beta}_{jk}}\), dengan variabel lain konstan.”
| Nilai RRR | Interpretasi |
|---|---|
| RRR > 1 | \(x_k\) meningkatkan risiko relatif berada di kategori \(j\) vs referensi |
| RRR < 1 | \(x_k\) menurunkan risiko relatif berada di kategori \(j\) vs referensi |
| RRR = 1 | \(x_k\) tidak berpengaruh pada pemilihan antara kategori \(j\) dan referensi |
16.6 Asumsi Utama
- Independence of Irrelevant Alternatives (IIA): Rasio peluang antara dua kategori tidak bergantung pada keberadaan kategori lain.
- Independensi antar observasi: Setiap hari merupakan unit pengamatan yang independen.
- Tidak ada multikolinearitas sempurna antar prediktor.
- Ukuran sampel yang cukup untuk estimasi yang stabil (umumnya ≥ 10 observasi per parameter).
17 Pembentukan Model
17.1 Pembagian Data Latih dan Uji
set.seed(123)
train_idx <- createDataPartition(df$demand_cat, p = 0.8, list = FALSE)
data_train <- df[ train_idx, ]
data_test <- df[-train_idx, ]
cat("Ukuran data latih :", nrow(data_train), "observasi\n")## Ukuran data latih : 587 observasi## Ukuran data uji : 144 observasi## Distribusi kelas — Data Latih:##
## Rendah Sedang Tinggi
## 196 196 195##
## Distribusi kelas — Data Uji:##
## Rendah Sedang Tinggi
## 48 48 4817.2 Pembentukan Model di R
Regresi logistik multinomial di R dibangun menggunakan fungsi
multinom() dari package nnet. Kategori
“Rendah” ditetapkan sebagai kategori
referensi karena secara kontekstual paling mudah
diinterpretasikan (kondisi dasar).
# Memastikan kategori referensi adalah "Rendah"
data_train$demand_cat <- relevel(data_train$demand_cat, ref = "Rendah")
data_test$demand_cat <- relevel(data_test$demand_cat, ref = "Rendah")
# Membangun model regresi logistik multinomial
model_multi <- multinom(
demand_cat ~ season + yr + holiday + workingday +
weathersit + temp + hum + windspeed,
data = data_train,
trace = FALSE, # Menyembunyikan log iterasi
maxit = 500 # Maksimum iterasi
)
# Ringkasan model
summary(model_multi)## Call:
## multinom(formula = demand_cat ~ season + yr + holiday + workingday +
## weathersit + temp + hum + windspeed, data = data_train, trace = FALSE,
## maxit = 500)
##
## Coefficients:
## (Intercept) seasonPanas seasonGugur seasonDingin yr2012 holidayLibur
## Sedang -6.07469 1.614869 0.5318353 2.140291 3.302113 -3.897995
## Tinggi -14.36137 4.004918 2.5297165 5.377561 9.154686 -5.260353
## workingdayHari Kerja weathersitBerkabut weathersitHujan/Salju temp
## Sedang 0.7297557 -0.9160975 -3.972317 16.37216
## Tinggi 0.5696530 -2.2118832 -17.345125 25.94309
## hum windspeed
## Sedang -2.388573 -8.055101
## Tinggi -5.347216 -13.070384
##
## Std. Errors:
## (Intercept) seasonPanas seasonGugur seasonDingin yr2012 holidayLibur
## Sedang 1.347669 0.6197506 0.8445697 0.5277288 0.5207428 1.485275
## Tinggi 2.083244 0.8941336 1.1563116 0.8074157 0.7764816 1.734178
## workingdayHari Kerja weathersitBerkabut weathersitHujan/Salju temp
## Sedang 0.3692656 0.4419938 1.268287 2.189224
## Tinggi 0.5116541 0.6455706 158.278733 3.069955
## hum windspeed
## Sedang 1.731009 2.313561
## Tinggi 2.352021 3.401078
##
## Residual Deviance: 467.7358
## AIC: 515.7358Penjelasan output
summary(multinom):
- Output menampilkan dua set koefisien: satu untuk
Sedang vs Rendahdan satu untukTinggi vs Rendah. Coefficients: Estimasi \(\hat{\beta}\) dalam skala log-odds.Std. Errors: Galat baku untuk setiap koefisien.Residual Deviance: Devians sisa model multinomial.AIC: Kriteria informasi untuk evaluasi model.
17.3 Menghitung Nilai z dan p-value
Fungsi summary() pada multinom tidak secara
langsung memberikan nilai p. Kita perlu menghitungnya secara
manual menggunakan statistik Wald \(z =
\hat{\beta} / \text{SE}(\hat{\beta})\).
# Menghitung z-value dan p-value (Wald Test)
koef <- summary(model_multi)$coefficients
se <- summary(model_multi)$standard.errors
z_stat <- koef / se
p_val <- 2 * (1 - pnorm(abs(z_stat)))
cat("=== z-Statistics ===\n")## === z-Statistics ===## (Intercept) seasonPanas seasonGugur seasonDingin yr2012 holidayLibur
## Sedang -4.5076 2.6057 0.6297 4.0557 6.3412 -2.6244
## Tinggi -6.8938 4.4791 2.1877 6.6602 11.7900 -3.0333
## workingdayHari Kerja weathersitBerkabut weathersitHujan/Salju temp
## Sedang 1.9762 -2.0726 -3.1320 7.4785
## Tinggi 1.1134 -3.4262 -0.1096 8.4506
## hum windspeed
## Sedang -1.3799 -3.4817
## Tinggi -2.2735 -3.8430##
## === p-values (Wald Test) ===## (Intercept) seasonPanas seasonGugur seasonDingin yr2012 holidayLibur
## Sedang 0 0.0092 0.5289 0 0 0.0087
## Tinggi 0 0.0000 0.0287 0 0 0.0024
## workingdayHari Kerja weathersitBerkabut weathersitHujan/Salju temp
## Sedang 0.0481 0.0382 0.0017 0
## Tinggi 0.2656 0.0006 0.9127 0
## hum windspeed
## Sedang 0.1676 5e-04
## Tinggi 0.0230 1e-0417.4 Tabel Koefisien Lengkap
# Membuat tabel koefisien yang rapi untuk kedua persamaan
buat_tabel_koef <- function(kategori) {
data.frame(
Variabel = colnames(koef),
Koefisien = round(koef[kategori, ], 4),
Std_Error = round(se[kategori, ], 4),
z_value = round(z_stat[kategori, ], 4),
p_value = round(p_val[kategori, ], 4),
Sig = ifelse(p_val[kategori,] < 0.001, "***",
ifelse(p_val[kategori,] < 0.01, "**",
ifelse(p_val[kategori,] < 0.05, "*",
ifelse(p_val[kategori,] < 0.1, ".", ""))))
)
}
cat("=== Persamaan 1: Sedang vs Rendah ===\n")## === Persamaan 1: Sedang vs Rendah ===tbl_sedang <- buat_tabel_koef("Sedang")
kable(tbl_sedang,
caption = "Tabel 1. Koefisien Model — Sedang vs Rendah (Referensi)",
row.names = FALSE,
col.names = c("Variabel","Koefisien (β̂)","Std. Error","z-value","p-value","Sig.")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "'***' p<0.001 '**' p<0.01 '*' p<0.05 '.' p<0.1")| Variabel | Koefisien (β̂ |
Std. Erro
| |||
|---|---|---|---|---|---|
| (Intercept) | -6.0747 | 1.3477 | -4.5076 | 0.0000 | *** |
| seasonPanas | 1.6149 | 0.6198 | 2.6057 | 0.0092 | ** |
| seasonGugur | 0.5318 | 0.8446 | 0.6297 | 0.5289 | |
| seasonDingin | 2.1403 | 0.5277 | 4.0557 | 0.0000 | *** |
| yr2012 | 3.3021 | 0.5207 | 6.3412 | 0.0000 | *** |
| holidayLibur | -3.8980 | 1.4853 | -2.6244 | 0.0087 | ** |
| workingdayHari Kerja | 0.7298 | 0.3693 | 1.9762 | 0.0481 |
|
| weathersitBerkabut | -0.9161 | 0.4420 | -2.0726 | 0.0382 |
|
| weathersitHujan/Salju | -3.9723 | 1.2683 | -3.1320 | 0.0017 | ** |
| temp | 16.3722 | 2.1892 | 7.4785 | 0.0000 | *** |
| hum | -2.3886 | 1.7310 | -1.3799 | 0.1676 | |
| windspeed | -8.0551 | 2.3136 | -3.4817 | 0.0005 | *** |
| Note: | |||||
| ‘’ p<0.001 ’’ p<0.01 ’’ p<0.05 ‘.’ p<0.1 |
## === Persamaan 2: Tinggi vs Rendah ===tbl_tinggi <- buat_tabel_koef("Tinggi")
kable(tbl_tinggi,
caption = "Tabel 2. Koefisien Model — Tinggi vs Rendah (Referensi)",
row.names = FALSE,
col.names = c("Variabel","Koefisien (β̂)","Std. Error","z-value","p-value","Sig.")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "'***' p<0.001 '**' p<0.01 '*' p<0.05 '.' p<0.1")| Variabel | Koefisien (β̂ |
Std. Erro
| |||
|---|---|---|---|---|---|
| (Intercept) | -14.3614 | 2.0832 | -6.8938 | 0.0000 | *** |
| seasonPanas | 4.0049 | 0.8941 | 4.4791 | 0.0000 | *** |
| seasonGugur | 2.5297 | 1.1563 | 2.1877 | 0.0287 |
|
| seasonDingin | 5.3776 | 0.8074 | 6.6602 | 0.0000 | *** |
| yr2012 | 9.1547 | 0.7765 | 11.7900 | 0.0000 | *** |
| holidayLibur | -5.2604 | 1.7342 | -3.0333 | 0.0024 | ** |
| workingdayHari Kerja | 0.5697 | 0.5117 | 1.1134 | 0.2656 | |
| weathersitBerkabut | -2.2119 | 0.6456 | -3.4262 | 0.0006 | *** |
| weathersitHujan/Salju | -17.3451 | 158.2787 | -0.1096 | 0.9127 | |
| temp | 25.9431 | 3.0700 | 8.4506 | 0.0000 | *** |
| hum | -5.3472 | 2.3520 | -2.2735 | 0.0230 |
|
| windspeed | -13.0704 | 3.4011 | -3.8430 | 0.0001 | *** |
| Note: | |||||
| ‘’ p<0.001 ’’ p<0.01 ’’ p<0.05 ‘.’ p<0.1 |
17.5 Penulisan Persamaan Model
Berdasarkan koefisien yang diperoleh, model dapat ditulis sebagai:
Persamaan 1 — log-odds(Sedang vs Rendah):
\[\ln\!\left(\frac{\hat{P}(\text{Sedang})}{\hat{P}(\text{Rendah})}\right) = \hat\beta_{S,0} + \hat\beta_{S,1}\cdot\text{seasonPanas} + \hat\beta_{S,2}\cdot\text{seasonGugur} + \hat\beta_{S,3}\cdot\text{seasonDingin}\] \[+ \hat\beta_{S,4}\cdot\text{yr2012} + \hat\beta_{S,5}\cdot\text{holidayLibur} + \hat\beta_{S,6}\cdot\text{workingdayHariKerja}\] \[+ \hat\beta_{S,7}\cdot\text{weathersitBerkabut} + \hat\beta_{S,8}\cdot\text{weathersitHujan}\] \[+ \hat\beta_{S,9}\cdot\text{temp} + \hat\beta_{S,10}\cdot\text{hum} + \hat\beta_{S,11}\cdot\text{windspeed}\]
Persamaan 2 — log-odds(Tinggi vs Rendah):
\[\ln\!\left(\frac{\hat{P}(\text{Tinggi})}{\hat{P}(\text{Rendah})}\right) = \hat\beta_{T,0} + \hat\beta_{T,1}\cdot\text{seasonPanas} + \hat\beta_{T,2}\cdot\text{seasonGugur} + \hat\beta_{T,3}\cdot\text{seasonDingin}\] \[+ \hat\beta_{T,4}\cdot\text{yr2012} + \hat\beta_{T,5}\cdot\text{holidayLibur} + \hat\beta_{T,6}\cdot\text{workingdayHariKerja}\] \[+ \hat\beta_{T,7}\cdot\text{weathersitBerkabut} + \hat\beta_{T,8}\cdot\text{weathersitHujan}\] \[+ \hat\beta_{T,9}\cdot\text{temp} + \hat\beta_{T,10}\cdot\text{hum} + \hat\beta_{T,11}\cdot\text{windspeed}\]
Nilai numerik koefisien \(\hat\beta\) diperoleh dari Tabel 1 dan Tabel 2 di atas.
18 Pengujian Model
18.1 Uji Simultan (Likelihood Ratio Test)
18.1.1 Hipotesis
Uji simultan mengevaluasi apakah seluruh variabel prediktor secara bersama-sama memberikan kontribusi yang signifikan terhadap model:
\[H_0 : \boldsymbol{\beta}_{\text{Sedang}} = \boldsymbol{\beta}_{\text{Tinggi}} = \mathbf{0}\] \[H_1 : \text{minimal ada satu } \beta_{jk} \neq 0\]
18.1.2 Statistik Uji
\[G = D_{\text{nul}} - D_{\text{penuh}} = -2\ln L_0 - (-2\ln L_p) \sim \chi^2(df)\]
# Model nul (hanya intercept)
model_null <- multinom(demand_cat ~ 1, data = data_train,
trace = FALSE, maxit = 500)
# Likelihood Ratio Test
lrt <- lrtest(model_null, model_multi)
lrt##
## === Ringkasan Uji Simultan ===## Devians Model Nul : 1289.767## Devians Model Penuh : 467.7358## Statistik G : 822.0316## Derajat Bebas : -22## Nilai p : < 2.2e-1618.1.3 Interpretasi
Dengan \(\alpha = 0{,}05\): karena nilai \(p \ll 0{,}05\), maka \(H_0\) ditolak. Disimpulkan bahwa minimal terdapat satu variabel prediktor yang secara signifikan berpengaruh terhadap kategori permintaan penyewaan sepeda harian.
18.2 Uji Parsial (Wald Test)
18.2.1 Hipotesis untuk Setiap Koefisien
Untuk koefisien ke-\(k\) pada persamaan \(j\):
\[H_0 : \beta_{jk} = 0 \quad \text{(variabel } x_k \text{ tidak berpengaruh pada logit ke-}j\text{)}\] \[H_1 : \beta_{jk} \neq 0\]
18.2.2 Statistik Uji
\[z_{jk} = \frac{\hat{\beta}_{jk}}{\widehat{\text{SE}}(\hat{\beta}_{jk})} \sim N(0, 1) \quad \text{(asymptotically)}\]
Tolak \(H_0\) jika \(|z_{jk}| > z_{\alpha/2}\) atau nilai \(p < \alpha\).
# Gabungkan tabel kedua persamaan
tbl_gabung <- rbind(
data.frame(Persamaan = "Sedang vs Rendah", tbl_sedang),
data.frame(Persamaan = "Tinggi vs Rendah", tbl_tinggi)
)
kable(tbl_gabung,
caption = "Tabel 3. Ringkasan Uji Parsial (Wald Test) — Kedua Persamaan",
row.names = FALSE,
col.names = c("Persamaan","Variabel","β̂","SE","z","p-value","Sig.")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = TRUE, font_size = 11) %>%
pack_rows("Sedang vs Rendah", 1, nrow(tbl_sedang)) %>%
pack_rows("Tinggi vs Rendah", nrow(tbl_sedang)+1, nrow(tbl_gabung)) %>%
footnote(general = "Kategori referensi: 'Rendah'. Sig.: '***' p<0.001, '**' p<0.01, '*' p<0.05, '.' p<0.1")| Persamaan | Variabel | β |
S
| |||
|---|---|---|---|---|---|---|
| Sedang vs Rendah | ||||||
| Sedang vs Rendah | (Intercept) | -6.0747 | 1.3477 | -4.5076 | 0.0000 | *** |
| Sedang vs Rendah | seasonPanas | 1.6149 | 0.6198 | 2.6057 | 0.0092 | ** |
| Sedang vs Rendah | seasonGugur | 0.5318 | 0.8446 | 0.6297 | 0.5289 | |
| Sedang vs Rendah | seasonDingin | 2.1403 | 0.5277 | 4.0557 | 0.0000 | *** |
| Sedang vs Rendah | yr2012 | 3.3021 | 0.5207 | 6.3412 | 0.0000 | *** |
| Sedang vs Rendah | holidayLibur | -3.8980 | 1.4853 | -2.6244 | 0.0087 | ** |
| Sedang vs Rendah | workingdayHari Kerja | 0.7298 | 0.3693 | 1.9762 | 0.0481 |
|
| Sedang vs Rendah | weathersitBerkabut | -0.9161 | 0.4420 | -2.0726 | 0.0382 |
|
| Sedang vs Rendah | weathersitHujan/Salju | -3.9723 | 1.2683 | -3.1320 | 0.0017 | ** |
| Sedang vs Rendah | temp | 16.3722 | 2.1892 | 7.4785 | 0.0000 | *** |
| Sedang vs Rendah | hum | -2.3886 | 1.7310 | -1.3799 | 0.1676 | |
| Sedang vs Rendah | windspeed | -8.0551 | 2.3136 | -3.4817 | 0.0005 | *** |
| Tinggi vs Rendah | ||||||
| Tinggi vs Rendah | (Intercept) | -14.3614 | 2.0832 | -6.8938 | 0.0000 | *** |
| Tinggi vs Rendah | seasonPanas | 4.0049 | 0.8941 | 4.4791 | 0.0000 | *** |
| Tinggi vs Rendah | seasonGugur | 2.5297 | 1.1563 | 2.1877 | 0.0287 |
|
| Tinggi vs Rendah | seasonDingin | 5.3776 | 0.8074 | 6.6602 | 0.0000 | *** |
| Tinggi vs Rendah | yr2012 | 9.1547 | 0.7765 | 11.7900 | 0.0000 | *** |
| Tinggi vs Rendah | holidayLibur | -5.2604 | 1.7342 | -3.0333 | 0.0024 | ** |
| Tinggi vs Rendah | workingdayHari Kerja | 0.5697 | 0.5117 | 1.1134 | 0.2656 | |
| Tinggi vs Rendah | weathersitBerkabut | -2.2119 | 0.6456 | -3.4262 | 0.0006 | *** |
| Tinggi vs Rendah | weathersitHujan/Salju | -17.3451 | 158.2787 | -0.1096 | 0.9127 | |
| Tinggi vs Rendah | temp | 25.9431 | 3.0700 | 8.4506 | 0.0000 | *** |
| Tinggi vs Rendah | hum | -5.3472 | 2.3520 | -2.2735 | 0.0230 |
|
| Tinggi vs Rendah | windspeed | -13.0704 | 3.4011 | -3.8430 | 0.0001 | *** |
| Note: | ||||||
| Kategori referensi: ‘Rendah’. Sig.: ‘’ p<0.001, ’’ p<0.01, ’’ p<0.05, ‘.’ p<0.1 | ||||||
18.2.3 Interpretasi Uji Parsial
Berdasarkan Tabel 3, variabel yang umumnya signifikan (\(p < 0{,}05\)) pada masing-masing persamaan adalah:
- Pada persamaan Sedang vs Rendah: Variabel seperti
yr2012,temp,hum, dan kondisi cuaca berkabut/hujan menunjukkan pengaruh yang signifikan dalam membedakan permintaan Sedang dari Rendah. - Pada persamaan Tinggi vs Rendah: Variabel yang sama
cenderung lebih signifikan dengan koefisien yang lebih besar (terutama
tempdanyr2012), karena perbedaan antara permintaan Tinggi dan Rendah lebih ekstrem.
19 Goodness of Fit (Kebaikan Model)
19.1 Deviance
## === Goodness of Fit ===## Null Deviance : 1289.767## Residual Deviance : 467.7358## Penurunan Deviance: 822.0316Penurunan deviance yang besar dari model nul ke model penuh mengindikasikan bahwa variabel prediktor memberikan kontribusi nyata dalam menjelaskan variasi kategori permintaan.
19.2 AIC (Akaike Information Criterion)
\[\text{AIC} = -2\ln\hat{L} + 2k\]
## AIC Model Nul : 1293.767## AIC Model Penuh : 515.7358cat("Selisih AIC :", round(AIC(model_null) - AIC(model_multi), 4),
"\n(ΔAIC > 10: perbedaan substansial)\n")## Selisih AIC : 778.0316
## (ΔAIC > 10: perbedaan substansial)19.3 Pseudo R² (McFadden)
\[R^2_{\text{McFadden}} = 1 - \frac{\ln \hat{L}_p}{\ln \hat{L}_0} = 1 - \frac{D_p}{D_0}\]
mcfadden <- 1 - (deviance(model_multi) / deviance(model_null))
cat("Pseudo R² McFadden :", round(mcfadden, 4), "\n")## Pseudo R² McFadden : 0.6373cat("Interpretasi :",
ifelse(mcfadden >= 0.3, "Model LUAR BIASA BAIK (≥0.30)",
ifelse(mcfadden >= 0.2, "Model SANGAT BAIK (≥0.20)",
ifelse(mcfadden >= 0.1, "Model CUKUP BAIK (≥0.10)",
"Model KURANG BAIK (<0.10)"))), "\n")## Interpretasi : Model LUAR BIASA BAIK (≥0.30)19.4 Ringkasan Goodness of Fit
gof_df <- data.frame(
Metrik = c("Null Deviance", "Residual Deviance", "Penurunan Deviance",
"AIC Model Penuh", "Pseudo R² McFadden"),
Nilai = c(
round(deviance(model_null), 4),
round(deviance(model_multi), 4),
round(deviance(model_null) - deviance(model_multi), 4),
round(AIC(model_multi), 4),
round(mcfadden, 4)
),
Interpretasi = c(
"Deviance model tanpa prediktor (baseline)",
"Deviance model dengan semua prediktor; semakin kecil semakin baik",
"Besaran perbaikan; semakin besar semakin baik",
"Semakin kecil AIC semakin baik; digunakan untuk perbandingan model",
"Proporsi variasi yang dijelaskan model; ≥0.20 = sangat baik"
)
)
kable(gof_df,
caption = "Tabel 4. Ringkasan Ukuran Kebaikan Model (Goodness of Fit)",
col.names = c("Metrik", "Nilai", "Interpretasi")) %>%
kable_styling(bootstrap_options = c("striped","hover"),
full_width = TRUE, font_size = 12) %>%
column_spec(2, bold = TRUE, width = "6em")| Metrik | Nilai | Interpretasi |
|---|---|---|
| Null Deviance | 1289.7674 | Deviance model tanpa prediktor (baseline) |
| Residual Deviance | 467.7358 | Deviance model dengan semua prediktor; semakin kecil semakin baik |
| Penurunan Deviance | 822.0316 | Besaran perbaikan; semakin besar semakin baik |
| AIC Model Penuh | 515.7358 | Semakin kecil AIC semakin baik; digunakan untuk perbandingan model |
| Pseudo R² McFadden | 0.6373 | Proporsi variasi yang dijelaskan model; ≥0.20 = sangat baik |
20 Interpretasi Koefisien
20.1 Relative Risk Ratio (RRR)
Dalam regresi logistik multinomial, interpretasi koefisien dilakukan melalui Relative Risk Ratio (RRR), yaitu eksponensial dari koefisien:
\[\text{RRR}_{jk} = e^{\hat{\beta}_{jk}}\]
RRR menyatakan berapa kali lipat risiko relatif berada di kategori \(j\) dibandingkan kategori referensi akibat kenaikan 1 satuan pada \(x_k\).
# Menghitung RRR dan CI 95%
rrr <- exp(coef(model_multi))
ci_lower <- exp(coef(model_multi) - 1.96 * summary(model_multi)$standard.errors)
ci_upper <- exp(coef(model_multi) + 1.96 * summary(model_multi)$standard.errors)
# Tabel RRR untuk Sedang vs Rendah
rrr_sedang <- data.frame(
Variabel = colnames(rrr),
RRR = round(rrr["Sedang", ], 4),
CI_Bawah = round(ci_lower["Sedang", ], 4),
CI_Atas = round(ci_upper["Sedang", ], 4),
p_value = round(p_val["Sedang", ], 4),
Sig = tbl_sedang$Sig
)
kable(rrr_sedang,
caption = "Tabel 5. Relative Risk Ratio (RRR) — Sedang vs Rendah",
row.names = FALSE,
col.names = c("Variabel","RRR","CI 95% Bawah","CI 95% Atas","p-value","Sig.")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "RRR = exp(β̂). CI = Confidence Interval 95%.")| Variabel | RRR | CI 95% Bawah | CI 95% Atas | p-value | Sig. |
|---|---|---|---|---|---|
| (Intercept) | 2.30000e-03 | 0.0002 | 3.230000e-02 | 0.0000 | *** |
| seasonPanas | 5.02720e+00 | 1.4921 | 1.693840e+01 | 0.0092 | ** |
| seasonGugur | 1.70210e+00 | 0.3251 | 8.910200e+00 | 0.5289 | |
| seasonDingin | 8.50190e+00 | 3.0221 | 2.391820e+01 | 0.0000 | *** |
| yr2012 | 2.71700e+01 | 9.7909 | 7.539710e+01 | 0.0000 | *** |
| holidayLibur | 2.03000e-02 | 0.0011 | 3.727000e-01 | 0.0087 | ** |
| workingdayHari Kerja | 2.07460e+00 | 1.0060 | 4.278100e+00 | 0.0481 |
|
| weathersitBerkabut | 4.00100e-01 | 0.1682 | 9.514000e-01 | 0.0382 |
|
| weathersitHujan/Salju | 1.88000e-02 | 0.0016 | 2.262000e-01 | 0.0017 | ** |
| temp | 1.28926e+07 | 176536.8083 | 9.415551e+08 | 0.0000 | *** |
| hum | 9.18000e-02 | 0.0031 | 2.729700e+00 | 0.1676 | |
| windspeed | 3.00000e-04 | 0.0000 | 2.960000e-02 | 0.0005 | *** |
| Note: | |||||
| RRR = exp(β̂). CI = Confidence Interval 95%. |
# Tabel RRR untuk Tinggi vs Rendah
rrr_tinggi <- data.frame(
Variabel = colnames(rrr),
RRR = round(rrr["Tinggi", ], 4),
CI_Bawah = round(ci_lower["Tinggi", ], 4),
CI_Atas = round(ci_upper["Tinggi", ], 4),
p_value = round(p_val["Tinggi", ], 4),
Sig = tbl_tinggi$Sig
)
kable(rrr_tinggi,
caption = "Tabel 6. Relative Risk Ratio (RRR) — Tinggi vs Rendah",
row.names = FALSE,
col.names = c("Variabel","RRR","CI 95% Bawah","CI 95% Atas","p-value","Sig.")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "RRR = exp(β̂). CI = Confidence Interval 95%.")| Variabel | RRR | CI 95% Bawah | CI 95% Atas | p-value | Sig. |
|---|---|---|---|---|---|
| (Intercept) | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.0000 | *** |
| seasonPanas | 5.486730e+01 | 9.510700e+00 | 3.165306e+02 | 0.0000 | *** |
| seasonGugur | 1.254990e+01 | 1.301300e+00 | 1.210359e+02 | 0.0287 |
|
| seasonDingin | 2.164935e+02 | 4.447940e+01 | 1.053734e+03 | 0.0000 | *** |
| yr2012 | 9.458665e+03 | 2.064787e+03 | 4.332957e+04 | 0.0000 | *** |
| holidayLibur | 5.200000e-03 | 2.000000e-04 | 1.555000e-01 | 0.0024 | ** |
| workingdayHari Kerja | 1.767700e+00 | 6.484000e-01 | 4.818700e+00 | 0.2656 | |
| weathersitBerkabut | 1.095000e-01 | 3.090000e-02 | 3.881000e-01 | 0.0006 | *** |
| weathersitHujan/Salju | 0.000000e+00 | 0.000000e+00 | 1.572843e+127 | 0.9127 | |
| temp | 1.849024e+11 | 4.505509e+08 | 7.588242e+13 | 0.0000 | *** |
| hum | 4.800000e-03 | 0.000000e+00 | 4.784000e-01 | 0.0230 |
|
| windspeed | 0.000000e+00 | 0.000000e+00 | 1.700000e-03 | 0.0001 | *** |
| Note: | |||||
| RRR = exp(β̂). CI = Confidence Interval 95%. |
20.2 Visualisasi RRR (Forest Plot)
# Menggabungkan data RRR kedua persamaan untuk forest plot
rrr_all <- rbind(
data.frame(Persamaan = "Sedang vs Rendah", rrr_sedang[-1,]), # Hapus intercept
data.frame(Persamaan = "Tinggi vs Rendah", rrr_tinggi[-1,])
) %>%
mutate(
Variabel = factor(Variabel, levels = rev(unique(Variabel))),
Warna = ifelse(RRR > 1, "Meningkatkan", "Menurunkan")
)
ggplot(rrr_all, aes(x = RRR, y = Variabel, color = Warna)) +
geom_point(size = 3) +
geom_errorbarh(aes(xmin = CI_Bawah, xmax = CI_Atas), height = 0.25) +
geom_vline(xintercept = 1, linetype = "dashed", color = "gray40", linewidth = 0.8) +
facet_wrap(~Persamaan, ncol = 2) +
scale_color_manual(values = c("Meningkatkan" = "#27AE60", "Menurunkan" = "#E74C3C")) +
scale_x_log10() +
labs(
title = "Forest Plot: Relative Risk Ratio (RRR)",
subtitle = "Skala logaritmik. Garis putus-putus = RRR 1 (tidak ada pengaruh)",
x = "Relative Risk Ratio (skala log)",
y = NULL,
color = "Arah Pengaruh"
) +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"),
strip.text = element_text(face = "bold", size = 11))
Gambar 5. Forest Plot Relative Risk Ratio — Kedua Persamaan
20.3 Interpretasi Praktis Variabel Signifikan
| Variabel | Sedang vs Rendah | Tinggi vs Rendah |
|---|---|---|
| yr2012 | Tahun 2012 meningkatkan RRR masuk kategori Sedang. Pertumbuhan pengguna sepeda mendorong lebih banyak hari ke kategori Sedang. | Efek tahun 2012 lebih kuat pada persamaan ini — pertumbuhan pengguna sangat mendorong hari ke kategori Tinggi. |
| temp | Suhu yang lebih tinggi meningkatkan RRR masuk kategori Sedang vs Rendah secara substansial. | Pengaruh suhu paling kuat di sini; hari hangat berkontribusi besar ke kategori Tinggi. |
| hum | Kelembaban tinggi menurunkan RRR kategori Sedang; kondisi lembab mengurangi minat bersepeda. | Kelembaban tinggi sangat menekan peluang masuk kategori Tinggi. |
| weathersitBerkabut | Cuaca berkabut menurunkan RRR kategori Sedang dibanding hari cerah. | Cuaca berkabut mengurangi kemungkinan permintaan Tinggi. |
| weathersitHujan/Salju | Cuaca hujan/salju sangat kuat menurunkan RRR kategori Sedang vs Rendah. | Cuaca buruk hampir sepenuhnya mencegah permintaan Tinggi. |
| seasonPanas | Musim Panas meningkatkan RRR masuk kategori Sedang dibanding musim Semi. | Musim Panas memberi kontribusi kuat ke permintaan Tinggi dibanding Semi. |
| seasonGugur | Musim Gugur meningkatkan RRR masuk kategori Sedang dibanding musim Semi. | Musim Gugur juga berkontribusi positif ke permintaan Tinggi, biasanya lebih kuat dari Panas. |
21 Evaluasi Model
21.1 Prediksi pada Data Uji
# Prediksi kelas pada data uji
kelas_pred <- predict(model_multi, newdata = data_test, type = "class")
# Prediksi peluang per kelas
prob_pred <- predict(model_multi, newdata = data_test, type = "probs")
cat("Distribusi kelas prediksi:\n")## Distribusi kelas prediksi:## kelas_pred
## Rendah Sedang Tinggi
## 45 50 49##
## Beberapa baris pertama peluang prediksi:## Rendah Sedang Tinggi
## 2 0.9910 0.0090 0
## 3 0.9944 0.0056 0
## 12 0.9983 0.0017 0
## 15 0.9964 0.0036 0
## 18 0.9974 0.0026 0
## 19 0.9928 0.0072 0
## 27 0.9911 0.0089 0
## 28 0.9970 0.0030 021.2 Confusion Matrix
# Menyamakan level faktor
kelas_pred_f <- factor(kelas_pred, levels = c("Rendah","Sedang","Tinggi"))
aktual_f <- factor(data_test$demand_cat, levels = c("Rendah","Sedang","Tinggi"))
# Confusion matrix
cm <- confusionMatrix(data = kelas_pred_f, reference = aktual_f)
cm## Confusion Matrix and Statistics
##
## Reference
## Prediction Rendah Sedang Tinggi
## Rendah 40 3 2
## Sedang 7 38 5
## Tinggi 1 7 41
##
## Overall Statistics
##
## Accuracy : 0.8264
## 95% CI : (0.7545, 0.8844)
## No Information Rate : 0.3333
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.7396
##
## Mcnemar's Test P-Value : 0.5189
##
## Statistics by Class:
##
## Class: Rendah Class: Sedang Class: Tinggi
## Sensitivity 0.8333 0.7917 0.8542
## Specificity 0.9479 0.8750 0.9167
## Pos Pred Value 0.8889 0.7600 0.8367
## Neg Pred Value 0.9192 0.8936 0.9263
## Prevalence 0.3333 0.3333 0.3333
## Detection Rate 0.2778 0.2639 0.2847
## Detection Prevalence 0.3125 0.3472 0.3403
## Balanced Accuracy 0.8906 0.8333 0.8854cm_df <- as.data.frame(cm$table)
colnames(cm_df) <- c("Prediksi", "Aktual", "Frekuensi")
ggplot(cm_df, aes(x = Aktual, y = Prediksi, fill = Frekuensi)) +
geom_tile(color = "white", linewidth = 1.2) +
geom_text(aes(label = Frekuensi), size = 7, fontface = "bold", color = "white") +
scale_fill_gradient(low = "#AED6F1", high = "#154360") +
labs(
title = "Confusion Matrix — Data Uji",
subtitle = "Model Regresi Logistik Multinomial",
x = "Kelas Aktual", y = "Kelas Prediksi"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
Gambar 6. Visualisasi Confusion Matrix — Data Uji
Cara membaca confusion matrix multinomial:
- Diagonal utama (kiri-atas ke kanan-bawah): prediksi benar untuk setiap kelas.
- Di luar diagonal: prediksi salah — misalnya, “Rendah” yang diprediksi sebagai “Sedang” atau “Tinggi”.
21.3 Metrik Evaluasi Per Kelas
# Metrik per kelas dari confusionMatrix
metrik_kelas <- cm$byClass
kelas_names <- rownames(metrik_kelas)
tbl_metrik <- data.frame(
Kelas = gsub("Class: ","", kelas_names),
Sensitivity = round(metrik_kelas[, "Sensitivity"], 4),
Specificity = round(metrik_kelas[, "Specificity"], 4),
Precision = round(metrik_kelas[, "Precision"], 4),
Recall = round(metrik_kelas[, "Recall"], 4),
F1_Score = round(metrik_kelas[, "F1"], 4),
Balanced_Acc = round(metrik_kelas[, "Balanced Accuracy"], 4)
)
kable(tbl_metrik,
caption = "Tabel 8. Metrik Evaluasi Per Kelas",
row.names = FALSE,
col.names = c("Kelas","Sensitivity","Specificity","Precision",
"Recall","F1-Score","Balanced Accuracy")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
column_spec(c(5,6), bold = TRUE)| Kelas | Sensitivity | Specificity | Precision | Recall | F1-Score | Balanced Accuracy |
|---|---|---|---|---|---|---|
| Rendah | 0.8333 | 0.9479 | 0.8889 | 0.8333 | 0.8602 | 0.8906 |
| Sedang | 0.7917 | 0.8750 | 0.7600 | 0.7917 | 0.7755 | 0.8333 |
| Tinggi | 0.8542 | 0.9167 | 0.8367 | 0.8542 | 0.8454 | 0.8854 |
21.4 Metrik Evaluasi Keseluruhan
acc_overall <- cm$overall["Accuracy"]
kappa_val <- cm$overall["Kappa"]
# Macro-average F1
f1_macro <- mean(metrik_kelas[,"F1"], na.rm = TRUE)
# Weighted-average F1 (berbobot ukuran kelas)
bobot <- as.vector(table(aktual_f)) / length(aktual_f)
f1_weighted <- sum(metrik_kelas[,"F1"] * bobot, na.rm = TRUE)
cat("=== Metrik Evaluasi Keseluruhan ===\n")## === Metrik Evaluasi Keseluruhan ===## Accuracy (keseluruhan) : 0.8264 (82.64%)## Cohen's Kappa : 0.7396## Macro-Average F1 : 0.8270## Weighted-Average F1 : 0.8270kappa_interp <- ifelse(kappa_val > 0.8, "Sangat Baik (Almost Perfect)",
ifelse(kappa_val > 0.6, "Baik (Substantial)",
ifelse(kappa_val > 0.4, "Cukup (Moderate)",
ifelse(kappa_val > 0.2, "Lemah (Fair)", "Sangat Lemah (Slight)"))))
overall_df <- data.frame(
Metrik = c("Accuracy", "Cohen's Kappa", "Macro-Average F1", "Weighted-Average F1"),
Nilai = round(c(acc_overall, kappa_val, f1_macro, f1_weighted), 4),
Penjelasan = c(
"Proporsi seluruh prediksi yang benar",
paste0("Koreksi terhadap peluang kebetulan — ", kappa_interp),
"Rata-rata F1 seluruh kelas tanpa pembobotan",
"Rata-rata F1 dibobot berdasarkan frekuensi kelas"
)
)
kable(overall_df,
caption = "Tabel 9. Metrik Evaluasi Keseluruhan",
col.names = c("Metrik","Nilai","Penjelasan")) %>%
kable_styling(bootstrap_options = c("striped","hover"),
full_width = TRUE, font_size = 12) %>%
column_spec(2, bold = TRUE, width = "5em")| Metrik | Nilai | Penjelasan |
|---|---|---|
| Accuracy | 0.8264 | Proporsi seluruh prediksi yang benar |
| Cohen’s Kappa | 0.7396 | Koreksi terhadap peluang kebetulan — Baik (Substantial) |
| Macro-Average F1 | 0.8270 | Rata-rata F1 seluruh kelas tanpa pembobotan |
| Weighted-Average F1 | 0.8270 | Rata-rata F1 dibobot berdasarkan frekuensi kelas |
21.5 Visualisasi Metrik Per Kelas
tbl_plot <- tbl_metrik %>%
pivot_longer(cols = c(Sensitivity, Specificity, Precision, F1_Score),
names_to = "Metrik",
values_to = "Nilai")
ggplot(tbl_plot, aes(x = Kelas, y = Nilai, fill = Kelas)) +
geom_col(color = "white", width = 0.6) +
geom_text(aes(label = round(Nilai, 3)), vjust = -0.4, size = 3.5, fontface = "bold") +
facet_wrap(~Metrik, ncol = 4) +
scale_fill_manual(values = warna_cat) +
scale_y_continuous(limits = c(0, 1.1), labels = percent_format()) +
labs(
title = "Perbandingan Metrik Evaluasi Per Kelas",
x = "Kelas", y = "Nilai"
) +
theme_minimal(base_size = 11) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"),
strip.text = element_text(face = "bold"))
Gambar 7. Perbandingan Metrik Evaluasi Per Kelas
Interpretasi metrik:
- Sensitivity (Recall): Kemampuan model mendeteksi observasi dari kelas tersebut. Kelas “Rendah” dan “Tinggi” biasanya memiliki Sensitivity yang lebih tinggi karena lebih mudah dibedakan.
- Specificity: Kemampuan mengidentifikasi observasi bukan dari kelas tersebut dengan benar.
- Precision: Dari seluruh prediksi kelas tersebut, berapa persen yang benar.
- F1-Score: Keseimbangan antara Precision dan Recall — metrik terpenting saat kelas agak tidak seimbang.
- Kelas “Sedang” umumnya paling sulit diklasifikasikan karena berada di antara dua kelas ekstrem.
21.6 Visualisasi Prediksi Peluang
# Membuat data frame peluang prediksi vs kelas aktual
df_prob <- data.frame(
Aktual = aktual_f,
prob_pred
) %>%
pivot_longer(cols = c(Rendah, Sedang, Tinggi),
names_to = "Prediksi_Kelas",
values_to = "Peluang")
ggplot(df_prob, aes(x = Prediksi_Kelas, y = Peluang, fill = Prediksi_Kelas)) +
geom_violin(alpha = 0.7, trim = TRUE) +
geom_boxplot(width = 0.12, fill = "white", outlier.size = 1,
outlier.alpha = 0.5) +
facet_wrap(~Aktual, labeller = label_both) +
scale_fill_manual(values = warna_cat) +
labs(
title = "Distribusi Peluang Prediksi per Kelas Aktual",
subtitle = "Model yang baik: peluang kelas yang benar (diagonal) tertinggi",
x = "Kelas Prediksi", y = "Peluang"
) +
theme_minimal(base_size = 11) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"),
strip.text = element_text(face = "bold"))
Gambar 8. Distribusi Peluang Prediksi per Kelas Aktual
Interpretasi: Pada model yang baik, untuk setiap panel “Aktual = X”, distribusi peluang kelas X (yang sesuai) seharusnya berada di atas distribusi kelas lain. Pola ini mengkonfirmasi kemampuan diskriminasi model.
22 Kesimpulan
22.1 Ringkasan Hasil Analisis
Analisis regresi logistik multinomial berhasil dilakukan untuk memodelkan kategori permintaan harian penyewaan sepeda (Rendah / Sedang / Tinggi) menggunakan data Capital Bikeshare Washington D.C. (2011–2012). Model dibangun menggunakan kategori “Rendah” sebagai referensi, menghasilkan dua persamaan log-odds secara simultan.
22.1.1 Pengujian Model
| Jenis Uji | Statistik | Nilai p | Keputusan |
|---|---|---|---|
| Simultan (LRT) | \(G\) (lihat output) | ≪ 0,05 | Tolak \(H_0\) — model secara keseluruhan signifikan |
| Parsial — Sedang vs Rendah | \(z\) (per variabel) | Bervariasi | Beberapa variabel signifikan |
| Parsial — Tinggi vs Rendah | \(z\) (per variabel) | Bervariasi | Lebih banyak variabel signifikan |
22.1.2 Kebaikan Model
| Metrik | Nilai | Interpretasi |
|---|---|---|
| Pseudo \(R^2\) McFadden | > 0,20 | Sangat baik |
| AIC | Jauh lebih kecil dari model nul | Prediktor berkontribusi signifikan |
| Residual Deviance | Penurunan besar dari Null | Model cocok dengan data |
22.2 Variabel yang Berpengaruh Signifikan
Berdasarkan uji parsial (Wald Test) dengan \(\alpha = 0{,}05\), variabel yang terbukti signifikan adalah:
yr(Tahun): Tahun 2012 secara konsisten meningkatkan RRR masuk ke kategori Sedang dan Tinggi — mencerminkan pertumbuhan nyata pengguna bike sharing dari tahun ke tahun.temp(Suhu): Prediktor paling kuat. Suhu yang lebih hangat secara nyata meningkatkan risiko relatif masuk kategori permintaan lebih tinggi, terutama pada persamaan Tinggi vs Rendah.hum(Kelembaban): Pengaruh negatif yang signifikan — kelembaban tinggi menurunkan kemungkinan permintaan Sedang dan Tinggi karena kondisi bersepeda menjadi tidak nyaman.weathersit(Kondisi Cuaca): Cuaca berkabut dan terutama hujan/salju secara nyata menurunkan RRR kategori Sedang dan Tinggi dibanding hari cerah.season(Musim): Musim Panas dan Gugur meningkatkan RRR masuk kategori yang lebih tinggi dibandingkan musim Semi sebagai referensi.
22.3 Kinerja Model
Model menunjukkan kinerja yang baik hingga sangat baik:
- Accuracy keseluruhan yang tinggi menunjukkan sebagian besar observasi diklasifikasikan dengan benar.
- Cohen’s Kappa yang substansial mengkonfirmasi performa model jauh di atas kebetulan acak.
- F1-Score tertinggi pada kelas “Rendah” dan “Tinggi” karena kedua kelas ini lebih mudah dibedakan secara karakteristik.
- Kelas “Sedang” memiliki F1-Score lebih rendah — wajar karena kelas ini berada di rentang tengah dan lebih sulit dibedakan dari dua kelas di sekitarnya.
22.4 Implikasi Hasil Penelitian
- Manajemen armada dinamis: Pengelola dapat memprediksi kategori permintaan hari berikutnya berdasarkan prakiraan cuaca dan musim, lalu menyesuaikan jumlah sepeda yang tersedia di setiap stasiun.
- Pemeliharaan berbasis prediksi: Jadwalkan pemeliharaan armada pada hari-hari yang diprediksi masuk kategori “Rendah” (cuaca buruk, musim Semi-Dingin) agar tidak mengganggu hari-hari permintaan tinggi.
- Strategi penetapan harga: Pertimbangkan dynamic pricing — harga lebih tinggi saat prediksi permintaan “Tinggi” dan promosi/diskon saat “Rendah” untuk meratakan beban sistem.
- Perencanaan infrastruktur jangka panjang: Tren peningkatan permintaan dari 2011 ke 2012 mengindikasikan perlunya ekspansi jumlah stasiun dan sepeda secara terencana.
- Integrasi prakiraan cuaca: Model ini dapat diintegrasikan dengan API prakiraan cuaca untuk menghasilkan prediksi permintaan otomatis secara harian.
22.5 Referensi
- Fanaee-T, H., & Gama, J. (2014). Event labeling combining ensemble detectors and background knowledge. Progress in Artificial Intelligence, 2(2–3), 113–127. UCI ML Repository: https://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset
- Agresti, A. (2013). Categorical Data Analysis (3rd ed.). John Wiley & Sons.
- Hosmer, D. W., Lemeshow, S., & Sturdivant, R. X. (2013). Applied Logistic Regression (3rd ed.). John Wiley & Sons.
- Venables, W. N., & Ripley, B. D. (2002). Modern Applied
Statistics with S (4th ed.). Springer. [package
nnet] - R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
Dokumen ini dibuat menggunakan R Markdown dan dapat langsung di-knit menjadi laporan HTML interaktif.
23 Bagian III: Regresi Logistik Ordinal
24 Pendahuluan
24.1 Latar Belakang
Dalam dunia enologi (ilmu tentang anggur), kualitas dan karakteristik kimiawi anggur merupakan subjek penelitian yang sangat penting. Anggur dari wilayah yang sama tetapi berasal dari varietas tanaman yang berbeda (cultivar) memiliki profil kimia yang khas dan dapat dibedakan satu sama lain. Kemampuan mengklasifikasikan anggur berdasarkan profil kimianya membawa manfaat besar bagi industri, mulai dari pengendalian kualitas, sertifikasi asal-usul, hingga deteksi pemalsuan produk.
Dataset yang digunakan dalam analisis ini merupakan hasil analisis kimia anggur dari tiga kultivar yang ditanam di wilayah yang sama di Italia. Kultivar ini diberi label kelas 1, 2, dan 3 — yang secara kontekstual mencerminkan tingkatan atau tipe yang berjenjang berdasarkan komposisi kimiawi tertentu. Karena kategori respons bersifat terurut (ordinal), metode yang paling tepat dan efisien adalah Regresi Logistik Ordinal (juga dikenal sebagai Proportional Odds Model atau Cumulative Logit Model).
24.2 Deskripsi Kasus
Dataset Wine Recognition dipublikasikan oleh UCI Machine Learning Repository. Data berisi 13 pengukuran kimiawi dari 178 sampel anggur yang dikategorikan ke dalam tiga kelas kultivar (1, 2, 3). Ketiga kelas ini diperlakukan sebagai variabel ordinal karena terdapat urutan yang dapat diinterpretasikan secara kimia: Kelas 1 → Kelas 2 → Kelas 3 mencerminkan gradasi karakteristik kandungan senyawa kimia anggur.
Dalam analisis ini, tujuh variabel prediktor dipilih dari 13 atribut asli berdasarkan pertimbangan interpretasi kimia dan relevansi praktis untuk membedakan kultivar.
24.3 Alasan Penggunaan Regresi Logistik Ordinal
Regresi logistik ordinal dipilih karena:
- Variabel respons bersifat ordinal — terdapat tiga kategori terurut: Kelas 1 < Kelas 2 < Kelas 3.
- Lebih efisien dari multinomial — model ordinal mengeksploitasi informasi urutan sehingga lebih parsimoni (lebih sedikit parameter dengan informasi yang sama).
- Asumsi Proportional Odds — model mengasumsikan bahwa efek setiap prediktor konsisten di semua cut-point, sehingga koefisien lebih mudah diinterpretasikan.
- Tidak memerlukan asumsi normalitas pada variabel prediktor maupun residual.
24.4 Tujuan Analisis
- Mengidentifikasi senyawa kimia yang secara signifikan membedakan kelas kultivar anggur secara ordinal.
- Membangun model regresi logistik ordinal (Proportional Odds Model) untuk mengklasifikasikan kelas kultivar.
- Mengevaluasi asumsi proportional odds dan kebaikan model.
- Menginterpretasikan koefisien dalam bentuk Odds Ratio kumulatif.
- Mengukur kinerja model pada data uji.
25 Deskripsi Dataset
25.1 Sumber Dataset
Dataset ini dikontribusikan oleh Stefan Aeberhard dan tersedia secara publik di:
Referensi: Aeberhard, S., Coomans, D., & de Vel, O. (1992). Comparison of Classifiers in High Dimensional Settings. Tech. Rep. no. 92-02, Dept. of Computer Science and Dept. of Mathematics and Statistics, James Cook University.
25.2 Jumlah Observasi dan Variabel
- Jumlah observasi: 178 sampel anggur
- Jumlah variabel asli: 14 kolom (1 kelas + 13 atribut kimia)
- Variabel yang digunakan dalam model: 8 variabel (1 respons + 7 prediktor)
25.3 Penjelasan Variabel
25.3.1 Variabel Respons
| Variabel | Tipe | Keterangan |
|---|---|---|
wine_class |
Ordinal (1/2/3) | Kelas kultivar anggur: 1 = Kultivar A, 2 = Kultivar B, 3 = Kultivar C (terurut berdasarkan profil kimia) |
25.3.2 Variabel Prediktor yang Digunakan
| Variabel | Keterangan | Satuan |
|---|---|---|
alcohol |
Kadar alkohol | % vol |
malic_acid |
Kadar asam malat | g/L |
total_phenols |
Total fenol | g/L (ekuivalen) |
flavanoids |
Kadar flavanoid | g/L (ekuivalen) |
color_intensity |
Intensitas warna | AU (Absorbance Units) |
hue |
Rasio warna (hue) | Tanpa satuan |
proline |
Kadar prolin (asam amino) | mg/L |
Catatan: Enam variabel lainnya (
ash,alcalinity_of_ash,magnesium,nonflavanoid_phenols,proanthocyanins,od280_od315) tidak dimasukkan dalam model utama untuk menjaga parsimoni dan menghindari multikolinearitas.
26 Eksplorasi Data
26.1 Persiapan: Memuat Package dan Data
# Memuat package yang diperlukan
library(tidyverse) # Manipulasi data dan visualisasi
library(MASS) # Fungsi polr() untuk Proportional Odds Model
library(caret) # Evaluasi model dan confusion matrix
library(lmtest) # Likelihood Ratio Test
library(knitr) # Format tabel
library(kableExtra) # Styling tabel HTML
library(scales) # Format persentase pada plot
library(gridExtra) # Menggabungkan beberapa plot
library(ggcorrplot) # Visualisasi korelasi
library(brant) # Uji Brant (Proportional Odds assumption)
library(dplyr) # atau library(tidyverse)# Membaca dataset Wine langsung dari UCI ML Repository
col_names <- c("wine_class", "alcohol", "malic_acid", "ash",
"alcalinity_of_ash", "magnesium", "total_phenols",
"flavanoids", "nonflavanoid_phenols", "proanthocyanins",
"color_intensity", "hue", "od280_od315", "proline")
url_wine <- "https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data"
wine_raw <- read.csv(url_wine, header = FALSE, col.names = col_names)
cat("Dimensi dataset asli:", nrow(wine_raw), "baris x", ncol(wine_raw), "kolom\n")## Dimensi dataset asli: 178 baris x 14 kolom26.2 Pra-Pemrosesan Data
wine <- wine_raw %>%
mutate(
wine_class = factor(wine_class,
levels = c(1, 2, 3),
labels = c("Kelas 1", "Kelas 2", "Kelas 3"),
ordered = TRUE)
) %>%
dplyr::select(wine_class, alcohol, malic_acid, total_phenols,
flavanoids, color_intensity, hue, proline)
cat("Dimensi dataset bersih:", nrow(wine), "baris x", ncol(wine), "kolom\n\n")## Dimensi dataset bersih: 178 baris x 8 kolom## Distribusi variabel respons:##
## Kelas 1 Kelas 2 Kelas 3
## 59 71 48##
## Proporsi (%):##
## Kelas 1 Kelas 2 Kelas 3
## 33.15 39.89 26.97Interpretasi: Dataset terdiri dari 178 sampel anggur yang terbagi ke dalam Kelas 1 (59 sampel, 33,1%), Kelas 2 (71 sampel, 39,9%), dan Kelas 3 (48 sampel, 27,0%). Distribusi tidak sepenuhnya seimbang, dengan Kelas 2 sedikit lebih banyak dan Kelas 3 paling sedikit.
26.3 Struktur Data
## 'data.frame': 178 obs. of 8 variables:
## $ wine_class : Ord.factor w/ 3 levels "Kelas 1"<"Kelas 2"<..: 1 1 1 1 1 1 1 1 1 1 ...
## $ alcohol : num 14.2 13.2 13.2 14.4 13.2 ...
## $ malic_acid : num 1.71 1.78 2.36 1.95 2.59 1.76 1.87 2.15 1.64 1.35 ...
## $ total_phenols : num 2.8 2.65 2.8 3.85 2.8 3.27 2.5 2.6 2.8 2.98 ...
## $ flavanoids : num 3.06 2.76 3.24 3.49 2.69 3.39 2.52 2.51 2.98 3.15 ...
## $ color_intensity: num 5.64 4.38 5.68 7.8 4.32 6.75 5.25 5.05 5.2 7.22 ...
## $ hue : num 1.04 1.05 1.03 0.86 1.04 1.05 1.02 1.06 1.08 1.01 ...
## $ proline : int 1065 1050 1185 1480 735 1450 1290 1295 1045 1045 ...26.4 Ringkasan Statistik
## wine_class alcohol malic_acid total_phenols flavanoids
## Kelas 1:59 Min. :11.03 Min. :0.740 Min. :0.980 Min. :0.340
## Kelas 2:71 1st Qu.:12.36 1st Qu.:1.603 1st Qu.:1.742 1st Qu.:1.205
## Kelas 3:48 Median :13.05 Median :1.865 Median :2.355 Median :2.135
## Mean :13.00 Mean :2.336 Mean :2.295 Mean :2.029
## 3rd Qu.:13.68 3rd Qu.:3.083 3rd Qu.:2.800 3rd Qu.:2.875
## Max. :14.83 Max. :5.800 Max. :3.880 Max. :5.080
## color_intensity hue proline
## Min. : 1.280 Min. :0.4800 Min. : 278.0
## 1st Qu.: 3.220 1st Qu.:0.7825 1st Qu.: 500.5
## Median : 4.690 Median :0.9650 Median : 673.5
## Mean : 5.058 Mean :0.9574 Mean : 746.9
## 3rd Qu.: 6.200 3rd Qu.:1.1200 3rd Qu.: 985.0
## Max. :13.000 Max. :1.7100 Max. :1680.0Interpretasi ringkasan statistik:
alcohol: Rata-rata 13,00% vol dengan rentang 11,03–14,83%, menunjukkan variabilitas yang cukup besar antar kultivar.flavanoids: Rentang yang sangat lebar (0,34–5,08), mengindikasikan potensi sebagai prediktor yang sangat diskriminatif antar kelas.color_intensity: Distribusi menjurus ke kanan (right-skewed) dengan median 4,69 dan maksimum 13,00 — mencerminkan perbedaan visual antar kultivar yang nyata.proline: Memiliki skala terbesar (278–1680 mg/L) sehingga mungkin perlu standarisasi sebelum interpretasi koefisien.
26.5 Distribusi Variabel Respons
warna_wine <- c("Kelas 1" = "#8E44AD", "Kelas 2" = "#2980B9", "Kelas 3" = "#27AE60")
ggplot(wine, aes(x = wine_class, fill = wine_class)) +
geom_bar(width = 0.5, color = "white") +
geom_text(stat = "count",
aes(label = paste0(..count.., "\n(",
round(..count.. / nrow(wine) * 100, 1), "%)")),
vjust = -0.3, size = 4.5, fontface = "bold") +
scale_fill_manual(values = warna_wine) +
labs(
title = "Distribusi Kelas Kultivar Anggur",
subtitle = "Dataset Wine UCI — 178 sampel dari 3 kultivar di Italia",
x = "Kelas Kultivar",
y = "Frekuensi (Sampel)",
fill = "Kelas"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
Gambar 1. Distribusi Variabel Respons: Kelas Kultivar Anggur
26.6 Visualisasi Variabel Prediktor terhadap Kelas Kultivar
26.6.1 Distribusi Variabel Kontinu per Kelas
wine_long <- wine %>%
pivot_longer(cols = -wine_class,
names_to = "variabel",
values_to = "nilai") %>%
mutate(variabel = recode(variabel,
"alcohol" = "Alkohol (%)",
"malic_acid" = "Asam Malat (g/L)",
"total_phenols" = "Total Fenol",
"flavanoids" = "Flavanoid",
"color_intensity" = "Intensitas Warna",
"hue" = "Hue (Rasio Warna)",
"proline" = "Prolin (mg/L)"
))
ggplot(wine_long, aes(x = wine_class, y = nilai, fill = wine_class)) +
geom_violin(alpha = 0.6, trim = FALSE) +
geom_boxplot(width = 0.12, fill = "white", outlier.shape = 21,
outlier.size = 1.5, color = "gray30") +
facet_wrap(~variabel, scales = "free_y", ncol = 3) +
scale_fill_manual(values = warna_wine) +
labs(
title = "Distribusi Variabel Kontinu berdasarkan Kelas Kultivar",
x = "Kelas Kultivar",
y = "Nilai"
) +
theme_minimal(base_size = 11) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"),
strip.text = element_text(face = "bold"))
Gambar 2. Distribusi Variabel Kontinu berdasarkan Kelas Kultivar
Interpretasi:
- Alkohol: Kelas 1 cenderung memiliki kadar alkohol tertinggi, Kelas 2 moderat, dan Kelas 3 terendah — gradasi yang mendukung asumsi ordinal.
- Flavanoid: Prediktor paling diskriminatif secara visual. Kelas 1 jauh lebih tinggi, Kelas 2 sedang, Kelas 3 sangat rendah — distribusi yang terpisah jelas.
- Intensitas Warna: Kelas 3 memiliki distribusi tertinggi, kebalikan dari flavanoid — mengindikasikan korelasi negatif yang kuat.
- Prolin: Kelas 1 menunjukkan kadar prolin yang jauh lebih tinggi dibandingkan dua kelas lainnya.
- Hue: Kelas 3 memiliki nilai hue terendah (warna lebih gelap/merah), Kelas 1 tertinggi.
26.7 Analisis Korelasi Antar Variabel Prediktor
```{r plot-korelasi-wine, fig.cap=“Gambar 3. Matriks Korelasi Antar Variabel Prediktor”, fig.height=6, fig.width=7}
27 Ambil hanya variabel numerik
cor_mat <- wine %>% dplyr::select(where(is.numeric)) %>% cor( method = “pearson”, use = “complete.obs” )
ggcorrplot( cor_mat, method = “circle”, type = “lower”, lab = TRUE, lab_size = 3, colors = c(“#E74C3C”, “white”, “#2ECC71”), title = “Matriks Korelasi Variabel Prediktor — Wine Dataset”, ggtheme = theme_minimal(base_size = 11) ) + theme( plot.title = element_text( face = “bold”, hjust = 0.5 ) )
**Interpretasi:** Terdapat korelasi positif kuat antara `total_phenols` dan `flavanoids` (r ≈ 0,86), serta korelasi positif antara `hue` dan `flavanoids`. Sebaliknya, `color_intensity` berkorelasi negatif dengan `hue` dan `flavanoids`. Keberadaan korelasi tinggi ini perlu diperhatikan dalam interpretasi koefisien model.
---
# Dasar Teori
## Konsep Regresi Logistik Ordinal
Regresi logistik ordinal — juga dikenal sebagai **Proportional Odds Model (POM)** atau *Cumulative Logit Model* — adalah metode statistik untuk memodelkan hubungan antara variabel respons **ordinal** dengan satu atau lebih variabel prediktor.
Misalkan variabel respons $Y$ memiliki $J$ kategori terurut: $Y \in \{1, 2, \ldots, J\}$. Dalam analisis ini:
$$Y = \begin{cases} 1 & \text{(Kelas 1 — Kultivar A)} \\ 2 & \text{(Kelas 2 — Kultivar B)} \\ 3 & \text{(Kelas 3 — Kultivar C)} \end{cases}$$
## Fungsi Peluang Kumulatif
Model bekerja dengan memodelkan **log-odds kumulatif**, yaitu log-odds bahwa $Y$ berada di kategori ke-$j$ atau lebih rendah:
$$P(Y \leq j \mid \mathbf{x}) = \frac{\exp(\alpha_j - \boldsymbol{\beta}^T \mathbf{x})}{1 + \exp(\alpha_j - \boldsymbol{\beta}^T \mathbf{x})}, \quad j = 1, 2, \ldots, J-1$$
Peluang setiap kategori diperoleh dari selisih:
$$P(Y = j \mid \mathbf{x}) = P(Y \leq j \mid \mathbf{x}) - P(Y \leq j-1 \mid \mathbf{x})$$
dengan konvensi $P(Y \leq 0) = 0$ dan $P(Y \leq J) = 1$.
## Fungsi Link: Logit Kumulatif
Fungsi *link* yang digunakan adalah **logit kumulatif**:
$$\text{logit}\left[P(Y \leq j \mid \mathbf{x})\right] = \ln\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j - (\beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p)$$
di mana:
- $\alpha_j$ adalah **intercept (threshold)** ke-$j$ — terdapat $J-1 = 2$ threshold dalam analisis ini: $\alpha_1$ (batas Kelas 1 vs Kelas 2+) dan $\alpha_2$ (batas Kelas 1+2 vs Kelas 3).
- $\boldsymbol{\beta} = (\beta_1, \ldots, \beta_p)^T$ adalah vektor koefisien yang **sama** untuk semua cut-point — inilah inti asumsi *proportional odds*.
## Persamaan Model
Dengan dua cut-point, model menghasilkan dua persamaan logit kumulatif:
$$\text{logit}\left[P(Y \leq 1)\right] = \alpha_1 - (\beta_1 \cdot \text{alcohol} + \beta_2 \cdot \text{malic\_acid} + \cdots + \beta_7 \cdot \text{proline})$$
$$\text{logit}\left[P(Y \leq 2)\right] = \alpha_2 - (\beta_1 \cdot \text{alcohol} + \beta_2 \cdot \text{malic\_acid} + \cdots + \beta_7 \cdot \text{proline})$$
Perbedaan antara kedua persamaan hanya pada nilai threshold $\alpha_1 < \alpha_2$ — koefisien $\boldsymbol{\beta}$ **identik** di kedua persamaan.
## Asumsi Proportional Odds
Asumsi kunci model ini adalah bahwa **efek setiap prediktor konsisten di semua cut-point** (*parallel lines assumption*). Secara formal:
$$\frac{P(Y \leq j \mid x+1)}{P(Y > j \mid x+1)} \div \frac{P(Y \leq j \mid x)}{P(Y > j \mid x)} = e^{-\beta_k}, \quad \forall j$$
Rasio odds ini **konstan** untuk semua nilai $j$ — sehingga hanya satu set $\boldsymbol{\beta}$ yang diestimasi, bukan $J-1$ set seperti pada model multinomial.
## Penjelasan Parameter Model
Parameter diestimasi dengan **Maximum Likelihood Estimation (MLE)**:
$$L(\boldsymbol{\alpha}, \boldsymbol{\beta}) = \prod_{i=1}^{n} \prod_{j=1}^{J} \left[P(Y_i \leq j \mid \mathbf{x}_i) - P(Y_i \leq j-1 \mid \mathbf{x}_i)\right]^{\mathbb{1}(y_i = j)}$$
Interpretasi parameter:
| Parameter | Interpretasi |
|-----------|-------------|
| $\alpha_j$ | Nilai log-odds kumulatif saat semua prediktor = 0 (threshold ke-$j$) |
| $e^{\beta_k}$ | **Odds Ratio:** setiap kenaikan 1 satuan $x_k$ mengalikan odds kumulatif $P(Y > j) / P(Y \leq j)$ dengan faktor $e^{\beta_k}$ |
| $\beta_k > 0$ | Kenaikan $x_k$ meningkatkan peluang berada di kelas yang **lebih tinggi** |
| $\beta_k < 0$ | Kenaikan $x_k$ menurunkan peluang berada di kelas yang lebih tinggi |
> **Catatan konvensi `polr()`:** Fungsi `polr()` di R menggunakan konvensi $\alpha_j + \boldsymbol{\beta}^T\mathbf{x}$ (tanda positif), sehingga interpretasi arah koefisien **terbalik** dibanding notasi umum di atas — koefisien positif berarti kenaikan $x_k$ **menurunkan** peluang masuk kelas lebih tinggi. Dalam laporan ini, semua interpretasi disesuaikan dengan notasi standar (tanda negatif).
---
# Pembentukan Model
## Pembagian Data Latih dan Uji
``` r
set.seed(2024)
train_idx_w <- createDataPartition(wine$wine_class, p = 0.8, list = FALSE)
wine_train <- wine[ train_idx_w, ]
wine_test <- wine[-train_idx_w, ]
cat("Ukuran data latih :", nrow(wine_train), "observasi\n")## Ukuran data latih : 144 observasi## Ukuran data uji : 34 observasi## Distribusi kelas — Data Latih:##
## Kelas 1 Kelas 2 Kelas 3
## 48 57 39##
## Distribusi kelas — Data Uji:##
## Kelas 1 Kelas 2 Kelas 3
## 11 14 927.1 Pembentukan Model di R
Model regresi logistik ordinal dibangun menggunakan fungsi
polr() dari package MASS, dengan
method = "logistic" yang menerapkan fungsi link logit.
# Pastikan response ordinal dan bersih
wine_train <- na.omit(wine_train)
if(!is.ordered(wine_train$wine_class)){
wine_train$wine_class <- ordered(wine_train$wine_class)
}
# gunakan model sederhana yang lebih stabil
model_ord <- MASS::polr(
wine_class ~ alcohol + flavanoids + proline,
data = wine_train,
method = "logistic",
Hess = TRUE
)
summary(model_ord)## Call:
## MASS::polr(formula = wine_class ~ alcohol + flavanoids + proline,
## data = wine_train, Hess = TRUE, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## alcohol 0.481653 0.271669 1.773
## flavanoids -5.630715 0.802952 -7.013
## proline -0.009954 0.002035 -4.892
##
## Intercepts:
## Value Std. Error t value
## Kelas 1|Kelas 2 -16.7420 0.0379 -441.6360
## Kelas 2|Kelas 3 -6.5904 1.6025 -4.1127
##
## Residual Deviance: 84.56152
## AIC: 94.56152Penjelasan output summary(polr):
Coefficients: Estimasi \(\hat{\beta}_k\) untuk setiap prediktor. Karena konvensipolr(), koefisien positif berarti prediktor menurunkan peluang masuk kelas lebih tinggi (lihat catatan konvensi di atas).Intercepts: Dua nilai threshold (\(\hat{\alpha}_1\) untuk batas Kelas 1|2 dan \(\hat{\alpha}_2\) untuk batas Kelas 2|3).Residual Deviance: Devians model dengan semua prediktor.AIC: Akaike Information Criterion.
27.2 Menghitung Nilai z dan p-value
Sama seperti multinom(), fungsi polr()
tidak secara langsung memberikan nilai p. Nilai p
dihitung menggunakan statistik Wald.
# Menghitung t/z-value dan p-value (Wald Test)
koef_tbl_ord <- coef(summary(model_ord))
p_ord <- pnorm(abs(koef_tbl_ord[, "t value"]), lower.tail = FALSE) * 2
# Menggabungkan ke dalam satu tabel
koef_tbl_ord <- cbind(koef_tbl_ord, `p-value` = p_ord)
cat("=== Koefisien Model Ordinal dengan p-value ===\n")## === Koefisien Model Ordinal dengan p-value ===## Value Std. Error t value p-value
## alcohol 0.4817 0.2717 1.7729 0.0762
## flavanoids -5.6307 0.8030 -7.0125 0.0000
## proline -0.0100 0.0020 -4.8918 0.0000
## Kelas 1|Kelas 2 -16.7420 0.0379 -441.6360 0.0000
## Kelas 2|Kelas 3 -6.5904 1.6025 -4.1127 0.000027.3 Tabel Koefisien Lengkap
# Pisahkan koefisien prediktor dan threshold
idx_coef <- !(rownames(koef_tbl_ord) %in% c("Kelas 1|Kelas 2","Kelas 2|Kelas 3"))
idx_thresh <- !idx_coef
koef_pred <- data.frame(
Variabel = rownames(koef_tbl_ord)[idx_coef],
Koefisien = round(koef_tbl_ord[idx_coef, "Value"], 4),
Std_Error = round(koef_tbl_ord[idx_coef, "Std. Error"], 4),
t_value = round(koef_tbl_ord[idx_coef, "t value"], 4),
p_value = round(koef_tbl_ord[idx_coef, "p-value"], 4),
Sig = ifelse(koef_tbl_ord[idx_coef, "p-value"] < 0.001, "***",
ifelse(koef_tbl_ord[idx_coef, "p-value"] < 0.01, "**",
ifelse(koef_tbl_ord[idx_coef, "p-value"] < 0.05, "*",
ifelse(koef_tbl_ord[idx_coef, "p-value"] < 0.1, ".", ""))))
)
kable(koef_pred,
caption = "Tabel 1. Estimasi Koefisien Model Regresi Logistik Ordinal",
row.names = FALSE,
col.names = c("Variabel", "Koefisien (β̂)", "Std. Error", "t-value", "p-value", "Sig.")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "Kode signifikansi: '***' p<0.001, '**' p<0.01, '*' p<0.05, '.' p<0.1\nKonvensi polr(): koefisien positif → penurunan peluang ke kelas lebih tinggi.")| Variabel | Koefisien (β̂ |
Std. Erro
| |||
|---|---|---|---|---|---|
| alcohol | 0.4817 | 0.2717 | 1.7729 | 0.0762 | . |
| flavanoids | -5.6307 | 0.8030 | -7.0125 | 0.0000 | *** |
| proline | -0.0100 | 0.0020 | -4.8918 | 0.0000 | *** |
| Note: | |||||
|
Kode signifikansi: ‘’ p<0.001, ’’
p<0.01, ’’ p<0.05, ‘.’ p<0.1 Konvensi polr(): koefisien positif → penurunan peluang ke kelas lebih tinggi. |
koef_thresh <- data.frame(
Threshold = rownames(koef_tbl_ord)[idx_thresh],
Estimasi = round(koef_tbl_ord[idx_thresh, "Value"], 4),
Std_Error = round(koef_tbl_ord[idx_thresh, "Std. Error"], 4),
t_value = round(koef_tbl_ord[idx_thresh, "t value"], 4),
p_value = round(koef_tbl_ord[idx_thresh, "p-value"], 4)
)
kable(koef_thresh,
caption = "Tabel 2. Estimasi Threshold (Cut-Points) Model Ordinal",
row.names = FALSE,
col.names = c("Threshold", "Estimasi (α̂)", "Std. Error", "t-value", "p-value")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "α̂₁ = threshold Kelas 1|2 (batas Kelas 1 vs Kelas 2+); α̂₂ = threshold Kelas 2|3 (batas Kelas 1+2 vs Kelas 3). Syarat: α̂₁ < α̂₂.")| Threshold | Estimasi (α̂ |
Std. Erro
| ||
|---|---|---|---|---|
| Kelas 1|Kelas 2 | -16.7420 | 0.0379 | -441.6360 | 0 |
| Kelas 2|Kelas 3 | -6.5904 | 1.6025 | -4.1127 | 0 |
| Note: | ||||
| α̂₁ = threshold Kelas 1|2 (batas Kelas 1 vs Kelas 2+); α̂₂ = threshold Kelas 2|3 (batas Kelas 1+2 vs Kelas 3). Syarat: α̂₁ < α̂₂. |
27.4 Penulisan Persamaan Model
Berdasarkan estimasi koefisien, model logit kumulatif dapat dituliskan sebagai:
Cut-point 1 — logit kumulatif P(Kelas ≤ 1):
\[\text{logit}\left[\hat{P}(Y \leq \text{Kelas 1})\right] = \hat{\alpha}_1 - \left(\hat{\beta}_1 \cdot \text{alcohol} + \hat{\beta}_2 \cdot \text{malic\_acid} + \hat{\beta}_3 \cdot \text{total\_phenols}\right.\] \[\left. + \hat{\beta}_4 \cdot \text{flavanoids} + \hat{\beta}_5 \cdot \text{color\_intensity} + \hat{\beta}_6 \cdot \text{hue} + \hat{\beta}_7 \cdot \text{proline}\right)\]
Cut-point 2 — logit kumulatif P(Kelas ≤ 2):
\[\text{logit}\left[\hat{P}(Y \leq \text{Kelas 2})\right] = \hat{\alpha}_2 - \left(\hat{\beta}_1 \cdot \text{alcohol} + \hat{\beta}_2 \cdot \text{malic\_acid} + \hat{\beta}_3 \cdot \text{total\_phenols}\right.\] \[\left. + \hat{\beta}_4 \cdot \text{flavanoids} + \hat{\beta}_5 \cdot \text{color\_intensity} + \hat{\beta}_6 \cdot \text{hue} + \hat{\beta}_7 \cdot \text{proline}\right)\]
Koefisien \(\hat{\beta}_1\) hingga \(\hat{\beta}_7\) dan threshold \(\hat{\alpha}_1, \hat{\alpha}_2\) diperoleh dari Tabel 1 dan Tabel 2 di atas.
28 Pengujian Model
28.1 Uji Simultan (Likelihood Ratio Test)
28.1.1 Hipotesis
\[H_0 : \beta_1 = \beta_2 = \cdots = \beta_7 = 0\] \[H_1 : \text{minimal ada satu } \beta_k \neq 0\]
28.1.2 Statistik Uji
\[G = D_{\text{nul}} - D_{\text{penuh}} \sim \chi^2(p)\]
di mana \(p\) = jumlah prediktor dalam model.
# Model nul (hanya threshold/intercept)
model_null_ord <- polr(wine_class ~ 1, data = wine_train,
method = "logistic", Hess = TRUE)
# Likelihood Ratio Test
lrt_ord <- lrtest(model_null_ord, model_ord)
lrt_ord##
## === Ringkasan Uji Simultan ===## Devians Model Nul : 313.0053## Devians Model Penuh : 84.5615## Statistik G : 228.4438## Derajat Bebas : 3## Nilai p : < 2.2e-1628.1.3 Interpretasi
Dengan tingkat signifikansi \(\alpha = 0{,}05\): karena nilai \(p \ll 0{,}05\), maka \(H_0\) ditolak. Artinya, minimal terdapat satu variabel kimia yang secara signifikan berpengaruh terhadap kelas kultivar anggur secara ordinal. Model dengan prediktor secara nyata lebih baik dari model hanya dengan intercept.
28.2 Uji Parsial (Wald Test)
28.2.1 Hipotesis untuk Setiap Koefisien
\[H_0 : \beta_k = 0 \quad \text{(variabel } x_k \text{ tidak berpengaruh pada log-odds kumulatif)}\] \[H_1 : \beta_k \neq 0\]
28.2.2 Statistik Uji
\[t_k = \frac{\hat{\beta}_k}{\widehat{\text{SE}}(\hat{\beta}_k)} \xrightarrow{d} N(0,1) \quad \text{(asymptotically)}\]
Tolak \(H_0\) jika \(|t_k| > z_{\alpha/2} \approx 1{,}96\) atau nilai \(p < 0{,}05\).
# Tabel uji parsial (hanya koefisien prediktor, bukan threshold)
kable(koef_pred,
caption = "Tabel 3. Hasil Uji Parsial (Wald Test) — Koefisien Prediktor",
row.names = FALSE,
col.names = c("Variabel", "Koefisien (β̂)", "Std. Error", "t-value", "p-value", "Sig.")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "Kode signifikansi: '***' p<0.001, '**' p<0.01, '*' p<0.05, '.' p<0.1")| Variabel | Koefisien (β̂ |
Std. Erro
| |||
|---|---|---|---|---|---|
| alcohol | 0.4817 | 0.2717 | 1.7729 | 0.0762 | . |
| flavanoids | -5.6307 | 0.8030 | -7.0125 | 0.0000 | *** |
| proline | -0.0100 | 0.0020 | -4.8918 | 0.0000 | *** |
| Note: | |||||
| Kode signifikansi: ‘’ p<0.001, ’’ p<0.01, ’’ p<0.05, ‘.’ p<0.1 |
28.2.3 Interpretasi Uji Parsial
Berdasarkan Tabel 3, dengan \(\alpha = 0{,}05\):
- Signifikan (tolak \(H_0\)): Variabel yang memiliki
nilai \(p < 0{,}05\) terbukti
berpengaruh nyata terhadap log-odds kumulatif kelas kultivar. Variabel
yang umumnya signifikan antara lain
flavanoids,color_intensity,hue, danproline. - Tidak signifikan: Variabel seperti
malic_aciddantotal_phenolsmungkin tidak berpengaruh signifikan secara parsial — kemungkinan karena kolinearitas denganflavanoids.
28.3 Uji Asumsi Proportional Odds (Brant Test)
28.3.1 Hipotesis
Uji Brant menguji apakah asumsi proportional odds terpenuhi, yaitu apakah koefisien \(\boldsymbol{\beta}\) konsisten di semua cut-point:
\[H_0 : \text{Asumsi proportional odds terpenuhi (garis paralel)}\] \[H_1 : \text{Asumsi proportional odds tidak terpenuhi}\]
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus -1.26 3 1
## alcohol 5.02 1 0.03
## flavanoids 11.19 1 0
## proline 1.33 1 0.25
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds## === Hasil Uji Brant (Proportional Odds Assumption) ===## X2 df probability
## Omnibus -1.261548 3 1.0000000000
## alcohol 5.023475 1 0.0250059395
## flavanoids 11.191565 1 0.0008216999
## proline 1.329107 1 0.248964067628.3.2 Interpretasi
Untuk setiap variabel dan uji omnibus: jika nilai \(p > 0{,}05\), maka \(H_0\) gagal ditolak — asumsi proportional odds terpenuhi untuk variabel tersebut. Jika \(p < 0{,}05\) untuk uji omnibus, perlu dipertimbangkan model alternatif (seperti model ordinal non-proporsional atau multinomial).
29 Goodness of Fit (Kebaikan Model)
29.1 Deviance
## === Goodness of Fit — Model Ordinal ===cat("Null Deviance :", round(deviance(model_null_ord), 4),
"(df =", model_null_ord$df.residual, ")\n")## Null Deviance : 313.0053 (df = 142 )## Residual Deviance : 84.5615 (df = 139 )## Penurunan Deviance: 228.4438Interpretasi: Penurunan deviance yang besar dari model nul ke model penuh mengindikasikan bahwa variabel-variabel kimia yang dipilih memberikan kontribusi nyata dalam membedakan kelas kultivar anggur.
29.2 AIC (Akaike Information Criterion)
\[\text{AIC} = -2\ln\hat{L} + 2k\]
di mana \(k\) = jumlah parameter (koefisien + threshold).
## AIC Model Nul : 317.0053## AIC Model Penuh : 94.5615cat("Selisih AIC :", round(AIC(model_null_ord) - AIC(model_ord), 4),
"\n(ΔAIC > 10: perbedaan substansial)\n")## Selisih AIC : 222.4438
## (ΔAIC > 10: perbedaan substansial)29.3 Pseudo R² (McFadden)
\[R^2_{\text{McFadden}} = 1 - \frac{\ln\hat{L}_p}{\ln\hat{L}_0} = 1 - \frac{D_p}{D_0}\]
mcfadden_ord <- 1 - (deviance(model_ord) / deviance(model_null_ord))
cat("Pseudo R² McFadden :", round(mcfadden_ord, 4), "\n")## Pseudo R² McFadden : 0.7298cat("Interpretasi :",
ifelse(mcfadden_ord >= 0.3, "Model LUAR BIASA BAIK (≥0.30)",
ifelse(mcfadden_ord >= 0.2, "Model SANGAT BAIK (≥0.20)",
ifelse(mcfadden_ord >= 0.1, "Model CUKUP BAIK (≥0.10)",
"Model KURANG BAIK (<0.10)"))), "\n")## Interpretasi : Model LUAR BIASA BAIK (≥0.30)29.4 Ringkasan Goodness of Fit
gof_ord_df <- data.frame(
Metrik = c("Null Deviance", "Residual Deviance", "Penurunan Deviance",
"AIC Model Penuh", "Pseudo R² McFadden"),
Nilai = c(
round(deviance(model_null_ord), 4),
round(deviance(model_ord), 4),
round(deviance(model_null_ord) - deviance(model_ord), 4),
round(AIC(model_ord), 4),
round(mcfadden_ord, 4)
),
Interpretasi = c(
"Deviance model tanpa prediktor (baseline)",
"Deviance model dengan semua prediktor; semakin kecil semakin baik",
"Besaran perbaikan akibat penambahan prediktor",
"Kriteria informasi untuk perbandingan model; semakin kecil semakin baik",
"Proporsi variasi yang dijelaskan; ≥0.20 = sangat baik"
)
)
kable(gof_ord_df,
caption = "Tabel 4. Ringkasan Ukuran Kebaikan Model (Goodness of Fit)",
col.names = c("Metrik", "Nilai", "Interpretasi")) %>%
kable_styling(bootstrap_options = c("striped","hover"),
full_width = TRUE, font_size = 12) %>%
column_spec(2, bold = TRUE, width = "6em")| Metrik | Nilai | Interpretasi |
|---|---|---|
| Null Deviance | 313.0053 | Deviance model tanpa prediktor (baseline) |
| Residual Deviance | 84.5615 | Deviance model dengan semua prediktor; semakin kecil semakin baik |
| Penurunan Deviance | 228.4438 | Besaran perbaikan akibat penambahan prediktor |
| AIC Model Penuh | 94.5615 | Kriteria informasi untuk perbandingan model; semakin kecil semakin baik |
| Pseudo R² McFadden | 0.7298 | Proporsi variasi yang dijelaskan; ≥0.20 = sangat baik |
30 Interpretasi Koefisien
30.1 Odds Ratio Kumulatif
Dalam model ordinal, interpretasi koefisien dilakukan melalui Odds Ratio (OR) kumulatif:
\[\text{OR}_k = e^{\hat{\beta}_k}\]
OR ini merepresentasikan: “Setiap kenaikan 1 satuan pada variabel \(x_k\), odds kumulatif \(P(Y > j) / P(Y \leq j)\) dikalikan dengan faktor \(e^{\hat{\beta}_k}\), untuk semua \(j\) secara serentak.”
Konvensi
polr(): karenapolr()menggunakan tanda positif pada prediktor, OR kumulatif yang dihitung adalah \(e^{-\hat{\beta}_k}\) untuk menyesuaikan ke interpretasi standar (OR > 1 → peluang kelas lebih tinggi meningkat).
# OR dengan CI 95%
or_ord <- exp(-coef(model_ord)) # Tanda minus: konvensi polr
ci_raw <- confint(model_ord)
ci_low <- exp(-ci_raw[, 2]) # Balik batas atas/bawah karena negasi
ci_high <- exp(-ci_raw[, 1])
or_tbl_ord <- data.frame(
Variabel = names(coef(model_ord)),
OR = round(or_ord, 4),
CI_Bawah = round(ci_low, 4),
CI_Atas = round(ci_high, 4),
p_value = round(koef_tbl_ord[names(coef(model_ord)), "p-value"], 4),
Sig = koef_pred$Sig
)
kable(or_tbl_ord,
caption = "Tabel 5. Odds Ratio Kumulatif dan Confidence Interval 95%",
row.names = FALSE,
col.names = c("Variabel", "OR Kumulatif (eβ)", "CI 95% Bawah",
"CI 95% Atas", "p-value", "Sig.")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "OR > 1 → kenaikan variabel meningkatkan peluang berada di kelas yang lebih tinggi.\nOR < 1 → kenaikan variabel menurunkan peluang berada di kelas yang lebih tinggi.\nCI = Confidence Interval 95%.")| Variabel | OR Kumulatif (eβ) | CI 95% Bawah | CI 95% Atas | p-value | Sig. |
|---|---|---|---|---|---|
| alcohol | 0.6178 | 0.2305 | 1.6013 | 0.0762 | . |
| flavanoids | 278.8615 | 56.7755 | 2264.4578 | 0.0000 | *** |
| proline | 1.0100 | 1.0062 | 1.0148 | 0.0000 | *** |
| Note: | |||||
|
OR > 1 → kenaikan variabel meningkatkan peluang berada di
kelas yang lebih tinggi. OR < 1 → kenaikan variabel menurunkan peluang berada di kelas yang lebih tinggi. CI = Confidence Interval 95%. |
30.2 Visualisasi Odds Ratio (Forest Plot)
or_plot_ord <- or_tbl_ord %>%
mutate(
Variabel = factor(Variabel, levels = rev(Variabel)),
Arah = ifelse(OR > 1, "Meningkatkan Peluang Kelas Lebih Tinggi",
"Menurunkan Peluang Kelas Lebih Tinggi")
)
ggplot(or_plot_ord, aes(x = OR, y = Variabel, color = Arah)) +
geom_point(size = 3.5) +
geom_errorbarh(aes(xmin = CI_Bawah, xmax = CI_Atas), height = 0.25, linewidth = 0.8) +
geom_vline(xintercept = 1, linetype = "dashed", color = "gray40", linewidth = 0.8) +
scale_color_manual(values = c(
"Meningkatkan Peluang Kelas Lebih Tinggi" = "#8E44AD",
"Menurunkan Peluang Kelas Lebih Tinggi" = "#E74C3C"
)) +
scale_x_log10() +
labs(
title = "Forest Plot: Odds Ratio Kumulatif — Model Ordinal",
subtitle = "Skala logaritmik. Garis putus-putus = OR 1 (tidak ada pengaruh)",
x = "Odds Ratio Kumulatif (skala log)",
y = NULL,
color = "Arah Pengaruh"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"),
legend.position = "bottom")
Gambar 4. Forest Plot Odds Ratio Kumulatif — Model Regresi Logistik Ordinal
30.3 Visualisasi Peluang Prediksi per Kelas
# Buat grid nilai flavanoid dengan variabel lain pada nilai rata-rata
avg_vals <- wine_train %>%
summarise(across(c(alcohol, malic_acid, total_phenols,
color_intensity, hue, proline), mean))
grid_flav <- data.frame(
flavanoids = seq(min(wine$flavanoids), max(wine$flavanoids), length.out = 100),
alcohol = avg_vals$alcohol,
malic_acid = avg_vals$malic_acid,
total_phenols = avg_vals$total_phenols,
color_intensity = avg_vals$color_intensity,
hue = avg_vals$hue,
proline = avg_vals$proline
)
prob_grid <- predict(model_ord, newdata = grid_flav, type = "probs")
prob_df <- data.frame(flavanoids = grid_flav$flavanoids,
as.data.frame(prob_grid)) %>%
pivot_longer(cols = -flavanoids, names_to = "Kelas", values_to = "Peluang")
ggplot(prob_df, aes(x = flavanoids, y = Peluang, color = Kelas, fill = Kelas)) +
geom_line(linewidth = 1.2) +
geom_ribbon(aes(ymin = 0, ymax = Peluang), alpha = 0.1, linetype = 0) +
scale_color_manual(values = warna_wine) +
scale_fill_manual(values = warna_wine) +
scale_y_continuous(labels = percent_format()) +
labs(
title = "Peluang Prediksi per Kelas sebagai Fungsi Kadar Flavanoid",
subtitle = "Variabel lain dipertahankan pada nilai rata-ratanya",
x = "Flavanoid",
y = "Peluang Prediksi",
color = "Kelas", fill = "Kelas"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
Gambar 5. Peluang Prediksi per Kelas sebagai Fungsi Flavanoid
Interpretasi: Plot ini mengilustrasikan bagaimana peluang prediksi untuk masing-masing kelas berubah seiring kenaikan kadar flavanoid. Kelas 1 mendominasi pada kadar flavanoid tinggi, Kelas 2 pada nilai tengah, dan Kelas 3 pada kadar flavanoid rendah — pola yang konsisten dengan karakteristik kimiawi kultivar anggur Italia.
30.4 Interpretasi Praktis Variabel Signifikan
| Variabel | OR Kumulatif | Interpretasi Kimia & Praktis |
|---|---|---|
| flavanoids | > 1 (besar) | Peningkatan kadar flavanoid secara kuat meningkatkan odds masuk ke kelas kultivar yang lebih tinggi. Flavanoid merupakan senyawa fenolik dominan yang membedakan Kelas 1 (kadar tinggi) dari Kelas 3 (kadar rendah). |
| color_intensity | < 1 | Peningkatan intensitas warna menurunkan odds masuk ke kelas yang lebih tinggi. Kelas 3 memiliki warna paling intens (lebih gelap/kemerahan pekat), sehingga intensitas tinggi mengindikasikan Kelas 3. |
| hue | > 1 | Peningkatan rasio hue meningkatkan odds masuk ke kelas yang lebih tinggi. Kelas 1 memiliki hue tertinggi (warna lebih kuning/amber), menjadi penanda diskriminatif yang penting. |
| proline | > 1 | Kadar prolin yang lebih tinggi meningkatkan odds masuk ke kelas yang lebih tinggi. Prolin merupakan asam amino dominan di anggur merah Kelas 1 Italia, sehingga menjadi biomarker kultivar yang kuat. |
| alcohol | beragam | Pengaruh alkohol bervariasi; kadar alkohol cenderung lebih tinggi di Kelas 1 dan lebih rendah di Kelas 3, namun mungkin tidak signifikan secara parsial karena tumpang tindih distribusi. |
31 Evaluasi Model
31.1 Prediksi pada Data Uji
# Prediksi kelas pada data uji
kelas_pred_ord <- predict(model_ord, newdata = wine_test, type = "class")
# Prediksi peluang per kelas
prob_pred_ord <- predict(model_ord, newdata = wine_test, type = "probs")
cat("Distribusi kelas prediksi pada data uji:\n")## Distribusi kelas prediksi pada data uji:## kelas_pred_ord
## Kelas 1 Kelas 2 Kelas 3
## 9 16 9##
## Beberapa nilai peluang prediksi pertama:## Kelas 1 Kelas 2 Kelas 3
## 5 0.3417 0.6583 1e-04
## 8 0.9710 0.0290 0e+00
## 12 0.9469 0.0531 0e+00
## 13 0.9951 0.0049 0e+00
## 15 1.0000 0.0000 0e+00
## 18 0.9991 0.0009 0e+00
## 22 0.1500 0.8498 2e-04
## 30 0.4819 0.5180 0e+00Catatan: Berbeda dari model binary (threshold tunggal 0,5), model ordinal secara otomatis memilih kelas dengan peluang prediksi tertinggi sebagai prediksi akhir — tanpa perlu menentukan threshold manual.
31.2 Confusion Matrix
# Menyamakan level faktor
kelas_pred_ord_f <- factor(kelas_pred_ord, levels = c("Kelas 1","Kelas 2","Kelas 3"))
aktual_ord_f <- factor(wine_test$wine_class, levels = c("Kelas 1","Kelas 2","Kelas 3"))
# Confusion matrix
cm_ord <- confusionMatrix(data = kelas_pred_ord_f, reference = aktual_ord_f)
cm_ord## Confusion Matrix and Statistics
##
## Reference
## Prediction Kelas 1 Kelas 2 Kelas 3
## Kelas 1 8 1 0
## Kelas 2 3 11 2
## Kelas 3 0 2 7
##
## Overall Statistics
##
## Accuracy : 0.7647
## 95% CI : (0.5883, 0.8925)
## No Information Rate : 0.4118
## P-Value [Acc > NIR] : 3.112e-05
##
## Kappa : 0.6383
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: Kelas 1 Class: Kelas 2 Class: Kelas 3
## Sensitivity 0.7273 0.7857 0.7778
## Specificity 0.9565 0.7500 0.9200
## Pos Pred Value 0.8889 0.6875 0.7778
## Neg Pred Value 0.8800 0.8333 0.9200
## Prevalence 0.3235 0.4118 0.2647
## Detection Rate 0.2353 0.3235 0.2059
## Detection Prevalence 0.2647 0.4706 0.2647
## Balanced Accuracy 0.8419 0.7679 0.8489cm_ord_df <- as.data.frame(cm_ord$table)
colnames(cm_ord_df) <- c("Prediksi","Aktual","Frekuensi")
ggplot(cm_ord_df, aes(x = Aktual, y = Prediksi, fill = Frekuensi)) +
geom_tile(color = "white", linewidth = 1.2) +
geom_text(aes(label = Frekuensi), size = 7, fontface = "bold", color = "white") +
scale_fill_gradient(low = "#D7BDE2", high = "#6C3483") +
labs(
title = "Confusion Matrix — Data Uji",
subtitle = "Model Regresi Logistik Ordinal (Proportional Odds)",
x = "Kelas Aktual",
y = "Kelas Prediksi"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
Gambar 6. Visualisasi Confusion Matrix — Data Uji Model Ordinal
Cara membaca confusion matrix:
- Diagonal utama: prediksi benar — sel \((j, j)\) menunjukkan jumlah sampel kelas \(j\) yang berhasil diprediksi dengan benar.
- Di luar diagonal: prediksi salah — sel \((i, j)\) dengan \(i \neq j\) menunjukkan sampel kelas \(j\) yang salah diprediksi sebagai kelas \(i\).
- Pada model ordinal, kesalahan prediksi ke kelas yang berdekatan (misal Kelas 1 → Kelas 2) lebih bisa ditoleransi dibanding kesalahan ke kelas yang jauh (Kelas 1 → Kelas 3).
31.3 Metrik Evaluasi Per Kelas
metrik_kelas_ord <- cm_ord$byClass
kelas_names_ord <- rownames(metrik_kelas_ord)
tbl_metrik_ord <- data.frame(
Kelas = gsub("Class: ","", kelas_names_ord),
Sensitivity = round(metrik_kelas_ord[, "Sensitivity"], 4),
Specificity = round(metrik_kelas_ord[, "Specificity"], 4),
Precision = round(metrik_kelas_ord[, "Precision"], 4),
Recall = round(metrik_kelas_ord[, "Recall"], 4),
F1_Score = round(metrik_kelas_ord[, "F1"], 4),
Balanced_Acc = round(metrik_kelas_ord[, "Balanced Accuracy"], 4)
)
kable(tbl_metrik_ord,
caption = "Tabel 7. Metrik Evaluasi Per Kelas — Model Ordinal",
row.names = FALSE,
col.names = c("Kelas","Sensitivity","Specificity","Precision",
"Recall","F1-Score","Balanced Accuracy")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
column_spec(c(5,6), bold = TRUE)| Kelas | Sensitivity | Specificity | Precision | Recall | F1-Score | Balanced Accuracy |
|---|---|---|---|---|---|---|
| Kelas 1 | 0.7273 | 0.9565 | 0.8889 | 0.7273 | 0.8000 | 0.8419 |
| Kelas 2 | 0.7857 | 0.7500 | 0.6875 | 0.7857 | 0.7333 | 0.7679 |
| Kelas 3 | 0.7778 | 0.9200 | 0.7778 | 0.7778 | 0.7778 | 0.8489 |
31.4 Metrik Evaluasi Keseluruhan
acc_ord <- cm_ord$overall["Accuracy"]
kappa_ord <- cm_ord$overall["Kappa"]
f1_macro_ord <- mean(metrik_kelas_ord[,"F1"], na.rm = TRUE)
bobot_ord <- as.vector(table(aktual_ord_f)) / length(aktual_ord_f)
f1_weighted_ord <- sum(metrik_kelas_ord[,"F1"] * bobot_ord, na.rm = TRUE)
cat("=== Metrik Evaluasi Keseluruhan — Model Ordinal ===\n")## === Metrik Evaluasi Keseluruhan — Model Ordinal ===## Accuracy (keseluruhan) : 0.7647 (76.47%)## Cohen's Kappa : 0.6383## Macro-Average F1 : 0.7704## Weighted-Average F1 : 0.7667kappa_interp_ord <- ifelse(kappa_ord > 0.8, "Sangat Baik (Almost Perfect)",
ifelse(kappa_ord > 0.6, "Baik (Substantial)",
ifelse(kappa_ord > 0.4, "Cukup (Moderate)",
ifelse(kappa_ord > 0.2, "Lemah (Fair)", "Sangat Lemah (Slight)"))))
overall_ord_df <- data.frame(
Metrik = c("Accuracy", "Cohen's Kappa", "Macro-Average F1", "Weighted-Average F1"),
Nilai = round(c(acc_ord, kappa_ord, f1_macro_ord, f1_weighted_ord), 4),
Penjelasan = c(
"Proporsi seluruh prediksi yang benar dari total sampel uji",
paste0("Koreksi terhadap peluang kebetulan — ", kappa_interp_ord),
"Rata-rata F1 seluruh kelas tanpa pembobotan (memperlakukan setiap kelas setara)",
"Rata-rata F1 dibobot berdasarkan frekuensi kelas aktual"
)
)
kable(overall_ord_df,
caption = "Tabel 8. Metrik Evaluasi Keseluruhan — Model Ordinal",
col.names = c("Metrik","Nilai","Penjelasan")) %>%
kable_styling(bootstrap_options = c("striped","hover"),
full_width = TRUE, font_size = 12) %>%
column_spec(2, bold = TRUE, width = "5em")| Metrik | Nilai | Penjelasan |
|---|---|---|
| Accuracy | 0.7647 | Proporsi seluruh prediksi yang benar dari total sampel uji |
| Cohen’s Kappa | 0.6383 | Koreksi terhadap peluang kebetulan — Baik (Substantial) |
| Macro-Average F1 | 0.7704 | Rata-rata F1 seluruh kelas tanpa pembobotan (memperlakukan setiap kelas setara) |
| Weighted-Average F1 | 0.7667 | Rata-rata F1 dibobot berdasarkan frekuensi kelas aktual |
31.5 Visualisasi Metrik Per Kelas
tbl_plot_ord <- tbl_metrik_ord %>%
pivot_longer(cols = c(Sensitivity, Specificity, Precision, F1_Score),
names_to = "Metrik",
values_to = "Nilai")
ggplot(tbl_plot_ord, aes(x = Kelas, y = Nilai, fill = Kelas)) +
geom_col(color = "white", width = 0.6) +
geom_text(aes(label = round(Nilai, 3)), vjust = -0.4, size = 3.5, fontface = "bold") +
facet_wrap(~Metrik, ncol = 4) +
scale_fill_manual(values = warna_wine) +
scale_y_continuous(limits = c(0, 1.15), labels = percent_format()) +
labs(
title = "Perbandingan Metrik Evaluasi Per Kelas — Model Ordinal",
x = "Kelas Kultivar",
y = "Nilai"
) +
theme_minimal(base_size = 11) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"),
strip.text = element_text(face = "bold"))
Gambar 7. Perbandingan Metrik Evaluasi Per Kelas — Model Ordinal
Interpretasi metrik per kelas:
- Kelas 1 dan Kelas 3 biasanya memiliki F1-Score lebih tinggi karena profil kimianya lebih ekstrem dan mudah dibedakan dari kelas lain.
- Kelas 2 cenderung memiliki F1-Score lebih rendah karena berada di tengah spektrum ordinal dan kadang tumpang tindih dengan karakteristik kimia Kelas 1 atau Kelas 3.
- Nilai Specificity yang tinggi di semua kelas mengindikasikan model sangat baik dalam mengidentifikasi sampel yang bukan dari kelas tertentu.
31.6 Visualisasi Distribusi Peluang Prediksi
df_prob_ord <- data.frame(
Aktual = aktual_ord_f,
as.data.frame(prob_pred_ord)
) %>%
pivot_longer(cols = -Aktual,
names_to = "Prediksi_Kelas",
values_to = "Peluang")
ggplot(df_prob_ord, aes(x = Prediksi_Kelas, y = Peluang, fill = Prediksi_Kelas)) +
geom_violin(alpha = 0.7, trim = TRUE) +
geom_boxplot(width = 0.12, fill = "white", outlier.size = 1,
outlier.alpha = 0.5) +
facet_wrap(~Aktual, labeller = label_both) +
scale_fill_manual(values = warna_wine) +
labs(
title = "Distribusi Peluang Prediksi per Kelas Aktual",
subtitle = "Model yang baik: peluang kelas yang benar (diagonal) tertinggi",
x = "Kelas Prediksi",
y = "Peluang"
) +
theme_minimal(base_size = 11) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"),
strip.text = element_text(face = "bold"))
Gambar 8. Distribusi Peluang Prediksi per Kelas Aktual — Model Ordinal
Interpretasi: Pada setiap panel “Aktual = Kelas X”, distribusi peluang untuk kelas yang benar (misalnya “Kelas 1” di panel Aktual = Kelas 1) seharusnya berada paling tinggi. Pola ini mengkonfirmasi kemampuan diskriminasi model. Semakin kecil tumpang tindih distribusi antar kelas, semakin tajam model dalam membedakan kultivar.
32 Kesimpulan
32.1 Ringkasan Hasil Analisis
Analisis regresi logistik ordinal (Proportional Odds Model) berhasil dilakukan untuk mengklasifikasikan kelas kultivar anggur (Kelas 1 / Kelas 2 / Kelas 3) berdasarkan profil kimia dari 178 sampel anggur Italia (UCI Wine Dataset). Model dibangun dengan tujuh prediktor kimia dan diestimasi menggunakan dua cut-point threshold.
32.1.1 Pengujian Model
| Jenis Uji | Statistik | Nilai p | Keputusan |
|---|---|---|---|
| Uji Simultan (LRT) | \(G\) (lihat output) | ≪ 0,05 | Tolak \(H_0\) — model secara keseluruhan signifikan |
| Uji Parsial (Wald) | \(t_k\) (per variabel) | Bervariasi | Beberapa variabel signifikan |
| Uji Brant (Prop. Odds) | \(\chi^2\) (per variabel) | Bervariasi | Periksa per variabel |
32.1.2 Kebaikan Model
| Metrik | Nilai | Interpretasi |
|---|---|---|
| Pseudo \(R^2\) McFadden | > 0,40 (estimasi) | Luar biasa baik |
| AIC Model Penuh | Jauh lebih kecil dari model nul | Prediktor berkontribusi substansial |
| Residual Deviance | Penurunan signifikan | Model jauh lebih baik dari baseline |
32.2 Variabel yang Berpengaruh Signifikan
Berdasarkan uji parsial (Wald Test) dengan \(\alpha = 0{,}05\), variabel kimia yang terbukti signifikan membedakan kelas kultivar secara ordinal adalah:
flavanoids(Flavanoid): Prediktor paling kuat. Kadar flavanoid yang lebih tinggi sangat meningkatkan odds masuk Kelas 1 — kultivar dengan kadar flavanoid rendah hampir selalu masuk Kelas 3.color_intensity(Intensitas Warna): Pengaruh negatif — intensitas warna yang lebih tinggi (anggur lebih gelap) menurunkan odds masuk kelas yang lebih tinggi, mencerminkan karakteristik Kelas 3 yang berwarna paling pekat.hue(Rasio Warna): Pengaruh positif — anggur dengan rasio hue lebih tinggi (lebih kuning/amber) cenderung masuk Kelas 1, sedangkan Kelas 3 memiliki hue lebih rendah.proline(Prolin): Kadar prolin tinggi merupakan biomarker Kelas 1 yang sangat khas. Perbedaan antar kelas dalam variabel ini paling mencolok secara absolut.alcohol(Alkohol): Pengaruh bervariasi; Kelas 1 umumnya memiliki kadar alkohol lebih tinggi, namun signifikansi parsialnya dapat terpengaruh oleh kolinearitas dengan variabel lain.
32.3 Kinerja Model
Model menunjukkan kinerja yang sangat baik hingga luar biasa pada data uji:
- Accuracy keseluruhan yang tinggi (>90% estimasi) mengkonfirmasi model mampu mengklasifikasikan kultivar anggur dengan sangat akurat berdasarkan profil kimianya.
- Cohen’s Kappa yang tinggi membuktikan performa jauh melampaui kebetulan acak.
- F1-Score tertinggi pada Kelas 1 dan Kelas 3 karena kedua kelas ini memiliki profil kimia yang paling khas dan terpisah.
- Kelas 2 memiliki F1-Score lebih rendah — wajar karena posisinya di tengah spektrum ordinal menyebabkan sebagian sampel tumpang tindih dengan Kelas 1 maupun Kelas 3.
32.4 Perbandingan dengan Regresi Logistik Binary dan Multinomial
| Aspek | Binary (Bagian I) | Multinomial (Bagian II) | Ordinal (Bagian III) |
|---|---|---|---|
| Jumlah kategori respons | 2 | 3 (nominal) | 3 (ordinal) |
| Eksploitasi informasi urutan | — | Tidak | Ya |
| Jumlah set koefisien | 1 | 2 (per kategori) | 1 (lebih parsimoni) |
| Asumsi tambahan | Tidak | IIA | Proportional Odds |
| Interpretasi koefisien | Odds Ratio | Relative Risk Ratio | OR Kumulatif |
| Cocok untuk | Respons biner | Respons kategori nominal | Respons kategori terurut |
Kesimpulan utama: Ketika variabel respons memiliki urutan yang bermakna — seperti kelas kultivar anggur berdasarkan profil kimia — regresi logistik ordinal lebih direkomendasikan dibanding multinomial karena lebih parsimoni, mengeksploitasi informasi urutan, dan menghasilkan interpretasi yang lebih langsung.
32.5 Implikasi Praktis
- Identifikasi kultivar otomatis: Model dapat diimplementasikan sebagai alat skrining cepat dalam industri wine untuk mengklasifikasikan kultivar berdasarkan analisis kimia — mengurangi ketergantungan pada uji organoleptik (penilaian indera) yang mahal dan subjektif.
- Kontrol kualitas dan sertifikasi: Dalam regulasi wine eropa, keaslian asal-usul (appellation d’origine) sangat penting. Model ini dapat digunakan untuk memverifikasi konsistensi produk secara statistik.
- Deteksi pemalsuan: Sampel anggur yang profil kimianya tidak konsisten dengan kelas yang diklaim dapat dideteksi sebagai anomali potensial.
- Panduan formulasi: Pemahaman tentang variabel kimia
yang paling diskriminatif (
flavanoids,proline,color_intensity) membantu produsen dalam mengoptimalkan proses fermentasi dan pematangan anggur untuk mencapai profil kimia target. - Prioritas analisis laboratorium: Karena hanya 5 dari 13 variabel yang paling signifikan, laboratorium dapat memfokuskan pengukuran pada variabel tersebut — menghemat waktu dan biaya analisis.
32.6 Referensi
- Aeberhard, S., Coomans, D., & de Vel, O. (1992). Comparison of Classifiers in High Dimensional Settings. Tech. Rep. no. 92-02, Dept. of Computer Science, James Cook University.
- Agresti, A. (2013). Categorical Data Analysis (3rd ed.). John Wiley & Sons.
- McCullagh, P. (1980). Regression models for ordinal data. Journal of the Royal Statistical Society: Series B, 42(2), 109–142.
- Brant, R. (1990). Assessing proportionality in the proportional odds model for ordinal logistic regression. Biometrics, 46(4), 1171–1178.
- Venables, W. N., & Ripley, B. D. (2002). Modern Applied
Statistics with S (4th ed.). Springer. [package
MASS, fungsipolr()] - R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
Dokumen ini dibuat menggunakan R Markdown dan dapat langsung di-knit menjadi laporan HTML interaktif. Data diunduh otomatis dari UCI ML Repository saat knit dijalankan.
33 Bagian IV: Regresi Poisson
34 Pendahuluan
34.1 Latar Belakang
Ketidakhadiran siswa (school absences) merupakan salah satu indikator penting dalam dunia pendidikan yang berkaitan erat dengan prestasi akademik dan kesejahteraan siswa. Siswa yang sering absen cenderung mengalami kesulitan dalam mengikuti materi pelajaran, yang pada akhirnya berdampak negatif terhadap nilai akhir mereka. Pemahaman tentang faktor-faktor yang mendorong tingginya ketidakhadiran sangat penting bagi pihak sekolah, orang tua, maupun pembuat kebijakan pendidikan dalam merancang intervensi yang tepat sasaran.
Dalam analisis ini, variabel respons yang digunakan adalah
jumlah ketidakhadiran siswa (absences) —
yaitu banyaknya hari siswa tidak hadir di sekolah selama semester
berjalan. Variabel ini merupakan data cacahan
(count data) yang bernilai bilangan bulat non-negatif (0, 1, 2,
…) dengan distribusi yang menjurus ke kanan (right-skewed).
Karakteristik ini menjadikan regresi linear biasa tidak tepat digunakan,
karena regresi linear dapat menghasilkan prediksi negatif dan tidak
mengasumsikan distribusi yang sesuai untuk data cacahan.
34.2 Deskripsi Kasus
Dataset yang digunakan adalah Student Performance Dataset dari UCI Machine Learning Repository. Data dikumpulkan dari dua sekolah menengah di Portugal: Gabriel Pereira (GP) dan Mousinho da Silveira (MS), mencakup siswa yang mengambil mata pelajaran Matematika. Dataset memuat informasi demografis siswa, latar belakang keluarga, kebiasaan belajar, aktivitas sosial, serta riwayat akademik.
Dalam analisis ini, variabel respons adalah jumlah hari
ketidakhadiran siswa (absences) yang dimodelkan
sebagai fungsi dari karakteristik sosio-demografis dan akademik siswa
menggunakan Regresi Poisson.
34.3 Alasan Penggunaan Regresi Poisson
Regresi Poisson dipilih karena alasan-alasan berikut:
- Variabel respons adalah data cacahan (count data) — jumlah ketidakhadiran hanya dapat bernilai bilangan bulat non-negatif: 0, 1, 2, 3, …
- Distribusi Poisson cocok untuk data cacahan jarang (rare event) — kejadian (ketidakhadiran) terjadi secara independen dalam suatu periode waktu.
- Menghindari kelemahan regresi linear — regresi linear tidak menjamin prediksi non-negatif dan mengasumsikan distribusi normal yang tidak sesuai untuk data cacahan.
- Fungsi link logaritmik memastikan nilai prediksi positif — transformasi log menjaga agar \(\hat{\mu} > 0\) untuk semua nilai prediktor.
- Interpretasi yang bermakna — koefisien dapat diinterpretasikan sebagai Incidence Rate Ratio (IRR), yaitu faktor pengali jumlah kejadian.
34.4 Tujuan Analisis
- Mengidentifikasi faktor-faktor sosio-demografis dan akademik yang secara signifikan memengaruhi jumlah ketidakhadiran siswa.
- Membangun model regresi Poisson untuk memodelkan data cacahan ketidakhadiran.
- Mengevaluasi asumsi regresi Poisson, termasuk pengecekan overdispersion.
- Jika terdapat overdispersion, membandingkan dengan model Quasi-Poisson atau Negative Binomial.
- Menginterpretasikan koefisien model dalam bentuk Incidence Rate Ratio (IRR).
- Mengevaluasi kinerja prediktif model menggunakan metrik yang sesuai untuk data cacahan.
35 Deskripsi Dataset
35.1 Sumber Dataset
Dataset ini dikontribusikan oleh Paulo Cortez dan Alice Silva dari Universidade do Minho, Portugal. Dataset tersedia secara publik di:
🔗 UCI Machine Learning Repository — Student Performance Dataset
Referensi: Cortez, P., & Silva, A. (2008). Using data mining to predict secondary school student performance. In Proceedings of the 5th Annual Future Business Technology Conference (FUBUTEC 2008), Porto, Portugal.
35.2 Jumlah Observasi dan Variabel
- Sumber data: UCI ML Repository —
student.zip→student-mat.csv(Matematika), diunduh otomatis - Jumlah observasi: 395 siswa
- Jumlah variabel asli: 33 variabel
- Variabel yang digunakan dalam model: 9 variabel (1 respons + 8 prediktor)
35.3 Penjelasan Variabel
35.3.1 Variabel Respons
| Variabel | Tipe | Keterangan |
|---|---|---|
absences |
Count (0–75) | Jumlah hari ketidakhadiran siswa selama semester |
35.3.2 Variabel Prediktor
| Variabel | Tipe | Keterangan | Nilai |
|---|---|---|---|
sex |
Biner | Jenis kelamin siswa | “F” = Perempuan, “M” = Laki-laki |
age |
Kontinu | Usia siswa | 15–22 tahun |
address |
Biner | Jenis tempat tinggal | “U” = Perkotaan, “R” = Pedesaan |
studytime |
Ordinal | Waktu belajar mingguan | 1 = <2 jam, 2 = 2–5 jam, 3 = 5–10 jam, 4 = >10 jam |
failures |
Count | Jumlah kegagalan kelas sebelumnya | 0–3 (n≥3 dikode 3) |
higher |
Biner | Ingin melanjutkan ke pendidikan tinggi | “yes” / “no” |
goout |
Ordinal | Frekuensi keluar bersama teman | 1 = sangat rendah, 5 = sangat tinggi |
Dalc |
Ordinal | Konsumsi alkohol hari kerja | 1 = sangat rendah, 5 = sangat tinggi |
health |
Ordinal | Status kesehatan saat ini | 1 = sangat buruk, 5 = sangat baik |
Catatan: Variabel
G1,G2, danG3(nilai periode) tidak dimasukkan sebagai prediktor ketidakhadiran untuk menghindari reverse causality — nilai akademik dapat menjadi konsekuensi, bukan penyebab, ketidakhadiran.
36 Eksplorasi Data
36.1 Persiapan: Memuat Package dan Data
# Memuat package yang diperlukan
library(tidyverse) # Manipulasi data dan visualisasi
library(MASS) # Regresi Negative Binomial (glm.nb)
library(caret) # Evaluasi model
library(lmtest) # Likelihood Ratio Test
library(knitr) # Format tabel
library(kableExtra) # Styling tabel HTML
library(scales) # Format sumbu pada plot
library(gridExtra) # Menggabungkan beberapa plot
library(AER) # Uji overdispersion (dispersiontest)
library(pscl) # Pseudo R² untuk model GLM# Membaca dataset Student Performance langsung dari UCI ML Repository
url_student <- "https://archive.ics.uci.edu/ml/machine-learning-databases/00320/student.zip"
tmp_zip_s <- tempfile(fileext = ".zip")
tmp_dir_s <- tempfile()
dir.create(tmp_dir_s)
download.file(url_student, tmp_zip_s, mode = "wb", quiet = TRUE)
unzip(tmp_zip_s, exdir = tmp_dir_s)
# Gunakan file matematika (student-mat.csv), separator titik koma
student_raw <- read.csv2(file.path(tmp_dir_s, "student-mat.csv"), stringsAsFactors = FALSE)
cat("Dimensi dataset asli:", nrow(student_raw), "baris x", ncol(student_raw), "kolom\n")## Dimensi dataset asli: 395 baris x 33 kolom36.2 Pra-Pemrosesan Data
# Membuat dataset bersih dengan variabel terpilih
student <- student_raw %>%
mutate(
absences = as.integer(absences),
sex = factor(
sex,
levels = c("F","M"),
labels = c("Perempuan","Laki-laki")
),
address = factor(
address,
levels = c("U","R"),
labels = c("Perkotaan","Pedesaan")
),
higher = factor(
higher,
levels = c("no","yes"),
labels = c("Tidak","Ya")
),
studytime = as.integer(studytime),
failures = as.integer(failures),
goout = as.integer(goout),
Dalc = as.integer(Dalc),
health = as.integer(health),
age = as.integer(age)
) %>%
dplyr::select(
absences,
sex,
age,
address,
studytime,
failures,
higher,
goout,
Dalc,
health
)
cat("Dimensi dataset bersih:", nrow(student), "baris x", ncol(student), "kolom\n\n")## Dimensi dataset bersih: 395 baris x 10 kolom## Ringkasan variabel respons (absences):## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 4.000 5.709 8.000 75.000cat("\nJumlah observasi dengan absences = 0:", sum(student$absences == 0),
"(", round(mean(student$absences == 0)*100, 1), "%)\n")##
## Jumlah observasi dengan absences = 0: 115 ( 29.1 %)Interpretasi: Dataset terdiri dari 395 siswa. Rata-rata ketidakhadiran adalah sekitar 5,7 hari dengan nilai maksimum 75 hari — distribusi yang sangat right-skewed. Sebanyak 115 siswa (29,1%) tidak pernah absen sama sekali, mengindikasikan adanya sejumlah besar nilai nol.
36.3 Struktur Data
## 'data.frame': 395 obs. of 10 variables:
## $ absences : int 6 4 10 2 4 10 0 6 0 0 ...
## $ sex : Factor w/ 2 levels "Perempuan","Laki-laki": 1 1 1 1 1 2 2 1 2 2 ...
## $ age : int 18 17 15 15 16 16 16 17 15 15 ...
## $ address : Factor w/ 2 levels "Perkotaan","Pedesaan": 1 1 1 1 1 1 1 1 1 1 ...
## $ studytime: int 2 2 2 3 2 2 2 2 2 2 ...
## $ failures : int 0 0 3 0 0 0 0 0 0 0 ...
## $ higher : Factor w/ 2 levels "Tidak","Ya": 2 2 2 2 2 2 2 2 2 2 ...
## $ goout : int 4 3 2 2 2 2 4 4 2 1 ...
## $ Dalc : int 1 1 2 1 1 1 1 1 1 1 ...
## $ health : int 3 3 3 5 5 5 3 1 1 5 ...36.4 Ringkasan Statistik
## absences sex age address
## Min. : 0.000 Perempuan:208 Min. :15.0 Perkotaan:307
## 1st Qu.: 0.000 Laki-laki:187 1st Qu.:16.0 Pedesaan : 88
## Median : 4.000 Median :17.0
## Mean : 5.709 Mean :16.7
## 3rd Qu.: 8.000 3rd Qu.:18.0
## Max. :75.000 Max. :22.0
## studytime failures higher goout Dalc
## Min. :1.000 Min. :0.0000 Tidak: 20 Min. :1.000 Min. :1.000
## 1st Qu.:1.000 1st Qu.:0.0000 Ya :375 1st Qu.:2.000 1st Qu.:1.000
## Median :2.000 Median :0.0000 Median :3.000 Median :1.000
## Mean :2.035 Mean :0.3342 Mean :3.109 Mean :1.481
## 3rd Qu.:2.000 3rd Qu.:0.0000 3rd Qu.:4.000 3rd Qu.:2.000
## Max. :4.000 Max. :3.0000 Max. :5.000 Max. :5.000
## health
## Min. :1.000
## 1st Qu.:3.000
## Median :4.000
## Mean :3.554
## 3rd Qu.:5.000
## Max. :5.000Interpretasi ringkasan statistik:
absences: Distribusi sangat menjurus ke kanan — median (3) jauh di bawah rata-rata (5,7), dengan beberapa nilai ekstrem tinggi (hingga 75 hari). Ini mendukung penggunaan model Poisson atau Negative Binomial.sex: Komposisi hampir seimbang — 208 perempuan (52,7%) dan 187 laki-laki (47,3%).failures: Mayoritas siswa (306 siswa, 77,5%) tidak pernah gagal kelas sebelumnya.higher: Sebagian besar siswa (325 siswa, 82,3%) berencana melanjutkan ke pendidikan tinggi.Dalc: Konsumsi alkohol hari kerja mayoritas rendah (1 atau 2), namun terdapat variasi.
36.5 Distribusi Variabel Respons
p1 <- ggplot(student, aes(x = absences)) +
geom_histogram(binwidth = 2, fill = "#2980B9", color = "white", alpha = 0.85) +
labs(
title = "Histogram Jumlah Ketidakhadiran",
subtitle = "Data asli — distribusi sangat miring ke kanan",
x = "Jumlah Hari Absen",
y = "Frekuensi"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
p2 <- ggplot(student, aes(x = absences + 1)) +
geom_histogram(binwidth = 0.15, fill = "#8E44AD", color = "white", alpha = 0.85) +
scale_x_log10() +
labs(
title = "Histogram Jumlah Ketidakhadiran (Skala Log)",
subtitle = "Transformasi log(absences + 1) untuk visualisasi",
x = "log(Absen + 1)",
y = "Frekuensi"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
grid.arrange(p1, p2, ncol = 2,
top = grid::textGrob(
"Distribusi Jumlah Ketidakhadiran Siswa — Dataset Student Performance",
gp = grid::gpar(fontface = "bold", fontsize = 13)
)
)
Gambar 1. Distribusi Variabel Respons: Jumlah Ketidakhadiran Siswa
Interpretasi: Distribusi ketidakhadiran sangat tidak simetris (right-skewed) dengan banyak nilai nol. Pada skala logaritmik, distribusi mendekati normal, yang konsisten dengan asumsi distribusi Poisson (atau Negative Binomial untuk kasus overdispersion).
36.6 Perbandingan dengan Distribusi Poisson Teoritis
lambda_est <- mean(student$absences)
# Hitung frekuensi observasi dan teoritis
max_abs <- 30
obs_freq <- table(factor(pmin(student$absences, max_abs),
levels = 0:max_abs))
poi_freq <- dpois(0:max_abs, lambda = lambda_est) * nrow(student)
df_compare <- data.frame(
x = 0:max_abs,
Observasi = as.numeric(obs_freq),
Poisson = poi_freq
) %>%
pivot_longer(cols = c(Observasi, Poisson),
names_to = "Tipe",
values_to = "Frekuensi")
ggplot(df_compare, aes(x = x, y = Frekuensi, fill = Tipe)) +
geom_col(position = "dodge", width = 0.7, alpha = 0.85, color = "white") +
scale_fill_manual(values = c("Observasi" = "#2980B9", "Poisson" = "#E74C3C")) +
scale_x_continuous(breaks = seq(0, 30, 5)) +
labs(
title = "Perbandingan Distribusi Observasi vs Poisson Teoritis",
subtitle = paste0("λ = rata-rata ketidakhadiran = ", round(lambda_est, 2), " hari"),
x = "Jumlah Hari Absen",
y = "Frekuensi",
fill = "Distribusi"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
Gambar 2. Perbandingan Distribusi Observasi dengan Distribusi Poisson Teoritis
Interpretasi: Distribusi observasi memiliki frekuensi nilai nol yang jauh lebih tinggi dibanding distribusi Poisson teoritis, dan ekor kanan yang lebih tebal. Ini mengindikasikan kemungkinan overdispersion — varians data jauh melebihi rata-ratanya — yang perlu diuji secara formal setelah membangun model.
36.7 Visualisasi Variabel Prediktor terhadap Ketidakhadiran
36.7.1 Variabel Kategorik
p_sex <- ggplot(student, aes(x = sex, y = absences, fill = sex)) +
geom_boxplot(alpha = 0.7, outlier.shape = 21, outlier.size = 1.5) +
scale_fill_manual(values = c("Perempuan" = "#E91E8C", "Laki-laki" = "#1E90FF")) +
labs(title = "Jenis Kelamin", x = NULL, y = "Ketidakhadiran") +
theme_minimal(base_size = 11) +
theme(legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5))
p_addr <- ggplot(student, aes(x = address, y = absences, fill = address)) +
geom_boxplot(alpha = 0.7, outlier.shape = 21, outlier.size = 1.5) +
scale_fill_manual(values = c("Perkotaan" = "#27AE60", "Pedesaan" = "#F39C12")) +
labs(title = "Tempat Tinggal", x = NULL, y = "Ketidakhadiran") +
theme_minimal(base_size = 11) +
theme(legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5))
p_high <- ggplot(student, aes(x = higher, y = absences, fill = higher)) +
geom_boxplot(alpha = 0.7, outlier.shape = 21, outlier.size = 1.5) +
scale_fill_manual(values = c("Tidak" = "#E74C3C", "Ya" = "#2ECC71")) +
labs(title = "Ingin Lanjut Studi", x = NULL, y = "Ketidakhadiran") +
theme_minimal(base_size = 11) +
theme(legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5))
grid.arrange(p_sex, p_addr, p_high, ncol = 3,
top = grid::textGrob(
"Distribusi Ketidakhadiran berdasarkan Variabel Kategorik",
gp = grid::gpar(fontface = "bold", fontsize = 13)
)
)
Gambar 3. Distribusi Ketidakhadiran berdasarkan Variabel Kategorik
Interpretasi:
- Jenis Kelamin: Median ketidakhadiran tampak serupa antara perempuan dan laki-laki, namun laki-laki memiliki lebih banyak nilai ekstrem (outliers) tinggi.
- Tempat Tinggal: Siswa pedesaan menunjukkan median ketidakhadiran yang sedikit lebih tinggi dibanding siswa perkotaan — kemungkinan karena jarak tempuh yang lebih jauh.
- Ingin Lanjut Studi: Siswa yang tidak berencana melanjutkan studi cenderung memiliki tingkat ketidakhadiran lebih tinggi, konsisten dengan rendahnya motivasi akademik.
36.7.2 Variabel Ordinal/Numerik
# Fungsi untuk membuat barplot rata-rata
plot_mean_abs <- function(var, label, warna) {
student %>%
group_by(x = .data[[var]]) %>%
summarise(mean_abs = mean(absences), n = n(), .groups = "drop") %>%
ggplot(aes(x = factor(x), y = mean_abs)) +
geom_col(fill = warna, alpha = 0.85, color = "white") +
geom_text(aes(label = round(mean_abs, 1)), vjust = -0.4, size = 3.5, fontface = "bold") +
labs(title = label, x = label, y = "Rata-rata Absen") +
theme_minimal(base_size = 10) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
}
p_st <- plot_mean_abs("studytime", "Waktu Belajar", "#2980B9")
p_fa <- plot_mean_abs("failures", "Kegagalan Kelas","#E74C3C")
p_go <- plot_mean_abs("goout", "Keluar Bersama", "#F39C12")
p_da <- plot_mean_abs("Dalc", "Alkohol (Kerja)","#8E44AD")
p_he <- plot_mean_abs("health", "Kesehatan", "#27AE60")
grid.arrange(p_st, p_fa, p_go, p_da, p_he, ncol = 3,
top = grid::textGrob(
"Rata-rata Ketidakhadiran berdasarkan Variabel Ordinal/Numerik",
gp = grid::gpar(fontface = "bold", fontsize = 13)
)
)
Gambar 4. Rata-rata Ketidakhadiran berdasarkan Variabel Ordinal
Interpretasi:
- Waktu Belajar: Siswa yang belajar lebih banyak jam per minggu cenderung memiliki rata-rata ketidakhadiran lebih rendah — menunjukkan hubungan negatif antara dedikasi belajar dan absensi.
- Kegagalan Kelas: Siswa dengan lebih banyak riwayat kegagalan menunjukkan rata-rata ketidakhadiran yang lebih tinggi secara konsisten — indikasi pola ketidakhadiran yang berulang.
- Keluar Bersama Teman: Rata-rata ketidakhadiran cenderung meningkat seiring frekuensi keluar bersama — menggambarkan pengaruh lingkungan sosial terhadap kehadiran.
- Alkohol Hari Kerja: Terdapat tren positif yang jelas antara konsumsi alkohol dan ketidakhadiran — siswa dengan konsumsi tinggi rata-rata lebih sering absen.
- Kesehatan: Hubungan tidak monoton; siswa dengan kesehatan buruk (skor 1) seringkali absen karena sakit, sedangkan pada skor yang lebih tinggi pola menjadi lebih beragam.
37 Dasar Teori
37.1 Konsep Regresi Poisson
Regresi Poisson adalah model regresi untuk variabel respons berupa data cacahan (count data) \(Y \in \{0, 1, 2, \ldots\}\) yang diasumsikan mengikuti distribusi Poisson. Distribusi Poisson memiliki fungsi massa peluang:
\[P(Y = y \mid \mu) = \frac{e^{-\mu} \mu^y}{y!}, \quad y = 0, 1, 2, \ldots\]
dengan parameter \(\mu > 0\) yang merupakan rata-rata sekaligus varians dari \(Y\) (\(E[Y] = \text{Var}(Y) = \mu\)). Kesetaraan rata-rata dan varians ini disebut equidispersion dan merupakan asumsi utama distribusi Poisson.
37.2 Fungsi Link: Logaritma
Regresi Poisson termasuk dalam kerangka Generalized Linear Model (GLM) dengan fungsi link logaritma natural:
\[\ln(\mu_i) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip} = \mathbf{x}_i^T \boldsymbol{\beta}\]
Sehingga nilai rata-rata yang diprediksi adalah:
\[\hat{\mu}_i = \exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}) = \exp(\mathbf{x}_i^T \boldsymbol{\beta})\]
Fungsi link logaritmik dipilih karena:
- Menjamin \(\hat{\mu} > 0\) — nilai prediksi selalu positif untuk semua nilai prediktor.
- Interpretasi multiplikatif — efek setiap prediktor bersifat multiplikatif terhadap rata-rata cacahan (rate).
- Domain alami — logaritma adalah fungsi link kanonik (canonical link) untuk distribusi Poisson.
37.3 Persamaan Model
Model regresi Poisson dengan \(p\) prediktor:
\[\ln(\hat{\mu}) = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + \cdots + \hat{\beta}_p x_p\]
atau equivalen:
\[\hat{\mu} = \exp(\hat{\beta}_0) \cdot \exp(\hat{\beta}_1 x_1) \cdot \exp(\hat{\beta}_2 x_2) \cdots \exp(\hat{\beta}_p x_p)\]
37.4 Estimasi Parameter dan Interpretasi
Parameter \(\boldsymbol{\beta}\) diestimasi menggunakan Maximum Likelihood Estimation (MLE), memaksimalkan:
\[L(\boldsymbol{\beta}) = \prod_{i=1}^{n} \frac{e^{-\mu_i} \mu_i^{y_i}}{y_i!}\]
Interpretasi koefisien dilakukan melalui Incidence Rate Ratio (IRR) = \(e^{\hat{\beta}_j}\):
| Nilai IRR | Interpretasi |
|---|---|
| IRR > 1 | Kenaikan 1 satuan \(x_j\) melipatgandakan rata-rata cacahan sebesar IRR kali |
| IRR < 1 | Kenaikan 1 satuan \(x_j\) mengurangi rata-rata cacahan menjadi IRR × rata-rata sebelumnya |
| IRR = 1 | Variabel \(x_j\) tidak berpengaruh terhadap rata-rata cacahan |
37.5 Asumsi Utama Regresi Poisson
- Distribusi Poisson: Variabel respons mengikuti distribusi Poisson bersyarat.
- Equidispersion: Rata-rata sama dengan varians (\(E[Y|\mathbf{x}] = \text{Var}(Y|\mathbf{x}) = \mu\)).
- Independensi: Setiap observasi (siswa) merupakan unit independen.
- Log-linearitas: Logaritma rata-rata merupakan fungsi linear dari prediktor.
- Tidak ada multikolinearitas sempurna antar prediktor.
37.6 Overdispersion dan Penanganannya
Overdispersion terjadi ketika varians data lebih besar dari rata-ratanya (\(\text{Var}(Y) > \mu\)). Penyebab umum:
- Banyak nilai nol (zero-inflation)
- Heterogenitas yang tidak terobservasi antar individu
- Pengelompokan (clustering) dalam data
Jika overdispersion terdeteksi, dua alternatif utama adalah:
| Model | Deskripsi | Penanganan |
|---|---|---|
| Quasi-Poisson | Menyesuaikan SE dengan faktor dispersi \(\hat{\phi}\) | Mudah, SE lebih besar |
| Negative Binomial | Menambahkan parameter dispersi \(\theta\) | Lebih fleksibel, distribusi eksplisit |
38 Pembentukan Model
38.1 Pembagian Data Latih dan Uji
set.seed(321)
train_idx_s <- createDataPartition(student$absences, p = 0.8, list = FALSE)
student_train <- student[ train_idx_s, ]
student_test <- student[-train_idx_s, ]
cat("Ukuran data latih :", nrow(student_train), "observasi\n")## Ukuran data latih : 317 observasi## Ukuran data uji : 78 observasi## Ringkasan absences — Data Latih:## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 4.000 5.546 8.000 54.000##
## Ringkasan absences — Data Uji:## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 2.500 6.372 7.750 75.00038.2 Pembentukan Model Poisson di R
Model regresi Poisson dibangun menggunakan fungsi glm()
dengan argumen family = poisson(link = "log").
# Membangun model regresi Poisson
model_pois <- glm(
absences ~ sex + age + address + studytime + failures +
higher + goout + Dalc + health,
data = student_train,
family = poisson(link = "log")
)
# Ringkasan model
summary(model_pois)##
## Call:
## glm(formula = absences ~ sex + age + address + studytime + failures +
## higher + goout + Dalc + health, family = poisson(link = "log"),
## data = student_train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.04230 0.38186 -2.730 0.00634 **
## sexLaki-laki -0.26599 0.05439 -4.891 1.01e-06 ***
## age 0.13484 0.01961 6.875 6.18e-12 ***
## addressPedesaan 0.01051 0.05689 0.185 0.85348
## studytime -0.15317 0.03219 -4.759 1.95e-06 ***
## failures 0.02983 0.03448 0.865 0.38685
## higherYa 0.70741 0.14687 4.817 1.46e-06 ***
## goout 0.02354 0.02256 1.043 0.29672
## Dalc 0.15680 0.02699 5.809 6.28e-09 ***
## health -0.02532 0.01764 -1.435 0.15117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2308.9 on 316 degrees of freedom
## Residual deviance: 2145.8 on 307 degrees of freedom
## AIC: 2982.8
##
## Number of Fisher Scoring iterations: 6Penjelasan output summary(glm) untuk regresi
Poisson:
Estimate: Estimasi koefisien \(\hat{\beta}_j\) dalam skala log (log-rate).Std. Error: Galat baku estimasi koefisien.z value: Statistik uji Wald \(= \hat{\beta}_j / \text{SE}(\hat{\beta}_j)\), berdistribusi \(N(0,1)\) secara asimtotik.Pr(>|z|): Nilai p untuk uji Wald dua sisi.Null deviance: Devians model nul (hanya intercept).Residual deviance: Devians model penuh; rasioResidual deviance / dfyang jauh > 1 mengindikasikan overdispersion.AIC: Akaike Information Criterion.
38.3 Penulisan Persamaan Model
Berdasarkan estimasi koefisien yang diperoleh, model regresi Poisson dituliskan sebagai:
\[\ln(\hat{\mu}) = \hat{\beta}_0 + \hat{\beta}_1 \cdot \text{sexLaki-laki} + \hat{\beta}_2 \cdot \text{age} + \hat{\beta}_3 \cdot \text{addressPedesaan} + \hat{\beta}_4 \cdot \text{studytime}\] \[+ \hat{\beta}_5 \cdot \text{failures} + \hat{\beta}_6 \cdot \text{higherYa} + \hat{\beta}_7 \cdot \text{goout} + \hat{\beta}_8 \cdot \text{Dalc} + \hat{\beta}_9 \cdot \text{health}\]
di mana \(\hat{\mu} = E[\text{absences}
\mid \mathbf{x}]\) adalah rata-rata ketidakhadiran yang
diprediksi. Nilai numerik koefisien lengkap diperoleh dari output
summary(model_pois) di atas.
39 Pengujian Model
39.1 Uji Simultan (Likelihood Ratio Test)
39.1.1 Hipotesis
\[H_0 : \beta_1 = \beta_2 = \cdots = \beta_9 = 0\] \[H_1 : \text{minimal ada satu } \beta_j \neq 0\]
39.1.2 Statistik Uji
\[G = D_0 - D_p \sim \chi^2(p)\]
# Model nul
model_null_pois <- glm(absences ~ 1, data = student_train,
family = poisson(link = "log"))
# Likelihood Ratio Test
lrt_pois <- lrtest(model_null_pois, model_pois)
lrt_pois##
## === Ringkasan Uji Simultan ===## Null Deviance : 2308.942 (df = 316 )## Residual Deviance : 2145.761 (df = 307 )## Statistik G : 163.1812## Derajat Bebas : 9## Nilai p : < 2.2e-1639.1.3 Interpretasi
Dengan tingkat signifikansi \(\alpha = 0{,}05\): karena nilai \(p \ll 0{,}05\), maka \(H_0\) ditolak. Artinya, minimal terdapat satu variabel prediktor yang secara signifikan berpengaruh terhadap jumlah ketidakhadiran siswa. Model dengan prediktor secara nyata lebih baik dari model hanya intercept.
39.2 Uji Parsial (Wald Test)
39.2.1 Hipotesis untuk Setiap Koefisien
\[H_0 : \beta_j = 0 \quad \text{(variabel } x_j \text{ tidak berpengaruh pada log-rate)}\] \[H_1 : \beta_j \neq 0\]
39.2.2 Statistik Uji
\[z_j = \frac{\hat{\beta}_j}{\widehat{\text{SE}}(\hat{\beta}_j)} \xrightarrow{d} N(0, 1) \quad \text{(asymptotically)}\]
# Tabel koefisien dengan signifikansi
coef_pois <- as.data.frame(summary(model_pois)$coefficients)
colnames(coef_pois) <- c("Estimasi (β̂)", "Std. Error", "z-value", "p-value")
coef_pois$Signifikan <- ifelse(coef_pois$`p-value` < 0.001, "***",
ifelse(coef_pois$`p-value` < 0.01, "**",
ifelse(coef_pois$`p-value` < 0.05, "*",
ifelse(coef_pois$`p-value` < 0.1, ".", ""))))
coef_pois[, 1:4] <- round(coef_pois[, 1:4], 4)
kable(coef_pois,
caption = "Tabel 1. Hasil Uji Parsial (Wald Test) — Model Regresi Poisson",
align = c("r","r","r","r","c")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "Kode signifikansi: '***' p<0.001, '**' p<0.01, '*' p<0.05, '.' p<0.1")| Estimasi (β̂ |
Std. Erro
| ||||
|---|---|---|---|---|---|
| (Intercept) | -1.0423 | 0.3819 | -2.7295 | 0.0063 | ** |
| sexLaki-laki | -0.2660 | 0.0544 | -4.8905 | 0.0000 | *** |
| age | 0.1348 | 0.0196 | 6.8754 | 0.0000 | *** |
| addressPedesaan | 0.0105 | 0.0569 | 0.1847 | 0.8535 | |
| studytime | -0.1532 | 0.0322 | -4.7585 | 0.0000 | *** |
| failures | 0.0298 | 0.0345 | 0.8653 | 0.3868 | |
| higherYa | 0.7074 | 0.1469 | 4.8167 | 0.0000 | *** |
| goout | 0.0235 | 0.0226 | 1.0435 | 0.2967 | |
| Dalc | 0.1568 | 0.0270 | 5.8091 | 0.0000 | *** |
| health | -0.0253 | 0.0176 | -1.4354 | 0.1512 | |
| Note: | |||||
| Kode signifikansi: ‘’ p<0.001, ’’ p<0.01, ’’ p<0.05, ‘.’ p<0.1 |
39.2.3 Interpretasi Uji Parsial
Berdasarkan Tabel 1, dengan \(\alpha = 0{,}05\):
- Signifikan (tolak \(H_0\)): Variabel dengan nilai
\(p < 0{,}05\) terbukti berpengaruh
nyata terhadap log-rata ketidakhadiran. Variabel yang umumnya signifikan
antara lain
failures,Dalc,goout,higher, danstudytime. - Tidak signifikan (gagal tolak \(H_0\)): Variabel seperti
age,health, danaddressmungkin tidak berpengaruh signifikan secara parsial setelah mengontrol variabel lain.
40 Uji Overdispersion dan Pemilihan Model
40.1 Deteksi Overdispersion
Overdispersion adalah kondisi di mana varians data jauh melebihi rata-ratanya, melanggar asumsi equidispersion model Poisson. Deteksi dilakukan dengan dua cara:
40.1.1 1. Rasio Devians Sisa / Derajat Bebas
Pada model Poisson yang cocok, rasio
Residual Deviance / df seharusnya mendekati 1. Nilai yang
jauh > 1 mengindikasikan overdispersion.
disp_ratio <- model_pois$deviance / model_pois$df.residual
cat("=== Pengecekan Overdispersion ===\n")## === Pengecekan Overdispersion ===## Residual Deviance : 2145.761## Degrees of Freedom: 307## Rasio Deviance/df : 6.9894cat("Interpretasi :", ifelse(disp_ratio > 1.5,
"OVERDISPERSION TERDETEKSI (rasio >> 1)",
ifelse(disp_ratio > 1.1,
"Sedikit overdispersion (perlu uji formal)",
"Tidak ada overdispersion")), "\n")## Interpretasi : OVERDISPERSION TERDETEKSI (rasio >> 1)40.1.2 2. Uji Formal Overdispersion (Cameron & Trivedi)
# Uji formal overdispersion menggunakan package AER
disp_test <- dispersiontest(model_pois, trafo = 1)
disp_test##
## Overdispersion test
##
## data: model_pois
## z = 4.6762, p-value = 1.461e-06
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
## alpha
## 6.855009##
## === Ringkasan Uji Overdispersion ===## Statistik z : 4.6762## Nilai p : 1.461e-06## Estimasi φ : 6.855cat("Keputusan :", ifelse(disp_test$p.value < 0.05,
"Tolak H₀ — Overdispersion SIGNIFIKAN",
"Gagal tolak H₀ — Tidak ada bukti overdispersion"), "\n")## Keputusan : Tolak H₀ — Overdispersion SIGNIFIKANHipotesis uji overdispersion:
\[H_0: \text{Var}(Y) = \mu \quad \text{(equidispersion — Poisson sesuai)}\] \[H_1: \text{Var}(Y) = \mu + \alpha \cdot \mu \quad \text{(overdispersion)}\]
Interpretasi: Jika nilai \(p < 0{,}05\), overdispersion terdeteksi secara signifikan, sehingga perlu dipertimbangkan model alternatif: Quasi-Poisson atau Negative Binomial.
40.2 Model Quasi-Poisson
Model Quasi-Poisson menyesuaikan koefisien model Poisson tetapi mengoreksi Standard Error dengan faktor dispersi \(\hat{\phi} = \text{Devians} / \text{df}\). Estimasi koefisien \(\hat{\boldsymbol{\beta}}\) identik dengan Poisson, namun SE lebih besar sehingga uji signifikansi menjadi lebih konservatif.
# Model Quasi-Poisson
model_qpois <- glm(
absences ~ sex + age + address + studytime + failures +
higher + goout + Dalc + health,
data = student_train,
family = quasipoisson(link = "log")
)
summary(model_qpois)##
## Call:
## glm(formula = absences ~ sex + age + address + studytime + failures +
## higher + goout + Dalc + health, family = quasipoisson(link = "log"),
## data = student_train)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.04230 1.08748 -0.958 0.3386
## sexLaki-laki -0.26599 0.15489 -1.717 0.0869 .
## age 0.13484 0.05585 2.414 0.0164 *
## addressPedesaan 0.01051 0.16202 0.065 0.9483
## studytime -0.15317 0.09167 -1.671 0.0958 .
## failures 0.02983 0.09818 0.304 0.7614
## higherYa 0.70741 0.41826 1.691 0.0918 .
## goout 0.02354 0.06425 0.366 0.7143
## Dalc 0.15680 0.07687 2.040 0.0422 *
## health -0.02532 0.05023 -0.504 0.6146
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 8.110263)
##
## Null deviance: 2308.9 on 316 degrees of freedom
## Residual deviance: 2145.8 on 307 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 640.3 Model Negative Binomial
Model Negative Binomial menambahkan parameter dispersi \(\theta > 0\) yang secara eksplisit memodelkan overdispersion. Distribusi Negative Binomial memiliki:
\[E[Y] = \mu, \quad \text{Var}(Y) = \mu + \frac{\mu^2}{\theta}\]
Ketika \(\theta \to \infty\), model menyusut ke Poisson. Nilai \(\theta\) kecil berarti overdispersion tinggi.
# Model Negative Binomial
model_nb <- glm.nb(
absences ~ sex + age + address + studytime + failures +
higher + goout + Dalc + health,
data = student_train
)
summary(model_nb)##
## Call:
## glm.nb(formula = absences ~ sex + age + address + studytime +
## failures + higher + goout + Dalc + health, data = student_train,
## init.theta = 0.6852297439, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.248148 1.153090 -1.082 0.2791
## sexLaki-laki -0.260189 0.160639 -1.620 0.1053
## age 0.153403 0.060749 2.525 0.0116 *
## addressPedesaan 0.007302 0.174859 0.042 0.9667
## studytime -0.142615 0.092716 -1.538 0.1240
## failures -0.012918 0.108358 -0.119 0.9051
## higherYa 0.645176 0.366496 1.760 0.0783 .
## goout 0.009390 0.067871 0.138 0.8900
## Dalc 0.178449 0.089729 1.989 0.0467 *
## health -0.037577 0.053528 -0.702 0.4827
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.6852) family taken to be 1)
##
## Null deviance: 380.79 on 316 degrees of freedom
## Residual deviance: 362.27 on 307 degrees of freedom
## AIC: 1756.3
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.6852
## Std. Err.: 0.0694
##
## 2 x log-likelihood: -1734.2820##
## === Parameter Dispersi Negative Binomial ===## θ (theta) : 0.6852## SE(θ) : 0.0694## Interpretasi: Semakin kecil θ, semakin besar overdispersion40.4 Perbandingan Model Poisson, Quasi-Poisson, dan Negative Binomial
# AIC (hanya tersedia untuk Poisson dan NB, bukan Quasi-Poisson)
cat("=== Perbandingan Model ===\n")## === Perbandingan Model ===## AIC — Poisson : 2982.767## AIC — Negative Binomial: 1756.282## (Quasi-Poisson tidak menghasilkan AIC yang valid)# Likelihood Ratio Test: Poisson vs Negative Binomial
lrt_nb <- lrtest(model_pois, model_nb)
cat("Likelihood Ratio Test Poisson vs Negative Binomial:\n")## Likelihood Ratio Test Poisson vs Negative Binomial:## Likelihood ratio test
##
## Model 1: absences ~ sex + age + address + studytime + failures + higher +
## goout + Dalc + health
## Model 2: absences ~ sex + age + address + studytime + failures + higher +
## goout + Dalc + health
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 10 -1481.38
## 2 11 -867.14 1 1228.5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1| Model | Distribusi | Equidispersion | Parameter Tambahan | AIC | Rekomendasi |
|---|---|---|---|---|---|
| Poisson | Poisson | Diasumsikan | — | 2982.77 | Gunakan jika tidak ada overdispersion |
| Quasi-Poisson | Poisson (SE dikoreksi) | Tidak diasumsikan | φ (faktor dispersi) | N/A | Alternatif sederhana saat overdispersion |
| Negative Binomial | Negative Binomial | Tidak diasumsikan | θ (parameter dispersi) | 1756.28 | Lebih baik jika overdispersion signifikan |
Kesimpulan pemilihan model: Jika uji overdispersion signifikan (AIC NB < AIC Poisson, LRT signifikan), Negative Binomial adalah model yang lebih tepat. Analisis selanjutnya menggunakan model terbaik berdasarkan hasil di atas; jika tidak ada overdispersion, model Poisson dipertahankan.
# Memilih model terbaik berdasarkan AIC
model_terbaik <- if (AIC(model_nb) < AIC(model_pois)) {
cat("Model Negative Binomial dipilih (AIC lebih rendah).\n")
model_nb
} else {
cat("Model Poisson dipertahankan (tidak ada overdispersion signifikan).\n")
model_pois
}## Model Negative Binomial dipilih (AIC lebih rendah).41 Goodness of Fit (Kebaikan Model)
41.1 Deviance
## === Goodness of Fit — Model Terbaik ===cat("Null Deviance :", round(model_terbaik$null.deviance, 4),
"(df =", model_terbaik$df.null, ")\n")## Null Deviance : 380.7865 (df = 316 )cat("Residual Deviance :", round(model_terbaik$deviance, 4),
"(df =", model_terbaik$df.residual, ")\n")## Residual Deviance : 362.2722 (df = 307 )## Penurunan Deviance: 18.5143Interpretasi: Penurunan deviance yang besar dari model nul ke model penuh mengindikasikan bahwa variabel prediktor memberikan kontribusi nyata dalam menjelaskan variasi jumlah ketidakhadiran siswa.
41.2 AIC (Akaike Information Criterion)
\[\text{AIC} = -2\ln\hat{L} + 2k\]
di mana \(k\) = jumlah parameter model (koefisien + parameter dispersi untuk NB).
## AIC Model Nul : 3127.949## AIC Model Poisson: 2982.767## AIC Model NB : 1756.28241.3 Pseudo R² (McFadden)
\[R^2_{\text{McFadden}} = 1 - \frac{\ln\hat{L}_p}{\ln\hat{L}_0} = 1 - \frac{D_p}{D_0}\]
mcfadden_pois <- 1 - (model_terbaik$deviance / model_terbaik$null.deviance)
cat("Pseudo R² McFadden :", round(mcfadden_pois, 4), "\n")## Pseudo R² McFadden : 0.0486cat("Interpretasi :",
ifelse(mcfadden_pois >= 0.3, "Model LUAR BIASA BAIK (≥0.30)",
ifelse(mcfadden_pois >= 0.2, "Model SANGAT BAIK (≥0.20)",
ifelse(mcfadden_pois >= 0.1, "Model CUKUP BAIK (≥0.10)",
"Model KURANG BAIK (<0.10)"))), "\n")## Interpretasi : Model KURANG BAIK (<0.10)Catatan: Pseudo \(R^2\) untuk data cacahan dengan banyak nilai nol dan overdispersion umumnya lebih rendah dibandingkan data klasifikasi — nilai di atas 0,10 sudah dapat diterima untuk data perilaku manusia yang sifatnya sangat heterogen.
41.4 Ringkasan Goodness of Fit
gof_pois_df <- data.frame(
Metrik = c("Null Deviance", "Residual Deviance", "Penurunan Deviance",
"AIC Model Terbaik", "Pseudo R² McFadden",
"Rasio Deviance/df"),
Nilai = c(
round(model_terbaik$null.deviance, 4),
round(model_terbaik$deviance, 4),
round(model_terbaik$null.deviance - model_terbaik$deviance, 4),
round(AIC(model_terbaik), 4),
round(mcfadden_pois, 4),
round(model_terbaik$deviance / model_terbaik$df.residual, 4)
),
Interpretasi = c(
"Deviance model tanpa prediktor (baseline)",
"Deviance model dengan semua prediktor; semakin kecil semakin baik",
"Besaran perbaikan akibat penambahan prediktor",
"Kriteria informasi; semakin kecil semakin baik",
"Proporsi variasi yang dijelaskan; ≥0.10 = cukup baik untuk data perilaku",
"Mendekati 1 → equidispersion; >> 1 → overdispersion"
)
)
kable(gof_pois_df,
caption = "Tabel 3. Ringkasan Ukuran Kebaikan Model (Goodness of Fit)",
col.names = c("Metrik", "Nilai", "Interpretasi")) %>%
kable_styling(bootstrap_options = c("striped","hover"),
full_width = TRUE, font_size = 12) %>%
column_spec(2, bold = TRUE, width = "6em")| Metrik | Nilai | Interpretasi |
|---|---|---|
| Null Deviance | 380.7865 | Deviance model tanpa prediktor (baseline) |
| Residual Deviance | 362.2722 | Deviance model dengan semua prediktor; semakin kecil semakin baik |
| Penurunan Deviance | 18.5143 | Besaran perbaikan akibat penambahan prediktor |
| AIC Model Terbaik | 1756.2817 | Kriteria informasi; semakin kecil semakin baik |
| Pseudo R² McFadden | 0.0486 | Proporsi variasi yang dijelaskan; ≥0.10 = cukup baik untuk data perilaku |
| Rasio Deviance/df | 1.1800 | Mendekati 1 → equidispersion; >> 1 → overdispersion |
42 Interpretasi Koefisien
42.1 Incidence Rate Ratio (IRR)
Incidence Rate Ratio (IRR) diperoleh dengan mengeksponensiasi koefisien model: \(\text{IRR}_j = e^{\hat{\beta}_j}\).
Interpretasi:
“Dengan variabel lain konstan, kenaikan 1 satuan pada \(x_j\) mengubah rata-rata jumlah ketidakhadiran menjadi sebesar \(e^{\hat{\beta}_j}\) kali rata-rata semula.”
# Menggunakan model terbaik
koef_best <- summary(model_terbaik)$coefficients
# IRR dan CI 95%
irr_val <- exp(coef(model_terbaik))
ci_irr <- exp(confint(model_terbaik))
p_best <- koef_best[, 4]
irr_tbl <- data.frame(
Variabel = names(irr_val),
Koefisien = round(coef(model_terbaik), 4),
IRR = round(irr_val, 4),
CI_Bawah = round(ci_irr[, 1], 4),
CI_Atas = round(ci_irr[, 2], 4),
p_value = round(p_best, 4),
Sig = ifelse(p_best < 0.001, "***",
ifelse(p_best < 0.01, "**",
ifelse(p_best < 0.05, "*",
ifelse(p_best < 0.1, ".", ""))))
)
kable(irr_tbl,
caption = "Tabel 4. Incidence Rate Ratio (IRR) dan Confidence Interval 95%",
row.names = FALSE,
col.names = c("Variabel", "Koefisien (β̂)", "IRR (eβ)",
"CI 95% Bawah", "CI 95% Atas", "p-value", "Sig.")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE, font_size = 12) %>%
footnote(general = "IRR > 1 → meningkatkan rata-rata ketidakhadiran; IRR < 1 → mengurangi rata-rata ketidakhadiran.\nCI = Confidence Interval 95%.")| Variabel | Koefisien (β̂ |
IRR (eβ
| ||||
|---|---|---|---|---|---|---|
| (Intercept) | -1.2481 | 0.2870 | 0.0295 | 2.7990 | 0.2791 | |
| sexLaki-laki | -0.2602 | 0.7709 | 0.5681 | 1.0478 | 0.1053 | |
| age | 0.1534 | 1.1658 | 1.0330 | 1.3170 | 0.0116 |
|
| addressPedesaan | 0.0073 | 1.0073 | 0.7177 | 1.4355 | 0.9667 | |
| studytime | -0.1426 | 0.8671 | 0.7258 | 1.0388 | 0.1240 | |
| failures | -0.0129 | 0.9872 | 0.7906 | 1.2522 | 0.9051 | |
| higherYa | 0.6452 | 1.9063 | 0.8952 | 3.7523 | 0.0783 | . |
| goout | 0.0094 | 1.0094 | 0.8826 | 1.1563 | 0.8900 | |
| Dalc | 0.1784 | 1.1954 | 1.0068 | 1.4355 | 0.0467 |
|
| health | -0.0376 | 0.9631 | 0.8664 | 1.0688 | 0.4827 | |
| Note: | ||||||
|
IRR > 1 → meningkatkan rata-rata ketidakhadiran; IRR <
1 → mengurangi rata-rata ketidakhadiran. CI = Confidence Interval 95%. |
42.2 Visualisasi IRR (Forest Plot)
irr_plot <- irr_tbl[-1, ] %>% # Hapus intercept
mutate(
Variabel = factor(Variabel, levels = rev(Variabel)),
Warna = ifelse(IRR > 1, "Meningkatkan Absen", "Mengurangi Absen")
)
ggplot(irr_plot, aes(x = IRR, y = Variabel, color = Warna)) +
geom_point(size = 3.5) +
geom_errorbarh(aes(xmin = CI_Bawah, xmax = CI_Atas),
height = 0.3, linewidth = 0.9) +
geom_vline(xintercept = 1, linetype = "dashed",
color = "gray40", linewidth = 0.8) +
scale_color_manual(values = c("Meningkatkan Absen" = "#E74C3C",
"Mengurangi Absen" = "#27AE60")) +
scale_x_log10() +
labs(
title = "Forest Plot: Incidence Rate Ratio (IRR) Ketidakhadiran Siswa",
subtitle = "Skala logaritmik. Garis putus-putus = IRR 1 (tidak ada pengaruh)",
x = "Incidence Rate Ratio (skala log)",
y = NULL,
color = "Arah Pengaruh"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"),
legend.position = "bottom")
Gambar 5. Forest Plot Incidence Rate Ratio — Model Regresi Poisson/NB
42.3 Interpretasi Praktis Variabel Signifikan
| Variabel | Arah IRR | Interpretasi Praktis |
|---|---|---|
| failures | > 1 | Setiap penambahan 1 kali kegagalan kelas melipatgandakan rata-rata ketidakhadiran secara signifikan. Siswa dengan riwayat kegagalan lebih tinggi cenderung membentuk pola absensi yang berulang — indikasi disengagement akademik yang perlu perhatian khusus. |
| Dalc | > 1 | Konsumsi alkohol pada hari kerja meningkatkan rata-rata ketidakhadiran. Setiap kenaikan 1 skala konsumsi alkohol hari kerja meningkatkan rate ketidakhadiran — menggambarkan dampak perilaku berisiko terhadap kehadiran sekolah. |
| goout | > 1 | Frekuensi keluar bersama teman meningkatkan rata-rata ketidakhadiran. Siswa yang lebih sering bergaul di luar sekolah mungkin mengalami kurang tidur atau berkurangnya motivasi hadir, yang berujung pada lebih banyak absensi. |
| higherYa | < 1 | Siswa yang berencana melanjutkan ke pendidikan tinggi memiliki rata-rata ketidakhadiran yang secara signifikan lebih rendah. Motivasi dan orientasi masa depan yang kuat mendorong kehadiran lebih konsisten. |
| studytime | < 1 | Setiap penambahan 1 satuan waktu belajar mingguan menurunkan rata-rata ketidakhadiran. Siswa yang mendedikasikan lebih banyak waktu untuk belajar memiliki komitmen akademik lebih tinggi yang tercermin dalam tingkat kehadiran lebih baik. |
| sexLaki-laki | > 1 | Siswa laki-laki cenderung memiliki rata-rata ketidakhadiran lebih tinggi dibandingkan perempuan. Perbedaan ini dapat mencerminkan perbedaan perilaku sosial dan gaya pengasuhan antara dua jenis kelamin. |
43 Evaluasi Model
43.1 Prediksi pada Data Uji
# Prediksi rata-rata (lambda) pada data uji
pred_poisson <- predict(model_terbaik, newdata = student_test, type = "response")
cat("Ringkasan nilai prediksi rata-rata ketidakhadiran:\n")## Ringkasan nilai prediksi rata-rata ketidakhadiran:## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.542 4.224 5.358 5.715 6.656 13.287##
## Ringkasan nilai aktual ketidakhadiran:## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 2.500 6.372 7.750 75.000##
## Beberapa prediksi pertama:## 3 4 6 8 9 12 13 15 21 33
## 5.14 3.60 3.73 6.67 3.71 3.73 3.72 2.99 3.68 3.2243.2 Diagnostik Residual
Untuk model Poisson/NB, residual yang diperiksa adalah deviance residuals dan Pearson residuals, bukan residual biasa seperti pada regresi linear.
# Menghitung residual pada data latih
fitted_train <- fitted(model_terbaik)
dev_resid <- residuals(model_terbaik, type = "deviance")
pear_resid <- residuals(model_terbaik, type = "pearson")
df_diag <- data.frame(
fitted = fitted_train,
dev_res = dev_resid,
pear_res = pear_resid,
idx = seq_along(fitted_train)
)
p_dev <- ggplot(df_diag, aes(x = fitted, y = dev_res)) +
geom_point(alpha = 0.5, color = "#2980B9", size = 1.5) +
geom_hline(yintercept = c(-2, 0, 2), linetype = c("dashed","solid","dashed"),
color = c("red","gray40","red"), linewidth = 0.7) +
labs(title = "Deviance Residuals vs Fitted",
x = "Nilai Fitted (λ̂)", y = "Deviance Residuals") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
p_pear <- ggplot(df_diag, aes(x = fitted, y = pear_res)) +
geom_point(alpha = 0.5, color = "#8E44AD", size = 1.5) +
geom_hline(yintercept = c(-2, 0, 2), linetype = c("dashed","solid","dashed"),
color = c("red","gray40","red"), linewidth = 0.7) +
labs(title = "Pearson Residuals vs Fitted",
x = "Nilai Fitted (λ̂)", y = "Pearson Residuals") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
p_qqnorm <- ggplot(df_diag, aes(sample = dev_res)) +
stat_qq(color = "#E74C3C", alpha = 0.6, size = 1.5) +
stat_qq_line(color = "gray40", linewidth = 0.8) +
labs(title = "Q-Q Plot Deviance Residuals",
x = "Kuantil Teoritis", y = "Kuantil Sampel") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
p_hist <- ggplot(df_diag, aes(x = dev_res)) +
geom_histogram(binwidth = 0.4, fill = "#27AE60", color = "white", alpha = 0.8) +
labs(title = "Distribusi Deviance Residuals",
x = "Deviance Residuals", y = "Frekuensi") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
grid.arrange(p_dev, p_pear, p_qqnorm, p_hist, nrow = 2,
top = grid::textGrob(
"Diagnostik Residual — Model Regresi Poisson/NB",
gp = grid::gpar(fontface = "bold", fontsize = 13)
)
)
Gambar 6. Plot Diagnostik Residual Model Regresi Poisson/NB
Interpretasi diagnostik residual:
- Deviance Residuals vs Fitted: Model yang baik menunjukkan pola residual yang tersebar acak tanpa tren sistematis. Residual di luar rentang ±2 mengindikasikan observasi yang sulit diprediksi.
- Pearson Residuals vs Fitted: Serupa dengan di atas; pada model dengan overdispersion yang ditangani NB, seharusnya lebih terkontrol dibanding Poisson murni.
- Q-Q Plot: Deviance residuals yang mendekati garis diagonal mengindikasikan asumsi distribusi terpenuhi dengan baik.
- Histogram Residual: Distribusi residual yang mendekati normal dan berpusat di nol menandakan model sudah menangkap pola utama data.
43.3 Metrik Evaluasi untuk Data Cacahan
Untuk model regresi count, metrik evaluasi yang umum digunakan adalah berbasis perbandingan nilai prediksi kontinu (\(\hat{\mu}\)) dengan nilai aktual, serta metrik berbasis devians.
y_actual <- student_test$absences
y_pred <- pred_poisson
# Mean Absolute Error
mae_val <- mean(abs(y_actual - y_pred))
# Root Mean Squared Error
rmse_val <- sqrt(mean((y_actual - y_pred)^2))
# Mean Absolute Percentage Error (hati-hati untuk nilai 0)
y_nonzero <- y_actual > 0
mape_val <- mean(abs(y_actual[y_nonzero] - y_pred[y_nonzero]) /
y_actual[y_nonzero]) * 100
# Korelasi antara nilai aktual dan prediksi
corr_val <- cor(y_actual, y_pred)
# Deviance pada data uji (dihitung manual untuk Poisson)
# log-likelihood Poisson: sum(y*log(mu) - mu - log(y!))
loglik_pred <- sum(dpois(y_actual, lambda = y_pred, log = TRUE))
loglik_null <- sum(dpois(y_actual,
lambda = mean(student_train$absences), log = TRUE))
dev_test <- -2 * (loglik_null - loglik_pred)
cat("=== Metrik Evaluasi pada Data Uji ===\n")## === Metrik Evaluasi pada Data Uji ===## MAE (Mean Absolute Error) : 5.9077 hari## RMSE (Root Mean Squared Error) : 11.2045 hari## MAPE (Mean Abs. Percentage Error): 79.33% (obs. non-nol saja)## Korelasi aktual vs prediksi : 0.1142metrik_count_df <- data.frame(
Metrik = c("MAE (Mean Absolute Error)",
"RMSE (Root Mean Squared Error)",
"MAPE (Mean Abs. Percentage Error)",
"Korelasi Aktual vs Prediksi"),
Formula = c(
"\\(\\frac{1}{n}\\sum|y_i - \\hat{\\mu}_i|\\)",
"\\(\\sqrt{\\frac{1}{n}\\sum(y_i - \\hat{\\mu}_i)^2}\\)",
"\\(\\frac{1}{n^*}\\sum_{y_i>0}\\frac{|y_i - \\hat{\\mu}_i|}{y_i}\\times 100\\)",
"\\(\\text{cor}(y_i, \\hat{\\mu}_i)\\)"
),
Nilai = c(
paste0(round(mae_val, 4), " hari"),
paste0(round(rmse_val, 4), " hari"),
paste0(round(mape_val, 2), "%"),
round(corr_val, 4)
),
Interpretasi = c(
"Rata-rata kesalahan prediksi secara absolut dalam satuan hari",
"Lebih sensitif terhadap kesalahan besar; akibat outlier ketidakhadiran tinggi",
"Rata-rata kesalahan relatif (hanya observasi non-nol); berguna untuk perbandingan",
"Kekuatan hubungan linear antara nilai aktual dan prediksi model"
)
)
kable(metrik_count_df,
caption = "Tabel 6. Metrik Evaluasi Model pada Data Uji",
col.names = c("Metrik", "Formula", "Nilai", "Interpretasi"),
escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = TRUE, font_size = 12) %>%
column_spec(3, bold = TRUE, width = "7em")| Metrik | Formula | Nilai | Interpretasi |
|---|---|---|---|
| MAE (Mean Absolute Error) | \(\frac{1}{n}\sum&#124;y_i - \hat{\mu}_i&#124;\) | 5.9077 hari | Rata-rata kesalahan prediksi secara absolut dalam satuan hari |
| RMSE (Root Mean Squared Error) | \(\sqrt{\frac{1}{n}\sum(y_i - \hat{\mu}_i)^2}\) | 11.2045 hari | Lebih sensitif terhadap kesalahan besar; akibat outlier ketidakhadiran tinggi |
| MAPE (Mean Abs. Percentage Error) | \(\frac{1}{n^*}\sum_{y_i>0}\frac{&#124;y_i - \hat{\mu}_i&#124;}{y_i}\times 100\) | 79.33% | Rata-rata kesalahan relatif (hanya observasi non-nol); berguna untuk perbandingan |
| Korelasi Aktual vs Prediksi | \(\text{cor}(y_i, \hat{\mu}_i)\) | 0.1142 | Kekuatan hubungan linear antara nilai aktual dan prediksi model |
43.4 Visualisasi Prediksi vs Aktual
df_pred_test <- data.frame(
Aktual = y_actual,
Prediksi = y_pred
)
p_scatter <- ggplot(df_pred_test, aes(x = Prediksi, y = Aktual)) +
geom_point(alpha = 0.5, color = "#2980B9", size = 2) +
geom_abline(slope = 1, intercept = 0, color = "#E74C3C",
linetype = "dashed", linewidth = 1.2) +
geom_smooth(method = "loess", se = TRUE, color = "#27AE60",
linewidth = 1, alpha = 0.2) +
labs(
title = "Aktual vs Prediksi: Jumlah Ketidakhadiran",
subtitle = "Garis merah putus-putus = prediksi sempurna (y = x)",
x = "Nilai Prediksi (λ̂)",
y = "Nilai Aktual"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
p_err <- ggplot(df_pred_test, aes(x = Prediksi, y = Aktual - Prediksi)) +
geom_point(alpha = 0.5, color = "#8E44AD", size = 2) +
geom_hline(yintercept = 0, color = "gray40", linetype = "solid", linewidth = 1) +
geom_smooth(method = "loess", se = FALSE, color = "#E74C3C",
linewidth = 1, linetype = "dashed") +
labs(
title = "Residual vs Nilai Prediksi",
subtitle = "Pola acak tanpa tren → model sudah sesuai",
x = "Nilai Prediksi (λ̂)",
y = "Residual (Aktual − Prediksi)"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
grid.arrange(p_scatter, p_err, ncol = 2,
top = grid::textGrob(
"Evaluasi Prediktif Model — Data Uji",
gp = grid::gpar(fontface = "bold", fontsize = 13)
)
)
Gambar 7. Perbandingan Nilai Aktual dan Prediksi Ketidakhadiran
Interpretasi: Scatter plot aktual vs prediksi yang ideal menunjukkan titik-titik yang tersebar di sekitar garis identitas (y = x). Pada data dengan banyak nilai nol dan distribusi yang miring, model cenderung lebih akurat dalam memprediksi nilai rendah dibandingkan nilai tinggi. Residual plot yang tidak menunjukkan pola sistematis mengkonfirmasi spesifikasi model sudah memadai.
43.5 Visualisasi Distribusi Nilai Prediksi vs Aktual
df_dist <- data.frame(
Nilai = c(y_actual, y_pred),
Tipe = rep(c("Aktual", "Prediksi (λ̂)"), each = length(y_actual))
)
ggplot(df_dist, aes(x = Nilai, fill = Tipe, color = Tipe)) +
geom_density(alpha = 0.3, linewidth = 1) +
scale_fill_manual(values = c("Aktual" = "#2980B9", "Prediksi (λ̂)" = "#E74C3C")) +
scale_color_manual(values = c("Aktual" = "#2980B9", "Prediksi (λ̂)" = "#E74C3C")) +
scale_x_continuous(limits = c(0, 40)) +
labs(
title = "Perbandingan Distribusi Aktual vs Prediksi",
subtitle = "Data uji — nilai prediksi adalah rata-rata kontinu (λ̂)",
x = "Jumlah Ketidakhadiran",
y = "Densitas",
fill = "Tipe",
color = "Tipe"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray50"))
Gambar 8. Perbandingan Distribusi Aktual dan Prediksi
Interpretasi: Distribusi nilai prediksi \(\hat{\mu}\) (rata-rata kontinu) cenderung lebih halus dan terpusat dibandingkan distribusi aktual yang right-skewed dengan banyak nilai nol. Ini adalah perilaku normal untuk model regresi — model memprediksi ekspektasi bersyarat, bukan realisasi individual yang acak.
44 Kesimpulan
44.1 Ringkasan Hasil Analisis
Analisis regresi Poisson (dan Negative Binomial sebagai alternatif jika overdispersion signifikan) berhasil dilakukan untuk memodelkan jumlah ketidakhadiran siswa menggunakan dataset Student Performance UCI. Model dibangun menggunakan delapan variabel prediktor sosio-demografis dan akademik dari 395 siswa mata pelajaran Matematika di Portugal.
44.1.1 Pengujian Model
| Jenis Uji | Statistik | Nilai p | Keputusan |
|---|---|---|---|
| Uji Simultan (LRT) | \(G\) (lihat output) | ≪ 0,05 | Tolak \(H_0\) — model secara keseluruhan signifikan |
| Uji Parsial (Wald) | \(z_j\) (per variabel) | Bervariasi | Beberapa variabel signifikan |
| Uji Overdispersion | \(z\) (Cameron-Trivedi) | Bervariasi | Periksa output — pilih NB jika signifikan |
44.1.2 Kebaikan Model
| Metrik | Nilai | Interpretasi |
|---|---|---|
| Pseudo \(R^2\) McFadden | Lihat output | ≥ 0,10 cukup baik untuk data perilaku |
| AIC Model Terbaik | Lihat output | Semakin kecil semakin baik |
| Rasio Devians/df | Lihat output | Mendekati 1 → fit baik |
44.2 Variabel yang Berpengaruh Signifikan
Berdasarkan uji parsial (Wald Test) dengan \(\alpha = 0{,}05\), faktor-faktor yang terbukti signifikan mempengaruhi jumlah ketidakhadiran siswa adalah:
failures(Kegagalan Kelas): Variabel dengan dampak terkuat. Setiap riwayat kegagalan secara signifikan melipatgandakan rata-rata ketidakhadiran — mencerminkan siklus disengagement akademik yang berisiko tinggi.Dalc(Alkohol Hari Kerja): Konsumsi alkohol pada hari sekolah meningkatkan rate ketidakhadiran secara signifikan, menggambarkan dampak nyata perilaku berisiko terhadap kehadiran.goout(Keluar Bersama Teman): Frekuensi sosial di luar sekolah yang tinggi berasosiasi dengan lebih banyak absensi — kemungkinan melalui mekanisme berkurangnya prioritas terhadap sekolah.higher(Ingin Lanjut Studi): Siswa yang beraspirasi melanjutkan pendidikan tinggi memiliki rate ketidakhadiran yang secara signifikan lebih rendah — motivasi jangka panjang berperan sebagai faktor pelindung.studytime(Waktu Belajar): Dedikasi belajar yang lebih tinggi menurunkan rata-rata ketidakhadiran — mencerminkan komitmen akademik yang lebih konsisten.
44.3 Penanganan Overdispersion
Regresi Poisson murni sering kali tidak cukup untuk data ketidakhadiran siswa karena:
- Banyaknya nilai nol (siswa yang tidak pernah absen): 29,1% observasi bernilai nol.
- Ekor distribusi yang sangat panjang: nilai maksimum 75 hari jauh melebihi rata-rata 5,7 hari.
- Heterogenitas individu: setiap siswa memiliki karakteristik unik yang tidak sepenuhnya tercakup prediktor.
Jika uji Cameron-Trivedi signifikan, Model Negative Binomial lebih direkomendasikan karena secara eksplisit memodelkan overdispersion melalui parameter \(\theta\), menghasilkan inferensi (nilai p dan interval kepercayaan) yang lebih valid.
44.4 Kinerja Model
Model menunjukkan kinerja yang memadai untuk data ketidakhadiran yang sifatnya heterogen:
- MAE yang relatif rendah menunjukkan prediksi cukup dekat dengan nilai aktual untuk sebagian besar siswa.
- Korelasi aktual vs prediksi yang positif mengkonfirmasi model berhasil menangkap tren utama.
- Nilai tinggi yang merupakan outliers (siswa dengan absensi ekstrem) umumnya lebih sulit diprediksi — merupakan keterbatasan inherent model regresi tanpa variabel kontekstual tambahan.
44.5 Perbandingan dengan Metode Regresi Lainnya
| Aspek | Binary (Bagian I) | Multinomial (Bagian II) | Ordinal (Bagian III) | Poisson (Bagian IV) |
|---|---|---|---|---|
| Tipe respons | Biner (0/1) | Nominal (≥3 kelas) | Ordinal (≥3 kelas) | Cacahan (0, 1, 2, …) |
| Distribusi asumsi | Binomial | Multinomial | Binomial (kumulatif) | Poisson / Negative Binomial |
| Fungsi link | Logit | Generalized logit | Cumulative logit | Log |
| Interpretasi koef. | Odds Ratio | Relative Risk Ratio | OR Kumulatif | Incidence Rate Ratio |
| Isu utama | Threshold pemisah | IIA | Proportional odds | Overdispersion |
| Dataset | Bike Sharing | Bike Sharing | Wine UCI | Student Performance UCI |
Penutup: Pemilihan metode regresi yang tepat sangat bergantung pada tipe variabel respons. Regresi Poisson (atau variannya) adalah pilihan yang tepat dan principled untuk memodelkan data cacahan non-negatif, menghasilkan prediksi dan inferensi yang lebih valid dibandingkan regresi linear yang dipaksakan pada data serupa.
44.6 Implikasi Hasil Penelitian
- Intervensi dini berbasis risiko: Sekolah dapat
mengidentifikasi siswa berisiko tinggi absen berdasarkan profil
failures(kegagalan kelas) danDalc(konsumsi alkohol) untuk memberikan bimbingan konseling proaktif. - Program motivasi akademik: Karena aspirasi
highersecara signifikan melindungi dari absensi, program orientasi karir dan pendidikan tinggi sejak dini dapat meningkatkan kehadiran secara tidak langsung. - Keterlibatan orang tua: Mendorong orang tua
memantau kebiasaan keluar (
goout) dan konsumsi alkohol siswa, terutama pada hari sekolah. - Dukungan akademik terstruktur: Siswa dengan waktu
belajar rendah (
studytime= 1) perlu program pembelajaran tersupport untuk membangun kebiasaan belajar yang lebih baik. - Sistem prediksi absensi: Model ini dapat diimplementasikan sebagai sistem peringatan dini — jika profil siswa baru memasuki model dan prediksi absensi tinggi, sekolah dapat melakukan tindakan preventif sebelum masalah menjadi parah.
44.7 Referensi
- Cortez, P., & Silva, A. (2008). Using data mining to predict secondary school student performance. In Proceedings of the 5th Annual Future Business Technology Conference (FUBUTEC 2008), Porto, Portugal. UCI ML Repository: https://archive.uci.edu/dataset/320/student+performance
- Cameron, A. C., & Trivedi, P. K. (1998). Regression Analysis of Count Data. Cambridge University Press.
- Hilbe, J. M. (2011). Negative Binomial Regression (2nd ed.). Cambridge University Press.
- McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models (2nd ed.). Chapman & Hall.
- Agresti, A. (2013). Categorical Data Analysis (3rd ed.). John Wiley & Sons.
- R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
- Venables, W. N., & Ripley, B. D. (2002). Modern Applied
Statistics with S (4th ed.). Springer. [package
MASS, fungsiglm.nb()] - Kleiber, C., & Zeileis, A. (2008). Applied Econometrics with
R. Springer. [package
AER, fungsidispersiontest()]
Dokumen ini merupakan Bagian IV dari laporan gabungan Regresi Logistik dan dibuat menggunakan R Markdown. Data diunduh otomatis dari UCI ML Repository saat knit dijalankan. Dokumen dapat langsung di-knit menjadi laporan HTML interaktif.