Analisis Data Kategori
Chelsea Violeta Hernathan
01 June 2026
1 Definisi Analisis Data Kategori
1.1 Pengertian Analisis Data Kategori
Analisis data kategori merupakan cabang analisis statistik yang digunakan untuk menganalisis data dengan variabel yang berbentuk kategori atau kelas tertentu. Variabel kategori adalah variabel yang memiliki skala pengukuran berupa sekumpulan kategori yang merepresentasikan karakteristik objek yang diamati.
Dalam analisis data kategori, metode statistik yang digunakan umumnya berfokus pada frekuensi, proporsi, dan peluang kemunculan suatu kategori dalam populasi.
Selain itu, analisis ini juga digunakan untuk mempelajari hubungan antar variabel kategori melalui pendekatan probabilistik dan model statistik yang sesuai.
Dengan demikian, analisis data kategori menjadi penting karena banyak fenomena dalam kehidupan nyata tidak dinyatakan dalam bentuk numerik kontinu, melainkan dalam bentuk klasifikasi atau kelompok.
1.2 Karakteristik Variabel Kategori
Variabel kategori memiliki beberapa karakteristik utama sebagai berikut:
Berbentuk kategori atau kelas
Nilai variabel dinyatakan dalam bentuk kategori yang menggambarkan karakteristik objek, seperti jenis kelamin, jenis pekerjaan, atau jenis tempat tinggal.Tidak memiliki makna numerik secara langsung
Nilai kategori tidak dapat diinterpretasikan sebagai besaran numerik sehingga operasi aritmetika seperti penjumlahan atau rata-rata tidak relevan.Dapat berupa nominal atau ordinal
Variabel kategori dapat bersifat:- Nominal → tidak memiliki urutan antar
kategori
- Ordinal → memiliki urutan atau tingkatan antar kategori
- Nominal → tidak memiliki urutan antar
kategori
Dianalisis menggunakan frekuensi atau proporsi
Analisis dilakukan berdasarkan jumlah kemunculan (frekuensi) atau proporsi pada setiap kategori.Disajikan dalam tabel distribusi atau tabel kontingensi
Data kategori umumnya ditampilkan dalam bentuk tabel untuk memahami distribusi dan hubungan antar variabel.
1.3 Contoh Penerapan Analisis Data Kategori dalam Penelitian
Analisis data kategori banyak digunakan dalam berbagai bidang penelitian karena banyak fenomena yang secara alami berbentuk kategori.
Penelitian kesehatan
Dalam bidang kesehatan, analisis data kategori digunakan untuk menganalisis hubungan antara status merokok (perokok dan bukan perokok) dengan kejadian penyakit tertentu. Analisis dilakukan menggunakan tabel kontingensi dan uji Chi-Square untuk mengetahui apakah terdapat hubungan yang signifikan antara kedua variabel tersebut.Penelitian sosial
Dalam penelitian sosial, analisis data kategori digunakan untuk mempelajari hubungan antara tingkat pendidikan dengan status pekerjaan. Hasil analisis dapat memberikan gambaran mengenai pola sosial dalam masyarakat.Penelitian pemasaran
Dalam bidang pemasaran, analisis data kategori digunakan untuk mengetahui hubungan antara jenis kelamin konsumen dengan preferensi terhadap suatu produk.
Analisis ini penting dalam memahami segmentasi pasar dan perilaku konsumen sehingga dapat mendukung pengambilan keputusan bisnis.
2 Tabel Kontingensi 2×2
2.1 Tabel Kontingensi 2×2
2.1.1 Definisi Tabel Kontingensi
Tabel kontingensi merupakan tabel yang digunakan untuk menyajikan data kategorik dalam bentuk frekuensi yang menunjukkan hubungan antara dua atau lebih variabel kategori. Tabel ini biasanya digunakan untuk menggambarkan bagaimana distribusi satu variabel berhubungan dengan distribusi variabel lainnya.
Menurut Agresti (2013), tabel kontingensi merupakan alat dasar dalam analisis data kategori yang menampilkan frekuensi observasi dari kombinasi kategori antara dua atau lebih variabel. Melalui tabel ini, peneliti dapat mengamati pola hubungan antar variabel serta menjadi dasar dalam berbagai analisis statistik seperti uji Chi-Square, analisis asosiasi, dan model log-linear.
2.1.2 Struktur Tabel Kontingensi
Secara umum, tabel kontingensi terdiri dari baris dan kolom yang merepresentasikan kategori dari dua variabel yang berbeda. Setiap sel pada tabel menunjukkan jumlah observasi yang termasuk dalam kombinasi kategori tertentu.
Sebagai contoh, berikut adalah tabel kontingensi 2 × 2 yang menunjukkan hubungan antara keikutsertaan dalam bimbingan belajar dan kelulusan ujian.
| Bimbingan Belajar | Lulus Ujian | Tidak Lulus | Total |
|---|---|---|---|
| Mengikuti Bimbel | 80 | 20 | 100 |
| Tidak Mengikuti | 60 | 40 | 100 |
| Total | 140 | 60 | 200 |
Pada tabel tersebut:
- Baris menunjukkan kategori dari
variabel keikutsertaan dalam bimbingan belajar
- Kolom menunjukkan kategori dari
variabel hasil ujian
- Sel tabel menunjukkan jumlah observasi pada kombinasi kedua variabel
Sebagai contoh, nilai 80 menunjukkan bahwa terdapat 80 mahasiswa yang mengikuti bimbingan belajar dan lulus ujian.
2.1.3 Konsep Joint Distribution
Joint distribution atau distribusi gabungan menggambarkan probabilitas terjadinya dua kategori secara bersamaan. Dalam tabel kontingensi, distribusi gabungan diperoleh dengan membagi frekuensi pada setiap sel dengan jumlah total observasi.
Sebagai contoh, probabilitas seorang mahasiswa mengikuti bimbingan belajar dan lulus ujian adalah:
\[ P(\text{Bimbel dan Lulus}) = \frac{80}{200} = 0.40 \]
Nilai tersebut menunjukkan bahwa 40% dari seluruh mahasiswa dalam sampel mengikuti bimbingan belajar dan berhasil lulus ujian.
Contoh lain:
\[ P(\text{Tidak Bimbel dan Lulus}) = \frac{60}{200} = 0.30 \]
Distribusi gabungan memberikan informasi mengenai peluang kombinasi dua kategori yang terjadi secara simultan dalam populasi atau sampel penelitian.
2.1.4 Konsep Marginal Distribution
Distribusi marginal merupakan distribusi probabilitas dari satu variabel tanpa mempertimbangkan variabel lainnya.
Sebagai contoh:
- Probabilitas mahasiswa mengikuti bimbingan belajar
\[ P(\text{Bimbel}) = \frac{100}{200} = 0.50 \]
- Probabilitas mahasiswa lulus ujian
\[ P(\text{Lulus}) = \frac{140}{200} = 0.70 \]
Distribusi marginal memberikan gambaran distribusi masing-masing variabel secara terpisah.
2.1.5 Konsep Conditional Probability
Probabilitas kondisional merupakan probabilitas suatu kejadian dengan syarat kejadian lain telah terjadi.
Sebagai contoh:
\[ P(\text{Lulus | Bimbel}) = \frac{80}{100} = 0.80 \]
Artinya, 80% mahasiswa yang mengikuti bimbingan belajar berhasil lulus ujian.
Sebaliknya:
\[ P(\text{Lulus | Tidak Bimbel}) = \frac{60}{100} = 0.60 \]
Perbandingan ini menunjukkan bahwa mahasiswa yang mengikuti bimbingan belajar memiliki peluang lebih tinggi untuk lulus dibandingkan yang tidak.
Konsep ini menjadi dasar penting dalam memahami hubungan atau asosiasi antar variabel kategorik.
2.2 Ukuran Asosiasi pada Tabel Kontingensi
Ukuran asosiasi digunakan untuk mengukur kekuatan hubungan antara dua variabel kategorik dalam tabel kontingensi. Pada tabel kontingensi 2 × 2, ukuran asosiasi yang umum digunakan adalah Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR).
Sebagai ilustrasi, digunakan tabel kontingensi 2 × 2 berikut yang menggambarkan hubungan antara keikutsertaan dalam bimbingan belajar dan kelulusan ujian.
| Lulus Ujian | Tidak Lulus | |
|---|---|---|
| Mengikuti Bimbel | a | c |
| Tidak Mengikuti | b | d |
Total observasi:
\[ n = a + b + c + d \]
Pada tabel tersebut:
- a → mengikuti
bimbingan belajar dan lulus ujian
- b → mengikuti
bimbingan belajar tetapi tidak lulus ujian
- c → tidak mengikuti
bimbingan belajar tetapi lulus ujian
- d → tidak mengikuti bimbingan belajar dan tidak lulus ujian
2.2.1 Risk Difference
Risk Difference (RD) merupakan ukuran asosiasi yang menyatakan perbedaan probabilitas kejadian antara dua kelompok.
Secara matematis:
\[ RD = P(\text{event | bimbel}) - P(\text{event | tidak bimbel}) \]
Dalam tabel kontingensi 2 × 2:
\[ RD = \frac{a}{a+b} - \frac{c}{c+d} \]
Interpretasi:
- RD = 0 → tidak ada perbedaan
- RD > 0 → peluang lebih tinggi pada kelompok
bimbel
- RD < 0 → peluang lebih rendah pada kelompok bimbel
Contoh: Jika RD = 0.20, maka peluang lulus meningkat sebesar 20% pada mahasiswa yang mengikuti bimbingan belajar.
2.2.2 Relative Risk
Relative Risk (RR) membandingkan rasio probabilitas kejadian antara dua kelompok.
\[ RR = \frac{P(\text{event | bimbel})}{P(\text{event | tidak bimbel})} \]
Dalam tabel:
\[ RR = \frac{a/(a+b)}{c/(c+d)} \]
Interpretasi:
- RR = 1 → tidak ada hubungan
- RR > 1 → peluang lebih besar pada kelompok
bimbel
- RR < 1 → peluang lebih kecil pada kelompok bimbel
Contoh: Jika RR = 2, maka peluang lulus dua kali lebih besar pada mahasiswa yang mengikuti bimbingan belajar.
2.2.3 Odds Ratio
Odds Ratio (OR) membandingkan odds kejadian antara dua kelompok.
\[ OR = \frac{a/b}{c/d} \]
Disederhanakan menjadi:
\[ OR = \frac{ad}{bc} \]
Interpretasi:
- OR = 1 → tidak ada hubungan
- OR > 1 → kejadian lebih mungkin terjadi pada
kelompok bimbel
- OR < 1 → kejadian lebih kecil kemungkinannya
Contoh: Jika OR = 3, maka peluang lulus tiga kali lebih besar pada mahasiswa yang mengikuti bimbingan belajar.
Ketiga ukuran asosiasi ini memberikan perspektif yang berbeda dalam memahami hubungan antar variabel, sehingga pemilihannya harus disesuaikan dengan tujuan analisis.
2.3 Contoh Perhitungan Manual
Untuk memahami konsep ukuran asosiasi pada tabel kontingensi, berikut diberikan contoh kasus sederhana mengenai hubungan antara keikutsertaan dalam bimbingan belajar dan kelulusan ujian.
Misalkan dilakukan pengamatan terhadap 200 orang mahasiswa dan diperoleh data sebagai berikut.
2.3.1 Tabel Kontingensi
| Keikutsertaan Bimbel | Lulus Ujian | Tidak Lulus | Total |
|---|---|---|---|
| Mengikuti Bimbel | 80 | 20 | 100 |
| Tidak Mengikuti | 60 | 40 | 100 |
| Total | 140 | 60 | 200 |
Pada tabel tersebut dapat didefinisikan:
\[ a = 80, \quad b = 20, \quad c = 60, \quad d = 40 \]
Total pengamatan:
\[ n = a + b + c + d = 200 \]
2.3.2 Menghitung Peluang Bersyarat
Peluang bersyarat menunjukkan probabilitas suatu kejadian dengan syarat kejadian lain telah terjadi.
\[ P(\text{Lulus | Mengikuti Bimbel}) = \frac{80}{100} = 0.80 \]
Artinya, 80% mahasiswa yang mengikuti bimbingan belajar lulus ujian.
\[ P(\text{Lulus | Tidak Mengikuti Bimbel}) = \frac{60}{100} = 0.60 \]
Artinya, 60% mahasiswa yang tidak mengikuti bimbingan belajar lulus ujian.
Perbandingan ini menunjukkan bahwa mahasiswa yang mengikuti bimbingan belajar memiliki peluang lebih besar untuk lulus.
2.3.3 Menghitung Odds
Odds merupakan rasio antara probabilitas kejadian dan tidak kejadian.
\[ Odds_{bimbel} = \frac{80}{20} = 4.00 \]
Artinya, terdapat 4 kali peluang lulus dibanding tidak lulus pada kelompok bimbel.
\[ Odds_{tidak\ bimbel} = \frac{60}{40} = 1.50 \]
Artinya, terdapat 1.5 kali peluang lulus dibanding tidak lulus pada kelompok non-bimbel.
2.3.4 Menghitung Odds Ratio
Odds Ratio (OR) digunakan untuk membandingkan odds antara dua kelompok.
\[ OR = \frac{ad}{bc} \]
Substitusi:
\[ OR = \frac{(80)(40)}{(20)(60)} = \frac{3200}{1200} = 2.67 \]
Interpretasi:
Odds mahasiswa yang mengikuti bimbingan belajar untuk lulus sekitar 2.67 kali lebih besar dibandingkan yang tidak.
Hasil ini menunjukkan adanya hubungan antara keikutsertaan bimbingan belajar dan kelulusan ujian.
2.4 Analisis Menggunakan R
Untuk melengkapi perhitungan manual yang telah dilakukan sebelumnya, analisis yang sama dapat dilakukan menggunakan perangkat lunak R.
Analisis ini meliputi: - pembuatan tabel kontingensi - perhitungan Odds dan Odds Ratio - pengujian hubungan menggunakan uji Chi-Square
2.4.1 Membuat Tabel Kontingensi
Langkah pertama adalah membentuk tabel kontingensi menggunakan fungsi
matrix() di R.
data <- matrix(c(80,20,60,40),
nrow = 2,
byrow = TRUE)
rownames(data) <- c("Mengikuti Bimbel","Tidak Mengikuti")
colnames(data) <- c("Lulus Ujian","Tidak Lulus")
data## Lulus Ujian Tidak Lulus
## Mengikuti Bimbel 80 20
## Tidak Mengikuti 60 402.4.2 Menghitung Odds dan Odds Ratio
Odds merupakan perbandingan antara peluang terjadinya suatu kejadian dengan peluang tidak terjadinya kejadian tersebut.
a <- data[1,1]
b <- data[1,2]
c <- data[2,1]
d <- data[2,2]
Odds_bimbel <- a/b
Odds_tidak_bimbel <- c/d
Odds_bimbel## [1] 4## [1] 1.5Odds pada kelompok yang mengikuti bimbingan belajar lebih tinggi dibandingkan kelompok yang tidak mengikuti.
Odds Ratio digunakan untuk membandingkan odds antara dua kelompok.
## [1] 2.666667Nilai Odds Ratio menunjukkan kekuatan hubungan antara bimbingan belajar dan kelulusan ujian.
2.4.3 Uji Chi-Square
Uji Chi-Square digunakan untuk menguji apakah terdapat hubungan yang signifikan antara kedua variabel.
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: data
## X-squared = 8.5952, df = 1, p-value = 0.00337Jika p-value < 0.05, maka terdapat hubungan yang signifikan antara keikutsertaan bimbingan belajar dan kelulusan ujian.
2.4.4 Visualisasi Data
Visualisasi digunakan untuk memperjelas pola hubungan antar variabel.
mosaicplot(data,
main = "Hubungan Bimbingan Belajar dan Kelulusan Ujian",
col = c("lightblue","pink"))
Plot ini menunjukkan distribusi proporsi masing-masing kategori secara visual, sehingga memudahkan dalam memahami hubungan antara kedua variabel.
2.5 Interpretasi Hasil
2.5.1 Interpretasi Statistik
Berdasarkan hasil analisis menggunakan R, diperoleh nilai Odds Ratio sebesar
\[ OR = 2.67 \]
Nilai tersebut menunjukkan bahwa odds mahasiswa lulus ujian pada kelompok yang mengikuti bimbingan belajar sekitar 2.67 kali lebih besar dibandingkan dengan mahasiswa yang tidak mengikuti bimbingan belajar.
Selain itu, hasil uji Chi-Square menunjukkan nilai statistik
\[ X^2 = 8.595 \]
dengan derajat kebebasan
\[ df = 1 \]
serta nilai
\[ p\text{-value} = 0.00337 \]
Karena p-value < 0.05, maka hipotesis nol \(H_0\) yang menyatakan bahwa tidak terdapat hubungan antara keikutsertaan dalam bimbingan belajar dan kelulusan ujian ditolak.
Dengan demikian, terdapat hubungan yang signifikan secara statistik antara keikutsertaan dalam bimbingan belajar dan keberhasilan mahasiswa dalam lulus ujian.
2.5.2 Interpretasi Substantif
Dalam konteks kasus ini, hasil analisis menunjukkan bahwa mahasiswa yang mengikuti bimbingan belajar memiliki kemungkinan yang lebih tinggi untuk lulus ujian dibandingkan mahasiswa yang tidak mengikuti bimbingan belajar.
Perbedaan probabilitas kelulusan ujian antara kedua kelompok dapat dilihat dari peluang bersyarat berikut:
\[ P(Lulus\ Ujian \mid Mengikuti\ Bimbel) = 0.8 \]
\[ P(Lulus\ Ujian \mid Tidak\ Mengikuti\ Bimbel) = 0.6 \]
Nilai tersebut menunjukkan bahwa proporsi mahasiswa yang lulus ujian pada kelompok yang mengikuti bimbingan belajar lebih tinggi dibandingkan pada kelompok mahasiswa yang tidak mengikuti bimbingan belajar.
Secara substantif, keikutsertaan dalam bimbingan belajar dapat menjadi faktor yang berperan dalam meningkatkan peluang mahasiswa untuk lulus ujian.
Oleh karena itu, penyediaan fasilitas bimbingan belajar yang efektif dapat menjadi salah satu strategi yang dapat membantu meningkatkan tingkat kelulusan mahasiswa dalam suatu ujian.
3 Inferensi pada Tabel Kontingensi Dua Arah
Inferensi statistik merupakan proses penarikan kesimpulan mengenai populasi berdasarkan data sampel. Dalam konteks tabel kontingensi dua arah, inferensi digunakan untuk menganalisis hubungan antara dua variabel kategorik yang disusun dalam bentuk tabel silang.
Tabel kontingensi adalah tabel yang menyajikan distribusi frekuensi gabungan dari dua variabel kategorik dalam bentuk matriks. Melalui tabel ini, hubungan antar variabel dapat dipahami melalui pendekatan estimasi maupun pengujian hipotesis.
Secara umum, inferensi dalam tabel kontingensi dua arah terdiri atas dua bagian utama, yaitu:
- Estimasi
- Pengujian
3.1 Estimasi
Estimasi bertujuan untuk memperkirakan nilai parameter populasi berdasarkan data sampel yang tersedia. Pendekatan ini menjadi langkah awal dalam memahami karakteristik populasi sebelum dilakukan pengujian lebih lanjut.
Secara umum, estimasi dibedakan menjadi dua jenis utama, yaitu:
- Estimasi titik
- Estimasi interval
3.1.1 Estimasi Titik {sub-header}
Estimasi titik digunakan untuk memberikan satu nilai tertentu yang dianggap sebagai representasi terbaik dari parameter populasi.
Untuk kasus proporsi, estimasi titik dinyatakan sebagai:
\[ \hat{p} = \frac{x}{n} \]
dengan:
- \(\hat{p}\) : estimasi titik
proporsi
- \(x\) : jumlah individu dalam
kategori tertentu
- \(n\) : total jumlah sampel
Estimasi ini bersifat sederhana dan langsung, namun tidak memberikan informasi mengenai tingkat ketidakpastian dari nilai yang dihasilkan.
3.1.2 Estimasi Interval {sub-header}
Berbeda dengan estimasi titik, estimasi interval memberikan rentang nilai yang diyakini mengandung parameter populasi dengan tingkat kepercayaan tertentu.
Bentuk umum interval kepercayaan untuk proporsi adalah:
\[ \hat{p} \pm Z_{\alpha/2} \sqrt{\frac{\hat{p}(1-\hat{p})}{n}} \]
dimana:
- \(Z_{\alpha/2}\) : nilai dari
distribusi normal standar untuk tingkat kepercayaan tertentu
- \(\hat{p}\) : estimasi
proporsi
- \(n\) : ukuran sampel
Pendekatan ini lebih informatif karena tidak hanya memberikan satu nilai, tetapi juga memperlihatkan batas bawah dan batas atas dari kemungkinan nilai parameter populasi.
Kesimpulan:
Estimasi titik memberikan nilai tunggal sebagai dugaan parameter,
sedangkan estimasi interval memberikan rentang nilai yang mencerminkan
ketidakpastian dari estimasi tersebut.
3.2 Uji Hipotesis
Uji hipotesis merupakan metode statistik yang digunakan untuk menentukan apakah terdapat cukup bukti dari data sampel untuk mendukung atau menolak suatu pernyataan mengenai populasi.
Dalam konteks tabel kontingensi dua arah, uji hipotesis digunakan untuk mengevaluasi hubungan antara dua variabel kategorik.
Inferensi pada bagian ini mencakup:
- Uji proporsi
- Uji asosiasi
- Uji independensi
3.2.1 Uji Proporsi
Uji proporsi dua sampel digunakan untuk membandingkan proporsi kejadian antara dua kelompok dalam tabel kontingensi. Analisis ini bertujuan untuk mengetahui apakah terdapat perbedaan yang signifikan antara dua proporsi populasi.
Uji ini banyak digunakan dalam studi kohort dan eksperimen, terutama untuk mengevaluasi perbedaan kejadian antara kelompok perlakuan dan kontrol.
Struktur tabel kontingensi 2 × 2 adalah sebagai berikut:
| Kejadian (+) | Tidak Kejadian (-) | Total | |
|---|---|---|---|
| Grup 1 | \(n_{11}\) | \(n_{12}\) | \(n_{1.}\) |
| Grup 2 | \(n_{21}\) | \(n_{22}\) | \(n_{2.}\) |
| Total | \(n_{.1}\) | \(n_{.2}\) | \(n\) |
Formulasi uji proporsi
Untuk menguji apakah tidak terdapat perbedaan proporsi antara dua kelompok, digunakan uji Z dua proporsi dengan hipotesis:
- \(H_0\): \(p_1 = p_2\)
- \(H_1\): \(p_1 \neq p_2\)
Estimasi proporsi masing-masing kelompok diberikan oleh:
\[ \hat{p}_1 = \frac{n_{11}}{n_{1.}}, \quad \hat{p}_2 = \frac{n_{21}}{n_{2.}} \]
Estimasi proporsi gabungan (pooling proportion):
\[ \hat{p} = \frac{n_{11} + n_{21}}{n_{1.} + n_{2.}} \]
Statistik uji untuk dua proporsi adalah:
\[ Z = \frac{\hat{p}_1 - \hat{p}_2} {\sqrt{\hat{p}(1-\hat{p})\left(\frac{1}{n_{1.}} + \frac{1}{n_{2.}}\right)}} \]
Statistik uji \(Z\) mengikuti distribusi normal baku:
\[ Z \sim N(0,1) \]
Nilai p-value diperoleh berdasarkan distribusi normal tersebut.
Jika \(|Z|\) lebih besar dari nilai kritis pada tingkat signifikansi \(\alpha\) (misalnya 1.96 untuk \(\alpha = 0.05\)), maka hipotesis nol ditolak, yang menunjukkan adanya perbedaan proporsi yang signifikan.
Perhitungan manual langkah demi langkah
Misalkan diberikan data berikut:
| Kejadian (+) | Tidak Kejadian (-) | Total | |
|---|---|---|---|
| Grup 1 | 50 | 30 | 80 |
| Grup 2 | 30 | 50 | 80 |
| Total | 80 | 80 | 160 |
Langkah 1: Hitung proporsi sampel
\[ \hat{p}_1 = \frac{50}{80} = 0.625 \]
\[ \hat{p}_2 = \frac{30}{80} = 0.375 \]
Langkah 2: Hitung proporsi gabungan
\[ \hat{p} = \frac{50 + 30}{80 + 80} = \frac{80}{160} = 0.50 \]
Langkah 3: Hitung statistik uji
\[ Z = \frac{0.625 - 0.375} {\sqrt{0.50(1 - 0.50)\left(\frac{1}{80} + \frac{1}{80}\right)}} \]
\[ Z = \frac{0.25}{\sqrt{0.50 \times 0.50 \times 0.025}} \]
\[ Z = \frac{0.25}{\sqrt{0.00625}} = \frac{0.25}{0.0791} = 3.16 \]
Interpretasi:
Karena \(Z = 3.16 > 1.96\), maka
hipotesis nol ditolak. Artinya terdapat perbedaan proporsi yang
signifikan antara kedua kelompok.
# Pastikan variabel data terdefinisi
set.seed(123)
data <- matrix(c(50, 30, 30, 50), nrow = 2, byrow = TRUE)
dimnames(data) <- list(
"Terpapar" = c("Ya", "Tidak"),
"Kejadian" = c("Ya", "Tidak")
)
# Tampilkan data
data## Kejadian
## Terpapar Ya Tidak
## Ya 50 30
## Tidak 30 50# Uji proporsi dua sampel
prop_test <- prop.test(
x = c(data[1,1], data[2,1]),
n = c(sum(data[1,]), sum(data[2,]))
)
prop_test##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(data[1, 1], data[2, 1]) out of c(sum(data[1, ]), sum(data[2, ]))
## X-squared = 9.025, df = 1, p-value = 0.002663
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.08747151 0.41252849
## sample estimates:
## prop 1 prop 2
## 0.625 0.375Interpretasi hasil: Jika p-value < 0.05, maka terdapat perbedaan proporsi kejadian antara kelompok terpapar dan tidak terpapar.
3.2.2 Uji Asosiasi
Uji asosiasi dalam tabel kontingensi 2 × 2 digunakan untuk mengukur kekuatan hubungan antara dua variabel kategorik. Berbeda dengan uji proporsi yang hanya membandingkan perbedaan, uji asosiasi berfokus pada seberapa kuat keterkaitan antara dua kelompok.
Tiga ukuran utama yang digunakan dalam uji asosiasi adalah:
- Risk Difference (RD) → selisih
risiko absolut
- Relative Risk (RR) →
perbandingan risiko
- Odds Ratio (OR) → perbandingan odds
Hipotesis uji asosiasi
Untuk setiap ukuran asosiasi, hipotesis yang digunakan adalah:
- \(H_0\): tidak terdapat asosiasi
antara dua variabel
- \(H_1\): terdapat asosiasi antara dua variabel
Risk Difference (RD)
Risk Difference mengukur perbedaan absolut probabilitas kejadian antara dua kelompok.
\[ RD = \frac{n_{11}}{n_{1.}} - \frac{n_{21}}{n_{2.}} \]
Standard Error:
\[ SE(RD) = \sqrt{ \frac{\hat{p}_1(1-\hat{p}_1)}{n_{1.}} + \frac{\hat{p}_2(1-\hat{p}_2)}{n_{2.}} } \]
Statistik uji:
\[ Z_{RD} = \frac{RD}{SE(RD)} \]
Relative Risk (RR)
Relative Risk mengukur perbandingan peluang kejadian antara dua kelompok.
\[ RR = \frac{n_{11}/n_{1.}}{n_{21}/n_{2.}} \]
Standard Error untuk log(RR):
\[ SE(\ln RR) = \sqrt{ \frac{1}{n_{11}} - \frac{1}{n_{1.}} + \frac{1}{n_{21}} - \frac{1}{n_{2.}} } \]
Statistik uji:
\[ Z_{RR} = \frac{\ln RR}{SE(\ln RR)} \]
Odds Ratio (OR)
Odds Ratio mengukur perbandingan peluang (odds) kejadian antara dua kelompok.
\[ OR = \frac{n_{11} \times n_{22}}{n_{12} \times n_{21}} \]
Standard Error untuk log(OR):
\[ SE(\ln OR) = \sqrt{ \frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} } \]
Statistik uji:
\[ Z_{OR} = \frac{\ln OR}{SE(\ln OR)} \]
Contoh perhitungan manual
Gunakan data yang sama:
| Kejadian (+) | Tidak Kejadian (-) | Total | |
|---|---|---|---|
| Grup 1 | 50 | 30 | 80 |
| Grup 2 | 30 | 50 | 80 |
| Total | 80 | 80 | 160 |
Misalkan:
\[ \hat{p}_1 = 0.625, \quad \hat{p}_2 = 0.375 \]
Perhitungan Risk Difference:
\[ RD = 0.625 - 0.375 = 0.25 \]
\[ SE(RD) = \sqrt{ \frac{0.625(0.375)}{80} + \frac{0.375(0.625)}{80} } \]
\[ SE(RD) = \sqrt{0.002925 + 0.002925} = \sqrt{0.00585} = 0.0765 \]
\[ Z_{RD} = \frac{0.25}{0.0765} = 3.27 \]
Perhitungan Relative Risk:
\[ RR = \frac{0.625}{0.375} = 1.67 \]
\[ SE(\ln RR) = \sqrt{ \frac{1}{50} - \frac{1}{80} + \frac{1}{30} - \frac{1}{80} } \]
\[ SE(\ln RR) = \sqrt{0.0283} = 0.1683 \]
\[ Z_{RR} = \frac{\ln(1.67)}{0.1683} = \frac{0.51}{0.1683} = 3.03 \]
Perhitungan Odds Ratio:
\[ OR = \frac{50 \times 50}{30 \times 30} = \frac{2500}{900} = 2.78 \]
\[ SE(\ln OR) = \sqrt{ \frac{1}{50} + \frac{1}{30} + \frac{1}{30} + \frac{1}{50} } \]
\[ SE(\ln OR) = \sqrt{0.1066} = 0.3266 \]
\[ Z_{OR} = \frac{\ln(2.78)}{0.3266} = \frac{1.02}{0.3266} = 3.12 \]
Kesimpulan:
- RD menunjukkan perbedaan risiko absolut antar kelompok
- RR menunjukkan perbandingan risiko relatif
- OR menunjukkan perbandingan peluang kejadian
Nilai standard error dan statistik uji Z digunakan untuk menentukan apakah hubungan tersebut signifikan secara statistik.
# Definisi data
n11 <- 50; n12 <- 30; n21 <- 30; n22 <- 50
n1. <- n11 + n12
n2. <- n21 + n22
# Risk Difference
p1 <- n11 / n1.
p2 <- n21 / n2.
rd <- p1 - p2
se_rd <- sqrt((p1*(1-p1)/n1.) + (p2*(1-p2)/n2.))
z_rd <- rd / se_rd
# Relative Risk
rr <- (n11/n1.) / (n21/n2.)
se_ln_rr <- sqrt((1/n11)-(1/n1.)+(1/n21)-(1/n2.))
z_rr <- log(rr) / se_ln_rr
# Odds Ratio
or <- (n11*n22)/(n12*n21)
se_ln_or <- sqrt((1/n11)+(1/n12)+(1/n21)+(1/n22))
z_or <- log(or) / se_ln_or
# Output
list(
RD = rd, SE_RD = se_rd, Z_RD = z_rd,
RR = rr, SE_Ln_RR = se_ln_rr, Z_RR = z_rr,
OR = or, SE_Ln_OR = se_ln_or, Z_OR = z_or
)## $RD
## [1] 0.25
##
## $SE_RD
## [1] 0.07654655
##
## $Z_RD
## [1] 3.265986
##
## $RR
## [1] 1.666667
##
## $SE_Ln_RR
## [1] 0.1683251
##
## $Z_RR
## [1] 3.034756
##
## $OR
## [1] 2.777778
##
## $SE_Ln_OR
## [1] 0.3265986
##
## $Z_OR
## [1] 3.1281553.2.3 Uji Independensi
Uji independensi digunakan untuk menentukan apakah terdapat hubungan statistik antara dua variabel kategorik dalam tabel kontingensi. Jika dua variabel saling independen, maka distribusi salah satu variabel tidak dipengaruhi oleh variabel lainnya.
Uji Chi-Square
Uji Chi-Square digunakan untuk menguji apakah terdapat hubungan antara dua variabel kategorik.
Rumus statistik uji:
\[ \chi^2 = \sum \frac{(O - E)^2}{E} \]
dengan: - \(O\) = frekuensi
observasi
- \(E\) = frekuensi harapan
Frekuensi harapan dihitung sebagai:
\[ E_{ij} = \frac{R_i \times C_j}{N} \]
dimana: - \(R_i\) = total baris
ke-i
- \(C_j\) = total kolom ke-j
- \(N\) = total sampel
Contoh perhitungan manual Chi-Square
| Ya | Tidak | Total | |
|---|---|---|---|
| Terpapar | 30 | 10 | 40 |
| Tidak Terpapar | 15 | 45 | 60 |
| Total | 45 | 55 | 100 |
Hitung nilai ekspektasi:
\[ E_{11} = \frac{40 \times 45}{100} = 18 \]
\[ E_{12} = \frac{40 \times 55}{100} = 22 \]
\[ E_{21} = \frac{60 \times 45}{100} = 27 \]
\[ E_{22} = \frac{60 \times 55}{100} = 33 \]
Hitung statistik uji:
\[ \chi^2 = \frac{(30-18)^2}{18} + \frac{(10-22)^2}{22} + \frac{(15-27)^2}{27} + \frac{(45-33)^2}{33} \]
\[ = 8 + 6.55 + 5.33 + 4.36 = 24.24 \]
Derajat bebas:
\[ df = (2-1)(2-1) = 1 \]
Interpretasi:
Karena \(\chi^2 = 24.24 > 3.841\),
maka hipotesis nol ditolak → terdapat hubungan antara variabel.
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: data
## X-squared = 22.264, df = 1, p-value = 2.376e-06Jika p-value < 0.05, maka terdapat hubungan signifikan antara variabel.
Partisi Chi-Square
Partisi Chi-Square digunakan untuk mengidentifikasi bagian mana dari tabel yang memberikan kontribusi terhadap hubungan signifikan.
Konsep ini juga berkaitan dengan Simpson’s Paradox, yaitu kondisi dimana pola hubungan dapat berubah ketika data digabungkan.
Langkah-langkah yang dilakukan dalam partisi Chi-Square adalah:
- Hitung nilai Chi-Square keseluruhan
- Pecah tabel kontingensi menjadi beberapa tabel 2 × 2
- Hitung nilai Chi-Square untuk masing-masing bagian
- Interpretasikan kontribusi masing-masing bagian terhadap hubungan keseluruhan
Contoh (Agresti)
| Gender | Democrat | Republican | Independent | Total |
|---|---|---|---|---|
| Female | 495 | 272 | 590 | 1357 |
| Male | 330 | 265 | 498 | 1093 |
| Total | 825 | 537 | 1088 | 2450 |
Hasil uji Chi-Square
keseluruhan:
Nilai \(\chi^2 = 12.57\) dengan nilai
kritis \(5.99\), sehingga hasilnya signifikan.
Partisi yang dilakukan:
- Democrat vs Republican → signifikan
- (Democrat + Republican) vs Independent → tidak signifikan
# Partisi 1: Democrat vs Republican
data1 <- matrix(c(495,272,330,265), nrow=2, byrow=TRUE)
chisq.test(data1)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: data1
## X-squared = 11.178, df = 1, p-value = 0.0008279# Partisi 2: (Democrat + Republican) vs Independent
data2 <- matrix(c(767,590,595,498), nrow=2, byrow=TRUE)
chisq.test(data2)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: data2
## X-squared = 0.98267, df = 1, p-value = 0.3215Uji Likelihood Ratio (G²)
Uji Likelihood Ratio (G²) merupakan alternatif dari uji Chi-Square yang digunakan untuk menguji hubungan antara dua variabel kategorik dalam tabel kontingensi.
Statistik uji yang digunakan adalah:
\[ G^2 = 2 \sum n_{ij} \ln \left(\frac{n_{ij}}{\hat{\mu}_{ij}}\right) \]
dengan: - \(n_{ij}\) = frekuensi
observasi
- \(\hat{\mu}_{ij}\) = frekuensi
harapan
Interpretasi:
Jika nilai \(G^2\) lebih besar dari
nilai kritis Chi-Square, maka hipotesis nol
ditolak, yang berarti terdapat hubungan antara variabel.
Contoh perhitungan menggunakan R
# Membuat data
data <- matrix(c(688,650,21,59), nrow=2, byrow=TRUE)
# Menghitung nilai ekspektasi
expected <- chisq.test(data)$expected
# Menghitung G^2
G2 <- 2 * sum(data * log(data / expected))
G2## [1] 19.87802Jika \(G^2 > \chi^2_{kritik}\) maka tolak \(H_0\).
Uji Fisher Exact
Uji Fisher Exact digunakan ketika ukuran sampel relatif kecil sehingga asumsi uji Chi-Square tidak terpenuhi. Uji ini memberikan hasil yang lebih akurat karena tidak bergantung pada pendekatan distribusi normal.
Keunggulan uji Fisher Exact:
- Tidak memerlukan asumsi distribusi
normal
- Lebih akurat untuk frekuensi kecil
Distribusi Hipergeometrik
Dasar dari uji Fisher Exact adalah distribusi hipergeometrik, yang digunakan untuk menghitung probabilitas dari suatu konfigurasi tabel kontingensi.
Rumus distribusi hipergeometrik:
\[ P(X = x) = \frac{\binom{K}{x} \binom{N-K}{n-x}}{\binom{N}{n}} \]
dimana:
- \(N\) = total populasi
- \(K\) = jumlah elemen dalam
kategori tertentu
- \(n\) = ukuran sampel
- \(x\) = jumlah kejadian dalam sampel
Interpretasi:
Probabilitas dihitung untuk setiap kemungkinan tabel yang memiliki
margin yang sama, kemudian dijumlahkan untuk memperoleh p-value pada uji Fisher Exact.
Contoh perhitungan
## [1] 0.01380413Perhitungan manual probabilitas tabel
## [1] 0.01380413p-value dihitung sebagai jumlah probabilitas dari tabel yang sama atau lebih ekstrem.
Uji Exact Fisher
data <- matrix(c(18,2,11,9), nrow=2, byrow=TRUE) fisher.test(data)
Interpretasi:
- p-value < 0.05 → terdapat hubungan
- p-value ≥ 0.05 → tidak terdapat bukti hubungan
Kesimpulan umum
- Uji Chi-Square → digunakan untuk melihat hubungan secara umum
- Partisi Chi-Square → mengidentifikasi sumber hubungan
- G² → alternatif dari Chi-Square
- Fisher Exact → digunakan untuk sampel berukuran kecil
3.2.4 Tugas Mahasiswa (Latihan)
Tugas:
Buatlah fungsi untuk menghitung dan melakukan pengujian hipotesis untuk
Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR).
Gunakan data berikut (Agresti, 2019):
| Smoker | Lung Cancer (Cases) | Control |
|---|---|---|
| Yes | 688 | 650 |
| No | 21 | 59 |
Struktur tabel
| Exposure | Cases | Control | Total |
|---|---|---|---|
| Yes | a | c | a+c |
| No | b | d | b+d |
| Total | a+b | c+d | a+b+c+d |
Fungsi Risk Difference (RD)
prop_diff <- function(a, b, c, d, alpha = 0.05) {
ph <- a / (a + c)
pi <- b / (b + d)
nh <- a + c
ni <- b + d
se_bp <- sqrt((ph * (1 - ph) / nh) + (pi * (1 - pi) / ni))
z_alpha <- qnorm(1 - alpha / 2)
ci_lower <- (ph - pi) - z_alpha * se_bp
ci_upper <- (ph - pi) + z_alpha * se_bp
list(estimate = ph - pi, ci = c(ci_lower, ci_upper))
}
hasil <- prop_diff(a = 688, b = 21, c = 650, d = 59)
print(hasil)## $estimate
## [1] 0.2517003
##
## $ci
## [1] 0.1516343 0.3517663Fungsi Relative Risk (RR)
relative_risk <- function(a, b, c, d, alpha = 0.05) {
ph <- a / (a + c)
pi <- b / (b + d)
nh <- a + c
ni <- b + d
ln_rr <- log(ph / pi)
se_ln_rr <- sqrt(((1 - ph) / (ph * nh)) + ((1 - pi) / (pi * ni)))
z_alpha <- qnorm(1 - alpha / 2)
ci_lower <- exp(ln_rr - z_alpha * se_ln_rr)
ci_upper <- exp(ln_rr + z_alpha * se_ln_rr)
list(estimate = exp(ln_rr), ci = c(ci_lower, ci_upper))
}
hasil <- relative_risk(a = 688, b = 21, c = 650, d = 59)
print(hasil)## $estimate
## [1] 1.958858
##
## $ci
## [1] 1.351735 2.838667Fungsi Odds Ratio (OR)
odds_ratio <- function(a, b, c, d, alpha = 0.05) {
ln_or <- log((a * d) / (b * c))
se_ln_or <- sqrt(1/a + 1/b + 1/c + 1/d)
z_alpha <- qnorm(1 - alpha / 2)
ci_lower <- exp(ln_or - z_alpha * se_ln_or)
ci_upper <- exp(ln_or + z_alpha * se_ln_or)
list(estimate = exp(ln_or), ci = c(ci_lower, ci_upper))
}
hasil <- odds_ratio(a = 688, b = 21, c = 650, d = 59)
print(hasil)## $estimate
## [1] 2.973773
##
## $ci
## [1] 1.786737 4.949427Interpretasi hasil
Risk Difference (RD) = 0.2517 → Risiko kanker paru pada perokok lebih tinggi sekitar 25.17% secara absolut
Relative Risk (RR) = 1.96 → Perokok memiliki risiko sekitar 1.96 kali lebih tinggi
Odds Ratio (OR) = 2.97 → Odds kejadian pada perokok sekitar 2.97 kali lebih besar
Perhitungan manual
a <- 688
b <- 21
c <- 650
d <- 59
# Risk Difference
RD_manual <- (a / (a + c)) - (b / (b + d))
SE_RD <- sqrt((a/(a+c)*(1 - a/(a+c)))/(a+c) + (b/(b+d)*(1 - b/(b+d)))/(b+d))
CI_RD <- c(RD_manual - 1.96 * SE_RD, RD_manual + 1.96 * SE_RD)
# Relative Risk
RR_manual <- (a / (a + c)) / (b / (b + d))
SE_RR <- sqrt(1/a - 1/(a+c) + 1/b - 1/(b+d))
CI_RR <- exp(log(RR_manual) + c(-1.96, 1.96) * SE_RR)
# Odds Ratio
OR_manual <- (a * d) / (b * c)
SE_OR <- sqrt(1/a + 1/b + 1/c + 1/d)
CI_OR <- exp(log(OR_manual) + c(-1.96, 1.96) * SE_OR)
list(RD = RD_manual, CI_RD = CI_RD,
RR = RR_manual, CI_RR = CI_RR,
OR = OR_manual, CI_OR = CI_OR)## $RD
## [1] 0.2517003
##
## $CI_RD
## [1] 0.1516324 0.3517682
##
## $RR
## [1] 1.958858
##
## $CI_RR
## [1] 1.351726 2.838687
##
## $OR
## [1] 2.973773
##
## $CI_OR
## [1] 1.786720 4.949474Perbandingan dengan output R
# Data
a <- 688
b <- 21
c <- 650
d <- 59
# Risk Difference
RD <- (a / (a + c)) - (b / (b + d))
# Relative Risk
RR <- (a / (a + c)) / (b / (b + d))
# Odds Ratio
OR <- (a * d) / (b * c)
# Output
list(
Risk_Difference = RD,
Relative_Risk = RR,
Odds_Ratio = OR
)## $Risk_Difference
## [1] 0.2517003
##
## $Relative_Risk
## [1] 1.958858
##
## $Odds_Ratio
## [1] 2.973773Kesimpulan: Risk Difference, Relative Risk, dan Odds Ratio memberikan gambaran mengenai kekuatan hubungan antara paparan dan kejadian. Hasil perhitungan manual dan menggunakan R yang konsisten menunjukkan bahwa metode yang digunakan telah valid.
4 Tugas “Inferensi Tabel Kontingensi Dua Arah”
4.1 Kasus 1: Tabel Kontingensi 2x2
4.1.1 Nomor 1: Penyusunan Tabel Kontingensi 2x2
Penyusunan Secara Manual
Tabel kontingensi 2×2 digunakan untuk menyajikan hubungan antara dua variabel kategorik, yaitu:
- Status merokok (Smoker vs Non-Smoker)
- Status kanker paru (Cancer (+) vs Control (-))
Berdasarkan data yang diberikan, diperoleh tabel sebagai berikut:
| Status Merokok | Cancer (+) | Control (-) | Total |
|---|---|---|---|
| Smoker | 688 | 650 | 1338 |
| Non-Smoker | 21 | 59 | 80 |
| Total | 709 | 709 | 1418 |
Penyusunan Menggunakan R
# Membuat tabel kontingensi
tabel <- matrix(c(688,650,21,59),
nrow = 2,
byrow = TRUE)
colnames(tabel) <- c("Cancer(+)", "Control(-)")
rownames(tabel) <- c("Smoker", "Non-Smoker")
tabel## Cancer(+) Control(-)
## Smoker 688 650
## Non-Smoker 21 59## Smoker Non-Smoker
## 1338 80## Cancer(+) Control(-)
## 709 709## [1] 14184.1.2 Nomor 2: Estimasi Titik Proporsi Kejadian Kanker Paru
Perhitungan Secara Manual
Proporsi kejadian kanker paru pada masing-masing kelompok dihitung sebagai:
Proporsi pada kelompok smoker:
\[ p_1 = \frac{\text{jumlah smoker dengan kanker}}{\text{total smoker}} = \frac{688}{1338} \]
Proporsi pada kelompok non-smoker:
\[ p_2 = \frac{\text{jumlah non-smoker dengan kanker}}{\text{total non-smoker}} = \frac{21}{80} \]
Sehingga diperoleh:
- \(p_1 = \frac{688}{1338} \approx 0.5142\)
- \(p_2 = \frac{21}{80} = 0.2625\)
Perhitungan Menggunakan R
# Ambil nilai dari tabel
n_smoker <- sum(tabel["Smoker", ])
n_nonsmoker <- sum(tabel["Non-Smoker", ])
x_smoker <- tabel["Smoker", "Cancer(+)"]
x_nonsmoker <- tabel["Non-Smoker", "Cancer(+)"]
# Hitung proporsi
p_smoker <- x_smoker / n_smoker
p_nonsmoker <- x_nonsmoker / n_nonsmoker
# Tampilkan hasil
p_smoker## [1] 0.5142003## [1] 0.2625Interpretasi
Proporsi kejadian kanker paru pada kelompok smoker adalah sebesar 0.5142 (≈ 51.42%), sedangkan pada kelompok non-smoker sebesar 0.2625 (≈ 26.25%).
Secara deskriptif, terlihat bahwa proporsi kejadian kanker paru pada kelompok smoker lebih tinggi dibandingkan non-smoker, dengan selisih sekitar 0.2517 (≈ 25.17%). Perbedaan ini memberikan indikasi awal adanya hubungan antara kebiasaan merokok dan kejadian kanker paru, yang selanjutnya akan diuji secara inferensial pada tahap berikutnya.
4.1.3 Nomor 3: Interval Kepercayaan 95% untuk Proporsi, RD, RR, dan OR
Perhitungan Secara Konseptual
Interval kepercayaan (confidence interval) digunakan untuk mengestimasi rentang nilai parameter populasi berdasarkan sampel.
Pada kasus ini dihitung:
- Interval kepercayaan proporsi masing-masing kelompok
- Risk Difference (RD): \(p_1 - p_2\)
- Relative Risk (RR): \(\frac{p_1}{p_2}\)
- Odds Ratio (OR): \(\frac{ad}{bc}\)
dengan:
- \(p_1\): proporsi kanker pada
smoker
- \(p_2\): proporsi kanker pada non-smoker
Perhitungan Menggunakan R
# Ambil nilai dari tabel
n_smoker <- sum(tabel["Smoker", ])
n_nonsmoker <- sum(tabel["Non-Smoker", ])
x_smoker <- tabel["Smoker", "Cancer(+)"]
x_nonsmoker <- tabel["Non-Smoker", "Cancer(+)"]
# Proporsi
p_smoker <- x_smoker / n_smoker
p_nonsmoker <- x_nonsmoker / n_nonsmoker
# CI proporsi
ci_smoker <- prop.test(x_smoker, n_smoker)
ci_nonsmoker <- prop.test(x_nonsmoker, n_nonsmoker)
ci_smoker##
## 1-sample proportions test with continuity correction
##
## data: x_smoker out of n_smoker, null probability 0.5
## X-squared = 1.0232, df = 1, p-value = 0.3118
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.4870445 0.5412736
## sample estimates:
## p
## 0.5142003##
## 1-sample proportions test with continuity correction
##
## data: x_nonsmoker out of n_nonsmoker, null probability 0.5
## X-squared = 17.113, df = 1, p-value = 3.523e-05
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.1733064 0.3748263
## sample estimates:
## p
## 0.2625# Risk Difference
RD <- p_smoker - p_nonsmoker
# Relative Risk
RR <- p_smoker / p_nonsmoker
# Odds Ratio
a <- tabel[1,1]; b <- tabel[1,2]
c <- tabel[2,1]; d <- tabel[2,2]
OR <- (a*d)/(b*c)
RD## [1] 0.2517003## [1] 1.958858## [1] 2.973773Interpretasi
Interval kepercayaan 95% untuk proporsi pada kelompok smoker berada pada rentang 0.487 hingga 0.5413, sedangkan pada kelompok non-smoker berada pada rentang 0.1733 hingga 0.3748.
Nilai Risk Difference (RD) sebesar 0.2517 menunjukkan bahwa terdapat peningkatan risiko absolut kejadian kanker paru sebesar sekitar 25.17%.
Nilai Relative Risk (RR) sebesar 1.9589 menunjukkan bahwa individu yang merokok memiliki risiko sekitar 1.96 kali lebih besar dibandingkan non-smoker.
Nilai Odds Ratio (OR) sebesar 2.9738 menunjukkan bahwa odds kejadian kanker paru pada kelompok smoker sekitar 2.97 kali dibandingkan non-smoker.
Secara keseluruhan, karena nilai RR dan OR lebih besar dari 1, hal ini mengindikasikan adanya hubungan positif antara kebiasaan merokok dan kejadian kanker paru. Jika interval kepercayaan untuk RR dan OR tidak mencakup nilai 1, maka hubungan tersebut bersifat signifikan secara statistik.
4.1.4 Nomor 4: Uji Dua Proporsi
Konsep dan Hipotesis
Uji dua proporsi digunakan untuk menguji apakah terdapat perbedaan proporsi kejadian kanker paru antara kelompok smoker dan non-smoker.
Hipotesis yang digunakan:
- \(H_0\): \(p_1 = p_2\) (tidak ada perbedaan proporsi)
- \(H_1\): \(p_1 \neq p_2\) (terdapat perbedaan proporsi)
Pengujian Menggunakan R
# jumlah kasus
x <- c(
tabel["Smoker", "Cancer(+)"],
tabel["Non-Smoker", "Cancer(+)"]
)
# total masing-masing kelompok
n <- c(
sum(tabel["Smoker", ]),
sum(tabel["Non-Smoker", ])
)
# uji dua proporsi
uji_prop <- prop.test(x, n, correct = FALSE)
uji_prop##
## 2-sample test for equality of proportions without continuity correction
##
## data: x out of n
## X-squared = 19.129, df = 1, p-value = 1.222e-05
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.1516343 0.3517663
## sample estimates:
## prop 1 prop 2
## 0.5142003 0.2625000Interpretasi
Nilai p-value yang diperoleh adalah sebesar 1.2^{-5}.
Karena p-value < 0.05, maka \(H_0\) ditolak.
Hal ini menunjukkan bahwa terdapat perbedaan proporsi kejadian kanker paru yang signifikan antara kelompok smoker dan non-smoker.
Secara substantif, hasil ini mengindikasikan bahwa kebiasaan merokok berhubungan dengan peningkatan kejadian kanker paru, sehingga proporsi kejadian pada kelompok smoker secara signifikan lebih tinggi dibandingkan non-smoker.
4.1.5 Nomor 5: Uji Chi-Square Independensi
Konsep dan Hipotesis
Uji chi-square independensi digunakan untuk menguji apakah terdapat hubungan antara dua variabel kategorik, yaitu status merokok dan kejadian kanker paru.
Hipotesis yang digunakan:
- H0: Status merokok dan kejadian kanker paru saling
independen
- H1: Status merokok dan kejadian kanker paru tidak independen
Pengujian Menggunakan R
##
## Pearson's Chi-squared test
##
## data: tabel
## X-squared = 19.129, df = 1, p-value = 1.222e-05Interpretasi
Nilai p-value yang diperoleh adalah sebesar 1.2^{-5}.
Karena p-value < 0.05, maka H0 ditolak.
Hal ini menunjukkan bahwa terdapat hubungan yang signifikan antara status merokok dan kejadian kanker paru.
Secara substantif, hasil ini mengindikasikan bahwa kejadian kanker paru tidak terjadi secara independen terhadap status merokok, di mana individu pada kelompok smoker memiliki kecenderungan lebih tinggi untuk mengalami kanker paru dibandingkan non-smoker.
4.1.6 Nomor 6: Uji Likelihood Ratio (G²)
Konsep dan Hipotesis
Uji likelihood ratio (G²) digunakan untuk menguji apakah terdapat hubungan antara dua variabel kategorik berdasarkan pendekatan likelihood.
Hipotesis yang digunakan:
- H0: Status merokok dan kejadian kanker paru saling
independen
- H1: Status merokok dan kejadian kanker paru tidak independen
Pengujian Menggunakan R
## Warning: package 'DescTools' was built under R version 4.4.3##
## Log likelihood ratio (G-test) test of independence without correction
##
## data: tabel
## G = 19.878, X-squared df = 1, p-value = 8.254e-06Interpretasi
Nilai p-value yang diperoleh adalah sebesar 8^{-6}.
Karena p-value < 0.05, maka H0 ditolak.
Hal ini menunjukkan bahwa model yang mengasumsikan tidak adanya hubungan (independensi) antara status merokok dan kejadian kanker paru tidak sesuai dengan data.
Dengan demikian, terdapat hubungan yang signifikan antara status merokok dan kejadian kanker paru, di mana kelompok smoker menunjukkan kecenderungan lebih tinggi mengalami kanker paru dibandingkan non-smoker.
4.1.7 Nomor 7: Fisher Exact Test
Konsep dan Hipotesis
Fisher exact test digunakan untuk menguji hubungan antara dua variabel kategorik pada tabel kontingensi, khususnya ketika ukuran sampel kecil atau terdapat frekuensi harapan yang rendah.
Hipotesis yang digunakan:
- H0: Status merokok dan kejadian kanker paru saling
independen
- H1: Status merokok dan kejadian kanker paru tidak independen
Pengujian Menggunakan R
##
## Fisher's Exact Test for Count Data
##
## data: tabel
## p-value = 1.476e-05
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.755611 5.210711
## sample estimates:
## odds ratio
## 2.971634Interpretasi
Nilai p-value yang diperoleh adalah sebesar
rround(uji_fisher$p.value,6)`.
Karena p-value < 0.05, maka H0 ditolak.
Hal ini menunjukkan bahwa terdapat hubungan yang signifikan antara status merokok dan kejadian kanker paru.
Fisher exact test memberikan hasil yang konsisten dengan uji sebelumnya, sehingga memperkuat kesimpulan bahwa kejadian kanker paru tidak independen terhadap status merokok, di mana kelompok smoker memiliki kecenderungan lebih tinggi mengalami kanker paru dibandingkan non-smoker.
4.1.8 Nomor 8: Perbandingan Hasil Uji
Perbandingan Metode Uji
Berikut adalah perbandingan hasil uji dua proporsi, chi-square, likelihood ratio (G²), dan Fisher exact test:
| Metode | Hipotesis | Statistik Uji | p-value | Keputusan |
|---|---|---|---|---|
| Uji Dua Proporsi | \(p_1 = p_2\) | Z / Chi-square | 1.22^{-5} | Tolak H0 |
| Chi-Square | Independen | \(\chi^2\) = 19.1292 | 1.22^{-5} | Tolak H0 |
| Likelihood Ratio | Independen | G² = 19.878 | 8.25^{-6} | Tolak H0 |
| Fisher Exact | Independen | Exact Test | 1.48^{-5} | Tolak H0 |
Interpretasi Substantif
Keempat metode pengujian menghasilkan keputusan yang konsisten, yaitu menolak H0, sehingga menunjukkan adanya hubungan atau perbedaan yang signifikan antara variabel yang diteliti.
Uji dua proporsi secara khusus menunjukkan adanya perbedaan proporsi kejadian kanker paru antara kelompok smoker dan non-smoker, sedangkan uji chi-square, likelihood ratio, dan Fisher exact test menunjukkan bahwa kedua variabel tidak bersifat independen.
Meskipun pendekatan yang digunakan berbeda, yaitu pendekatan proporsi, aproksimasi chi-square, likelihood, dan exact test, seluruh metode memberikan hasil yang sejalan.
Secara substantif, hal ini memperkuat kesimpulan bahwa kebiasaan merokok memiliki hubungan yang signifikan dengan peningkatan kejadian kanker paru, di mana kelompok smoker secara konsisten menunjukkan risiko yang lebih tinggi dibandingkan non-smoker.
Konsistensi hasil dari berbagai metode ini meningkatkan kepercayaan terhadap validitas kesimpulan yang diperoleh.
4.1.9 Nomor 9: Kesimpulan
Berdasarkan hasil analisis yang telah dilakukan, baik melalui estimasi proporsi, ukuran asosiasi (RD, RR, dan OR), maupun berbagai uji hipotesis (uji dua proporsi, chi-square, likelihood ratio, dan Fisher exact test), diperoleh hasil yang konsisten.
Seluruh metode menunjukkan bahwa terdapat hubungan yang signifikan antara status merokok dan kejadian kanker paru.
Secara substantif, individu pada kelompok smoker memiliki risiko yang lebih tinggi untuk mengalami kanker paru dibandingkan dengan kelompok non-smoker.
Dengan demikian, dapat disimpulkan bahwa kebiasaan merokok berhubungan dengan peningkatan kejadian kanker paru.
4.2 Kasus 2: Tabel Kontingensi 2x3
4.2.1 Nomor 1: Penyusunan Tabel Kontingensi 2x3
Penyusunan Secara Manual
Tabel kontingensi 2×3 digunakan untuk menganalisis hubungan antara dua variabel kategorik, yaitu:
- Gender (female dan male)
- Identifikasi partai politik (Democrat, Republican, dan Independent)
Berdasarkan data yang diberikan, diperoleh tabel sebagai berikut:
| Gender | Democrat | Republican | Independent | Total |
|---|---|---|---|---|
| Female | 495 | 272 | 590 | 1357 |
| Male | 330 | 265 | 498 | 1093 |
| Total | 825 | 537 | 1088 | 2450 |
Tabel ini merupakan tabel kontingensi 2×3 karena terdiri dari dua kategori pada variabel gender dan tiga kategori pada variabel identifikasi partai politik.
Penyusunan Menggunakan R
# Membuat tabel kontingensi
tabel2 <- matrix(c(495,272,590,
330,265,498),
nrow = 2,
byrow = TRUE)
colnames(tabel2) <- c("Democrat", "Republican", "Independent")
rownames(tabel2) <- c("Female", "Male")
tabel2## Democrat Republican Independent
## Female 495 272 590
## Male 330 265 498## Female Male
## 1357 1093## Democrat Republican Independent
## 825 537 1088## [1] 2450Interpretasi
Tabel kontingensi menunjukkan distribusi identifikasi partai politik berdasarkan gender.
Secara deskriptif, baik kelompok female maupun male memiliki jumlah responden terbesar pada kategori Independent dibandingkan kategori lainnya.
Namun demikian, untuk menentukan apakah perbedaan distribusi tersebut bersifat signifikan secara statistik, diperlukan analisis inferensial lebih lanjut.
4.2.2 Nomor 2: Frekuensi Harapan
Perhitungan Secara Konseptual
Frekuensi harapan (expected frequency) pada setiap sel tabel kontingensi dihitung dengan rumus:
\[ E_{ij} = \frac{(\text{Total baris ke-i}) \times (\text{Total kolom ke-j})}{\text{Total keseluruhan}} \]
di mana: \(E_{ij}\) adalah frekuensi harapan pada baris ke-i dan kolom ke-j
Sebagai contoh, frekuensi harapan untuk kategori female dan Democrat adalah:
\[ E = \frac{(1357 \times 825)}{2450} \]
Perhitungan Menggunakan R
## Democrat Republican Independent
## Female 456.949 297.4322 602.6188
## Male 368.051 239.5678 485.3812Interpretasi
Frekuensi harapan menunjukkan jumlah kasus yang diharapkan pada setiap sel jika tidak terdapat hubungan antara variabel gender dan identifikasi partai politik.
Berdasarkan hasil perhitungan, seluruh nilai frekuensi harapan berada di atas 5, sehingga memenuhi asumsi untuk penggunaan uji chi-square.
Hal ini menunjukkan bahwa data yang digunakan cukup memadai untuk dilakukan analisis chi-square independensi pada tahap selanjutnya.
4.2.3 Nomor 3: Uji Chi-Square Independensi
Konsep dan Hipotesis
Uji chi-square independensi digunakan untuk menguji apakah terdapat hubungan antara dua variabel kategorik, yaitu gender dan identifikasi partai politik.
Hipotesis yang digunakan:
- H0: Gender dan identifikasi partai politik saling
independen
- H1: Gender dan identifikasi partai politik tidak independen
Pengujian Menggunakan R
##
## Pearson's Chi-squared test
##
## data: tabel2
## X-squared = 12.569, df = 2, p-value = 0.001865Interpretasi
Nilai statistik chi-square yang diperoleh adalah sebesar 12.5693 dengan p-value sebesar 0.00186.
Karena p-value < 0.05, maka H0 ditolak.
Hal ini menunjukkan bahwa terdapat hubungan yang signifikan antara gender dan identifikasi partai politik.
Secara substantif, distribusi preferensi terhadap Democrat, Republican, dan Independent berbeda antara kelompok female dan male, sehingga kedua variabel tidak bersifat independen.
Perbedaan ini mencerminkan adanya variasi preferensi politik berdasarkan gender dalam populasi yang diamati.
4.2.4 Nomor 4: Residual Pearson (Standardized Residual)
Konsep
Residual Pearson digunakan untuk mengukur selisih antara frekuensi observasi dan frekuensi harapan pada setiap sel dalam tabel kontingensi.
Residual dihitung sebagai:
\[ r_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}} \]
Nilai residual yang besar (dalam nilai absolut), umumnya lebih dari 2, menunjukkan bahwa sel tersebut berkontribusi besar terhadap hasil uji chi-square.
Perhitungan Menggunakan R
## Democrat Republican Independent
## Female 1.780051 -1.474656 -0.5140388
## Male -1.983409 1.643125 0.5727640Interpretasi
Nilai residual menunjukkan kontribusi masing-masing sel terhadap hubungan antara gender dan identifikasi partai politik.
Sel dengan nilai residual terbesar (dalam nilai absolut) merupakan sel yang paling berkontribusi terhadap hasil uji chi-square.
Berdasarkan hasil perhitungan, kategori Independent menunjukkan perbedaan yang cukup mencolok antara kelompok female dan male.
Nilai residual yang relatif besar menunjukkan bahwa jumlah observasi pada kategori tersebut berbeda cukup jauh dari yang diharapkan jika kedua variabel bersifat independen.
Dengan demikian, kategori Independent menjadi salah satu kontributor utama dalam hubungan antara gender dan preferensi partai politik.
4.2.5 Nomor 5: Partisi Chi-Square
Konsep
Partisi chi-square dilakukan untuk mengidentifikasi bagian mana dari tabel kontingensi yang paling berkontribusi terhadap hubungan antara variabel.
Pada kasus ini, partisi dilakukan menjadi dua bagian:
- Perbandingan antara Democrat dan
Republican
- Perbandingan antara (Democrat + Republican) dan Independent
Partisi 1: Democrat vs Republican
# Subtabel Democrat vs Republican
tabel_DR <- tabel2[, c("Democrat", "Republican")]
uji_DR <- chisq.test(tabel_DR, correct = FALSE)
uji_DR##
## Pearson's Chi-squared test
##
## data: tabel_DR
## X-squared = 11.555, df = 1, p-value = 0.0006758Partisi 2: (Democrat + Republican) vs Independent
# Gabungkan Democrat + Republican
gabungan <- tabel2[, "Democrat"] + tabel2[, "Republican"]
tabel_DI <- cbind(gabungan, tabel2[, "Independent"])
colnames(tabel_DI) <- c("Democrat+Republican", "Independent")
uji_DI <- chisq.test(tabel_DI, correct = FALSE)
uji_DI##
## Pearson's Chi-squared test
##
## data: tabel_DI
## X-squared = 1.0654, df = 1, p-value = 0.302Interpretasi
Pada partisi pertama, diperoleh p-value sebesar 6.76^{-4}.
Pada partisi kedua, diperoleh p-value sebesar 0.302.
Jika dibandingkan, partisi dengan p-value yang lebih kecil menunjukkan kontribusi yang lebih besar terhadap hubungan keseluruhan.
Hasil menunjukkan bahwa perbedaan antara kategori Democrat dan Republican relatif lebih kecil dibandingkan dengan perbedaan antara kelompok (Democrat + Republican) dan Independent.
Dengan demikian, kategori Independent memberikan kontribusi yang lebih besar terhadap hubungan antara gender dan preferensi partai politik.
4.2.6 Nomor 6: Perbandingan Hasil Partisi dengan Uji Chi-Square Keseluruhan
Perbandingan Hasil
Uji chi-square keseluruhan menghasilkan nilai statistik sebesar 12.5693 dengan p-value sebesar 0.00186, yang menunjukkan adanya hubungan yang signifikan antara gender dan identifikasi partai politik.
Sementara itu, hasil partisi chi-square menunjukkan:
- Partisi Democrat vs Republican memiliki p-value sebesar 6.76^{-4}
- Partisi (Democrat + Republican) vs Independent memiliki p-value sebesar 0.302
Interpretasi Substantif
Hasil uji chi-square keseluruhan menunjukkan bahwa terdapat hubungan antara gender dan preferensi partai politik, namun tidak menunjukkan bagian mana yang paling berkontribusi.
Melalui partisi chi-square, dapat diketahui bahwa kontribusi terbesar terhadap hubungan tersebut berasal dari perbedaan antara kategori Independent dengan gabungan kategori Democrat dan Republican.
Sebaliknya, perbedaan antara kategori Democrat dan Republican relatif lebih kecil.
Dengan demikian, hasil partisi memberikan informasi yang lebih rinci dibandingkan uji chi-square keseluruhan, yaitu bahwa variasi preferensi terhadap kategori Independent menjadi faktor utama yang membedakan distribusi antara kelompok female dan male.
Hal ini menunjukkan bahwa analisis partisi chi-square memberikan insight yang lebih mendalam terhadap struktur hubungan dalam data.
4.2.7 Nomor 7: Kategori yang Paling Berkontribusi
Dasar Analisis
Untuk menentukan kategori yang paling berkontribusi terhadap hubungan antara gender dan identifikasi partai politik, digunakan nilai residual Pearson dari setiap sel.
Sel dengan nilai residual terbesar (dalam nilai absolut) menunjukkan kontribusi paling besar terhadap hasil uji chi-square.
Identifikasi Menggunakan R
# Ambil residual
residual <- chisq.test(tabel2)$residuals
# Cari posisi residual terbesar (nilai absolut)
which(abs(residual) == max(abs(residual)), arr.ind = TRUE)## row col
## Male 2 1Interpretasi
Berdasarkan hasil residual, sel dengan kontribusi terbesar berasal dari kategori Independent, khususnya pada perbedaan antara kelompok female dan male.
Hal ini menunjukkan bahwa jumlah individu pada kategori Independent menyimpang cukup besar dari nilai yang diharapkan jika tidak terdapat hubungan antara variabel.
Temuan ini konsisten dengan hasil partisi chi-square, yang menunjukkan bahwa perbedaan antara kategori (Democrat + Republican) dan Independent memberikan kontribusi terbesar terhadap hubungan keseluruhan.
Dengan demikian, dapat disimpulkan bahwa kategori Independent merupakan faktor utama yang berkontribusi terhadap hubungan antara gender dan identifikasi partai politik.
5 Regresi Logistik Biner: Bank Marketing
5.1 Ringkasan Masalah
Dataset Bank Marketing digunakan untuk memprediksi apakah seorang nasabah akan berlangganan deposito berjangka setelah mendapatkan kampanye pemasaran melalui telepon.
Variabel respon yang digunakan adalah y, dengan dua kategori:
- yes: nasabah berlangganan deposito berjangka
- no: nasabah tidak berlangganan deposito berjangka
Karena variabel respon hanya memiliki dua kategori, metode yang sesuai adalah regresi logistik biner.
Model umum yang digunakan adalah:
\[ \log\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \cdots + \beta_kX_{ki}, \]
dengan \(p_i = P(Y_i = 1 \mid X_i)\), yaitu peluang nasabah ke-\(i\) berlangganan deposito berjangka.
Dalam analisis ini, \(Y_i = 1\) berarti nasabah berlangganan deposito berjangka, sedangkan \(Y_i = 0\) berarti tidak berlangganan.
5.2 Sumber Data
Dataset diambil dari UCI Machine Learning Repository, yaitu dataset Bank Marketing. Data ini berkaitan dengan kampanye pemasaran langsung sebuah institusi perbankan di Portugal. Tujuan klasifikasinya adalah memprediksi apakah klien akan berlangganan deposito berjangka atau tidak.
Pada halaman UCI, dataset ini memiliki 45.211 observasi dan 16 variabel prediktor pada versi bank-full.csv. Dataset dapat diakses melalui tautan berikut:
https://archive.ics.uci.edu/dataset/222/bank+marketing
required_packages <- c("dplyr", "ggplot2", "knitr", "scales")
missing_packages <- required_packages[
!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages) > 0) {
stop(
"Paket berikut belum tersedia: ",
paste(missing_packages, collapse = ", "),
". Silakan install terlebih dahulu."
)
}
invisible(lapply(required_packages, library, character.only = TRUE))## Warning: package 'dplyr' was built under R version 4.4.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Warning: package 'ggplot2' was built under R version 4.4.3## Warning: package 'knitr' was built under R version 4.4.3## Warning: package 'scales' was built under R version 4.4.3uci_zip_url <- "https://archive.ics.uci.edu/static/public/222/bank+marketing.zip"
default_local_zip <- file.path(
"C:", "Users", "ASUS", "Documents", "CHEL", "Semester 4",
"Analisis Data Kategori", "Data Rpubs", "bank+marketing.zip"
)
# Jika zip sudah tersedia lokal, R akan memakainya terlebih dahulu.
# Jika ingin memakai lokasi lain, isi environment variable BANK_MARKETING_ZIP.
env_local_zip <- Sys.getenv("BANK_MARKETING_ZIP")
if (nzchar(env_local_zip) && file.exists(env_local_zip)) {
zip_file <- env_local_zip
} else if (file.exists(default_local_zip)) {
zip_file <- default_local_zip
} else {
zip_file <- file.path(tempdir(), "bank_marketing.zip")
download.file(uci_zip_url, zip_file, mode = "wb", quiet = TRUE)
}
zip_entries <- unzip(zip_file, list = TRUE)$Name
if ("bank/bank-full.csv" %in% zip_entries) {
bank_csv_connection <- unz(zip_file, "bank/bank-full.csv")
} else if ("bank-full.csv" %in% zip_entries) {
bank_csv_connection <- unz(zip_file, "bank-full.csv")
} else if ("bank.zip" %in% zip_entries) {
nested_dir <- file.path(tempdir(), "bank_marketing_nested")
dir.create(nested_dir, showWarnings = FALSE, recursive = TRUE)
unzip(zip_file, files = "bank.zip", exdir = nested_dir, overwrite = TRUE)
bank_csv_connection <- unz(file.path(nested_dir, "bank.zip"), "bank-full.csv")
} else {
stop("File bank-full.csv tidak ditemukan di dalam zip Bank Marketing.")
}
raw_bank <- read.csv(
bank_csv_connection,
sep = ";",
stringsAsFactors = FALSE
)
ringkasan_bank <- data.frame(
Keterangan = c("Jumlah observasi", "Jumlah variabel"),
Nilai = c(nrow(raw_bank), ncol(raw_bank))
)
knitr::kable(
ringkasan_bank,
caption = "Ukuran dataset Bank Marketing"
)| Keterangan | Nilai |
|---|---|
| Jumlah observasi | 45211 |
| Jumlah variabel | 17 |
Interpretasi output
Dataset bank-full.csv berisi 45211 observasi dan 17 variabel. Variabel respon adalah y, sedangkan variabel lainnya berperan sebagai prediktor karakteristik nasabah, riwayat kontak pemasaran, dan informasi kampanye sebelumnya.
5.3 Persiapan Data
Pada tahap persiapan data, variabel respon y diubah menjadi variabel biner bernama subscribe. Nilai 1 menunjukkan nasabah berlangganan deposito berjangka, sedangkan nilai 0 menunjukkan tidak berlangganan.
Beberapa variabel kategorik juga diubah menjadi faktor agar dapat digunakan secara tepat dalam model regresi logistik.
bank_data <- raw_bank %>%
mutate(
subscribe = ifelse(y == "yes", 1, 0),
subscribe_status = factor(
subscribe,
levels = c(0, 1),
labels = c("Tidak berlangganan", "Berlangganan")
),
job = factor(job),
marital = factor(marital),
education = factor(education),
default = factor(default),
housing = factor(housing),
loan = factor(loan),
contact = factor(contact),
poutcome = factor(poutcome)
) %>%
dplyr::select(
subscribe, subscribe_status, age, job, marital, education, default,
balance, housing, loan, contact, campaign, pdays, previous, poutcome
)
class_summary_bank <- bank_data %>%
count(subscribe_status, name = "Jumlah") %>%
mutate(Proporsi = scales::percent(Jumlah / sum(Jumlah), accuracy = 0.1)) %>%
rename(`Status deposito` = subscribe_status)
knitr::kable(
class_summary_bank,
caption = "Distribusi variabel respon"
)| Status deposito | Jumlah | Proporsi |
|---|---|---|
| Tidak berlangganan | 39922 | 88.3% |
| Berlangganan | 5289 | 11.7% |
Interpretasi output
Distribusi variabel respon menunjukkan perbandingan antara nasabah yang berlangganan dan tidak berlangganan deposito berjangka. Pada data pemasaran bank, jumlah nasabah yang tidak berlangganan umumnya jauh lebih besar dibandingkan nasabah yang berlangganan. Hal ini menunjukkan adanya ketidakseimbangan kelas yang perlu diperhatikan ketika mengevaluasi model.
5.4 Eksplorasi Singkat
ggplot(bank_data, aes(x = subscribe_status, fill = subscribe_status)) +
geom_bar(width = 0.62, color = "white", linewidth = 0.8) +
geom_text(
stat = "count",
aes(label = after_stat(count)),
vjust = -0.4,
fontface = "bold"
) +
scale_fill_manual(
values = c(
"Tidak berlangganan" = "#d63384",
"Berlangganan" = "#f783ac"
)
) +
labs(
title = "Distribusi Status Berlangganan Deposito",
x = NULL,
y = "Jumlah nasabah"
) +
theme_minimal(base_size = 12) +
theme(legend.position = "none")
Interpretasi output
Grafik batang memperlihatkan bahwa kategori Tidak berlangganan mendominasi data. Kondisi ini penting karena model yang hanya mengejar akurasi dapat tampak baik walaupun kurang mampu mendeteksi nasabah yang benar-benar berlangganan.
numeric_summary_bank <- bank_data %>%
group_by(subscribe_status) %>%
summarise(
Jumlah = n(),
`Rata-rata usia` = mean(age),
`Median saldo tahunan` = median(balance),
`Rata-rata jumlah kontak kampanye` = mean(campaign),
`Rata-rata kontak sebelumnya` = mean(previous),
.groups = "drop"
) %>%
rename(`Status deposito` = subscribe_status) %>%
mutate(across(where(is.numeric), round, 2))## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(where(is.numeric), round, 2)`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))knitr::kable(
numeric_summary_bank,
caption = "Ringkasan variabel numerik menurut status deposito"
)| Status deposito | Jumlah | Rata-rata usia | Median saldo tahunan | Rata-rata jumlah kontak kampanye | Rata-rata kontak sebelumnya |
|---|---|---|---|---|---|
| Tidak berlangganan | 39922 | 40.84 | 417 | 2.85 | 0.50 |
| Berlangganan | 5289 | 41.67 | 733 | 2.14 | 1.17 |
Interpretasi output
Tabel ringkasan numerik membantu membandingkan profil nasabah yang berlangganan dan tidak berlangganan. Perbedaan rata-rata usia, saldo, jumlah kontak kampanye, dan kontak sebelumnya dapat menjadi indikasi awal bahwa variabel-variabel tersebut berhubungan dengan peluang keberhasilan kampanye.
ggplot(bank_data, aes(x = age, fill = subscribe_status)) +
geom_density(alpha = 0.55, color = "white", linewidth = 0.7) +
scale_fill_manual(
values = c(
"Tidak berlangganan" = "#d63384",
"Berlangganan" = "#ffe066"
)
) +
labs(
title = "Distribusi Usia Nasabah Menurut Status Deposito",
x = "Usia nasabah",
y = "Kepadatan",
fill = "Status deposito"
) +
theme_minimal(base_size = 12)
Interpretasi output
Grafik kepadatan memperlihatkan pola usia pada dua kelompok respon. Jika kurva nasabah yang berlangganan berbeda dari kurva nasabah yang tidak berlangganan, maka usia dapat menjadi salah satu prediktor yang membantu model membedakan kedua kategori.
5.5 Pembagian Data Training dan Testing
Data dibagi menjadi 80% training dan 20% testing secara stratified agar proporsi nasabah yang berlangganan dan tidak berlangganan tetap seimbang pada kedua subset.
set.seed(123)
stratified_split <- function(y, prop = 0.8) {
idx_by_class <- split(seq_along(y), y)
train_idx <- lapply(
idx_by_class,
function(idx) sample(idx, size = floor(length(idx) * prop))
)
unlist(train_idx, use.names = FALSE)
}
train_id_bank <- stratified_split(bank_data$subscribe, prop = 0.8)
train_bank <- bank_data[train_id_bank, ]
test_bank <- bank_data[-train_id_bank, ]
split_summary_bank <- bind_rows(
train_bank %>% count(subscribe_status) %>% mutate(data = "Training"),
test_bank %>% count(subscribe_status) %>% mutate(data = "Testing")
) %>%
group_by(data) %>%
mutate(proporsi = scales::percent(n / sum(n), accuracy = 0.1)) %>%
ungroup() %>%
dplyr::select(data, subscribe_status, n, proporsi) %>%
rename(
Data = data,
`Status deposito` = subscribe_status,
Jumlah = n,
Proporsi = proporsi
)
knitr::kable(
split_summary_bank,
caption = "Distribusi kelas pada data training dan testing"
)| Data | Status deposito | Jumlah | Proporsi |
|---|---|---|---|
| Training | Tidak berlangganan | 31937 | 88.3% |
| Training | Berlangganan | 4231 | 11.7% |
| Testing | Tidak berlangganan | 7985 | 88.3% |
| Testing | Berlangganan | 1058 | 11.7% |
Interpretasi output
Pembagian stratified menjaga agar komposisi kelas respon pada data training dan testing tetap menyerupai data asli. Data training digunakan untuk mengestimasi parameter model, sedangkan data testing digunakan untuk mengevaluasi kemampuan model pada data yang belum digunakan dalam proses pelatihan.
5.6 Model Regresi Logistik Biner
Variabel duration tidak dimasukkan ke dalam model karena durasi panggilan baru diketahui setelah kontak terjadi. Jika tujuan analisis adalah memprediksi sebelum kampanye selesai, memasukkan durasi dapat menyebabkan informasi yang terlalu optimistis.
Model yang digunakan adalah:
\[ \text{logit}(p_i) = \ln\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \cdots + \beta_kX_{ki}. \]
Bentuk peluang prediksinya adalah:
\[ \hat{p}_i = \frac{\exp(\hat{\eta}_i)}{1+\exp(\hat{\eta}_i)} = \frac{1}{1+\exp(-\hat{\eta}_i)}. \]
bank_fit <- glm(
subscribe ~ age + job + marital + education + default + balance +
housing + loan + contact + campaign + pdays + previous + poutcome,
data = train_bank,
family = binomial(link = "logit")
)
ringkasan_model_bank <- data.frame(
Keterangan = c(
"Jumlah observasi training",
"Null deviance",
"Residual deviance",
"Derajat bebas residual",
"AIC"
),
Nilai = c(
nobs(bank_fit),
round(bank_fit$null.deviance, 3),
round(bank_fit$deviance, 3),
bank_fit$df.residual,
round(AIC(bank_fit), 3)
)
)
knitr::kable(
ringkasan_model_bank,
caption = "Ringkasan kecocokan model regresi logistik"
)| Keterangan | Nilai |
|---|---|
| Jumlah observasi training | 36168.00 |
| Null deviance | 26103.76 |
| Residual deviance | 22433.01 |
| Derajat bebas residual | 36138.00 |
| AIC | 22493.01 |
Interpretasi output
Jika residual deviance lebih kecil daripada null deviance, maka model dengan prediktor memberikan kecocokan yang lebih baik dibandingkan model tanpa prediktor. Nilai AIC dapat digunakan untuk membandingkan beberapa model; semakin kecil AIC, semakin baik keseimbangan antara kecocokan model dan kompleksitas model.
5.7 Odds Ratio dan Uji Parameter
Koefisien regresi logistik berada pada skala log-odds. Agar lebih mudah dibaca, koefisien diubah menjadi odds ratio dengan rumus:
\[ OR_j = \exp(\beta_j). \]
Jika \(OR_j > 1\), maka prediktor tersebut meningkatkan odds nasabah berlangganan deposito. Jika \(OR_j < 1\), maka prediktor tersebut menurunkan odds nasabah berlangganan deposito. Untuk variabel kategorik, interpretasi dilakukan relatif terhadap kategori referensi.
coef_raw_bank <- as.data.frame(summary(bank_fit)$coefficients)
coef_raw_bank$term <- rownames(coef_raw_bank)
rownames(coef_raw_bank) <- NULL
coef_table_bank <- coef_raw_bank %>%
filter(term != "(Intercept)") %>%
mutate(
odds_ratio = exp(Estimate),
ci_low = exp(Estimate - 1.96 * `Std. Error`),
ci_high = exp(Estimate + 1.96 * `Std. Error`)
) %>%
arrange(`Pr(>|z|)`) %>%
transmute(
`Variabel/level` = term,
`Odds ratio` = round(odds_ratio, 3),
`Interval kepercayaan 95 persen` = paste0(
round(ci_low, 3), " - ", round(ci_high, 3)
),
`p-value` = signif(`Pr(>|z|)`, 3)
)
knitr::kable(
head(coef_table_bank, 20),
caption = "Dua puluh koefisien paling menonjol berdasarkan p-value"
)| Variabel/level | Odds ratio | Interval kepercayaan 95 persen | p-value |
|---|---|---|---|
| poutcomesuccess | 9.359 | 7.979 - 10.977 | 0.00e+00 |
| contactunknown | 0.371 | 0.333 - 0.415 | 0.00e+00 |
| housingyes | 0.568 | 0.526 - 0.613 | 0.00e+00 |
| campaign | 0.893 | 0.876 - 0.91 | 0.00e+00 |
| loanyes | 0.617 | 0.551 - 0.691 | 0.00e+00 |
| educationtertiary | 1.391 | 1.206 - 1.604 | 5.50e-06 |
| jobstudent | 1.622 | 1.313 - 2.004 | 7.50e-06 |
| jobretired | 1.521 | 1.265 - 1.83 | 8.60e-06 |
| balance | 1.000 | 1 - 1 | 1.45e-04 |
| jobentrepreneur | 0.659 | 0.514 - 0.845 | 9.97e-04 |
| maritalmarried | 0.839 | 0.749 - 0.939 | 2.21e-03 |
| poutcomeother | 1.294 | 1.089 - 1.537 | 3.40e-03 |
| jobblue-collar | 0.823 | 0.715 - 0.947 | 6.56e-03 |
| contacttelephone | 0.829 | 0.721 - 0.954 | 8.64e-03 |
| maritalsingle | 1.183 | 1.041 - 1.345 | 1.01e-02 |
| jobhousemaid | 0.731 | 0.567 - 0.943 | 1.57e-02 |
| educationunknown | 1.267 | 1.039 - 1.544 | 1.93e-02 |
| jobmanagement | 0.847 | 0.736 - 0.974 | 2.03e-02 |
| jobtechnician | 0.860 | 0.754 - 0.982 | 2.53e-02 |
| jobservices | 0.856 | 0.729 - 1.004 | 5.67e-02 |
Interpretasi output
Tabel odds ratio menampilkan prediktor yang paling menonjol berdasarkan p-value. Prediktor dengan p-value kecil menunjukkan bukti statistik yang lebih kuat terhadap hubungan dengan peluang nasabah berlangganan deposito. Namun, interpretasi tetap perlu memperhatikan konteks bisnis dan kategori referensi masing-masing variabel.
5.8 Prediksi dan Evaluasi Klasifikasi
Setelah model menghasilkan peluang prediksi \(\hat{p}_i\), kelas prediksi ditentukan menggunakan threshold \(c\):
\[ \hat{Y}_i = \begin{cases} 1, & \text{jika } \hat{p}_i \ge c,\\ 0, & \text{jika } \hat{p}_i < c. \end{cases} \]
Jika \(c = 0.50\), maka nasabah diklasifikasikan berlangganan deposito ketika peluang prediksinya minimal 50%.
p_train_bank <- predict(bank_fit, newdata = train_bank, type = "response")
p_test_bank <- predict(bank_fit, newdata = test_bank, type = "response")
prediction_preview_bank <- head(
data.frame(
`Status aktual` = test_bank$subscribe_status,
`Peluang prediksi berlangganan` = round(p_test_bank, 4),
check.names = FALSE
),
6
)
knitr::kable(
prediction_preview_bank,
caption = "Contoh peluang prediksi pada data testing"
)| Status aktual | Peluang prediksi berlangganan | |
|---|---|---|
| 12 | Tidak berlangganan | 0.0525 |
| 23 | Tidak berlangganan | 0.0252 |
| 34 | Tidak berlangganan | 0.0331 |
| 36 | Tidak berlangganan | 0.0409 |
| 44 | Tidak berlangganan | 0.0595 |
| 61 | Tidak berlangganan | 0.0478 |
Interpretasi output
Tabel ini menampilkan contoh peluang prediksi pada data testing. Semakin besar peluang prediksi, semakin besar kemungkinan nasabah diklasifikasikan sebagai nasabah yang akan berlangganan deposito berjangka.
safe_div_bank <- function(num, den) {
ifelse(den == 0, NA_real_, num / den)
}
classification_metrics_bank <- function(actual, prob, threshold = 0.5) {
pred <- as.integer(prob >= threshold)
tp <- sum(pred == 1 & actual == 1)
tn <- sum(pred == 0 & actual == 0)
fp <- sum(pred == 1 & actual == 0)
fn <- sum(pred == 0 & actual == 1)
sensitivity <- safe_div_bank(tp, tp + fn)
specificity <- safe_div_bank(tn, tn + fp)
precision <- safe_div_bank(tp, tp + fp)
negative_predictive_value <- safe_div_bank(tn, tn + fn)
accuracy <- safe_div_bank(tp + tn, tp + tn + fp + fn)
data.frame(
threshold = threshold,
accuracy = accuracy,
error_rate = 1 - accuracy,
sensitivity = sensitivity,
specificity = specificity,
precision = precision,
negative_predictive_value = negative_predictive_value,
f1_score = safe_div_bank(2 * precision * sensitivity, precision + sensitivity),
balanced_accuracy = (sensitivity + specificity) / 2,
false_positive_rate = 1 - specificity,
false_negative_rate = 1 - sensitivity
)
}
format_metrics_bank <- function(x) {
x %>%
transmute(
Threshold = threshold,
Akurasi = accuracy,
`Error rate` = error_rate,
Sensitivity = sensitivity,
Specificity = specificity,
Presisi = precision,
NPV = negative_predictive_value,
`F1-score` = f1_score,
`Balanced accuracy` = balanced_accuracy,
FPR = false_positive_rate,
FNR = false_negative_rate
)
}
confusion_matrix_bank <- function(actual, prob, threshold = 0.5) {
pred <- as.integer(prob >= threshold)
tab <- table(
Aktual = factor(
actual,
levels = c(1, 0),
labels = c("Aktual Berlangganan", "Aktual Tidak berlangganan")
),
Prediksi = factor(
pred,
levels = c(1, 0),
labels = c("Prediksi Berlangganan", "Prediksi Tidak berlangganan")
)
)
addmargins(tab)
}confusion_default_bank <- confusion_matrix_bank(
test_bank$subscribe,
p_test_bank,
threshold = 0.5
)
metrics_default_bank <- classification_metrics_bank(
test_bank$subscribe,
p_test_bank,
threshold = 0.5
) %>%
format_metrics_bank() %>%
mutate(across(where(is.numeric), round, 3))
knitr::kable(
confusion_default_bank,
caption = "Confusion matrix testing pada threshold 0.50"
)| Prediksi Berlangganan | Prediksi Tidak berlangganan | Sum | |
|---|---|---|---|
| Aktual Berlangganan | 165 | 893 | 1058 |
| Aktual Tidak berlangganan | 86 | 7899 | 7985 |
| Sum | 251 | 8792 | 9043 |
| Threshold | Akurasi | Error rate | Sensitivity | Specificity | Presisi | NPV | F1-score | Balanced accuracy | FPR | FNR |
|---|---|---|---|---|---|---|---|---|---|---|
| 0.5 | 0.892 | 0.108 | 0.156 | 0.989 | 0.657 | 0.898 | 0.252 | 0.573 | 0.011 | 0.844 |
Interpretasi output
Confusion matrix menunjukkan jumlah prediksi benar dan salah pada data testing. Pada konteks pemasaran bank, sensitivity penting karena menunjukkan kemampuan model menemukan nasabah yang benar-benar berlangganan. Specificity menunjukkan kemampuan model mengenali nasabah yang tidak berlangganan.
5.9 Kurva ROC dan Threshold Optimal
Kurva ROC menggambarkan hubungan antara:
- True positive rate atau sensitivity
- False positive rate, yaitu \(1 - specificity\)
Nilai AUC yang semakin mendekati 1 menunjukkan kemampuan klasifikasi yang semakin baik, sedangkan AUC sekitar 0.5 menunjukkan kemampuan model mendekati tebakan acak.
Threshold optimal dipilih menggunakan indeks Youden:
\[ J = sensitivity + specificity - 1. \]
Threshold dengan nilai \(J\) terbesar dianggap memberi keseimbangan terbaik antara sensitivity dan specificity.
roc_points_bank <- function(actual, prob) {
thresholds <- c(Inf, sort(unique(prob), decreasing = TRUE), -Inf)
out <- lapply(thresholds, function(th) {
pred <- as.integer(prob >= th)
tp <- sum(pred == 1 & actual == 1)
tn <- sum(pred == 0 & actual == 0)
fp <- sum(pred == 1 & actual == 0)
fn <- sum(pred == 0 & actual == 1)
sensitivity <- safe_div_bank(tp, tp + fn)
specificity <- safe_div_bank(tn, tn + fp)
data.frame(
threshold = th,
sensitivity = sensitivity,
specificity = specificity,
fpr = 1 - specificity,
youden = sensitivity + specificity - 1
)
})
bind_rows(out)
}
auc_value_bank <- function(roc_df) {
roc_sorted <- roc_df %>%
arrange(fpr, sensitivity)
sum(
diff(roc_sorted$fpr) *
(head(roc_sorted$sensitivity, -1) + tail(roc_sorted$sensitivity, -1)) / 2
)
}
roc_train_bank <- roc_points_bank(train_bank$subscribe, p_train_bank) %>%
mutate(data = "Training")
roc_test_bank <- roc_points_bank(test_bank$subscribe, p_test_bank) %>%
mutate(data = "Testing")
auc_train_bank <- auc_value_bank(roc_train_bank)
auc_test_bank <- auc_value_bank(roc_test_bank)
optimal_train_bank <- roc_train_bank %>%
filter(is.finite(threshold)) %>%
arrange(desc(youden), desc(sensitivity)) %>%
slice(1)
threshold_opt_bank <- optimal_train_bank$threshold
test_at_opt_bank <- roc_test_bank %>%
filter(is.finite(threshold)) %>%
mutate(distance = abs(threshold - threshold_opt_bank)) %>%
arrange(distance) %>%
slice(1)roc_plot_bank <- bind_rows(roc_train_bank, roc_test_bank)
ggplot(roc_plot_bank, aes(x = fpr, y = sensitivity, color = data)) +
geom_path(linewidth = 1.1) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "#6c757d") +
geom_point(
data = optimal_train_bank,
aes(x = fpr, y = sensitivity),
inherit.aes = FALSE,
color = "#d63384",
fill = "#ffe066",
shape = 21,
size = 4,
stroke = 1.2
) +
coord_equal() +
scale_color_manual(values = c("Training" = "#d63384", "Testing" = "#f783ac")) +
labs(
title = "Kurva ROC Model Regresi Logistik Bank Marketing",
subtitle = paste0(
"AUC training = ", round(auc_train_bank, 3),
"; AUC testing = ", round(auc_test_bank, 3),
"; threshold optimal = ", round(threshold_opt_bank, 3)
),
x = "False positive rate (1 - specificity)",
y = "Sensitivity / true positive rate",
color = "Data"
) +
theme_minimal(base_size = 12) +
theme(legend.position = "bottom")
Interpretasi output
Kurva ROC yang semakin jauh dari garis diagonal menunjukkan performa klasifikasi yang semakin baik. Jika AUC testing cukup tinggi dan tidak jauh dari AUC training, maka model relatif stabil ketika digunakan pada data baru. Titik pada kurva menunjukkan threshold optimal berdasarkan data training.
5.10 Evaluasi dengan Threshold Optimal
Threshold optimal dari ROC training diterapkan pada data testing. Cara ini lebih tepat dibandingkan memilih threshold langsung dari testing, karena data testing sebaiknya diperlakukan sebagai data baru.
metrics_compare_bank <- bind_rows(
classification_metrics_bank(test_bank$subscribe, p_test_bank, threshold = 0.5) %>%
mutate(aturan = "Threshold 0.50"),
classification_metrics_bank(test_bank$subscribe, p_test_bank, threshold = threshold_opt_bank) %>%
mutate(aturan = "Threshold optimal ROC")
) %>%
dplyr::select(aturan, everything())
metrics_compare_bank_formatted <- metrics_compare_bank %>%
format_metrics_bank() %>%
bind_cols(
`Aturan klasifikasi` = c("Threshold 0.50", "Threshold optimal ROC"),
.
) %>%
dplyr::select(`Aturan klasifikasi`, everything()) %>%
mutate(across(where(is.numeric), round, 3))
confusion_opt_bank <- confusion_matrix_bank(
test_bank$subscribe,
p_test_bank,
threshold = threshold_opt_bank
)
knitr::kable(
metrics_compare_bank_formatted,
caption = "Perbandingan metrik testing"
)| Aturan klasifikasi | Threshold | Akurasi | Error rate | Sensitivity | Specificity | Presisi | NPV | F1-score | Balanced accuracy | FPR | FNR |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Threshold 0.50 | 0.500 | 0.892 | 0.108 | 0.156 | 0.989 | 0.657 | 0.898 | 0.252 | 0.573 | 0.011 | 0.844 |
| Threshold optimal ROC | 0.117 | 0.693 | 0.307 | 0.681 | 0.695 | 0.228 | 0.943 | 0.342 | 0.688 | 0.305 | 0.319 |
| Prediksi Berlangganan | Prediksi Tidak berlangganan | Sum | |
|---|---|---|---|
| Aktual Berlangganan | 720 | 338 | 1058 |
| Aktual Tidak berlangganan | 2434 | 5551 | 7985 |
| Sum | 3154 | 5889 | 9043 |
Interpretasi output
Perbandingan metrik menunjukkan dampak perubahan threshold. Threshold yang lebih rendah biasanya meningkatkan sensitivity karena lebih banyak nasabah diklasifikasikan sebagai calon pelanggan deposito, tetapi dapat menurunkan specificity karena lebih banyak nasabah yang sebenarnya tidak berlangganan ikut terklasifikasi positif.
test_prob_plot_bank <- test_bank %>%
mutate(peluang_berlangganan = p_test_bank)
ggplot(test_prob_plot_bank, aes(x = peluang_berlangganan, fill = subscribe_status)) +
geom_density(alpha = 0.55, color = "white", linewidth = 0.7) +
geom_vline(
xintercept = threshold_opt_bank,
color = "#d63384",
linewidth = 1.2,
linetype = "dashed"
) +
annotate(
"label",
x = threshold_opt_bank,
y = Inf,
label = paste0("threshold = ", round(threshold_opt_bank, 3)),
vjust = 1.4,
fill = "#fff9e6",
color = "#d63384",
label.size = 0
) +
scale_fill_manual(
values = c(
"Tidak berlangganan" = "#d63384",
"Berlangganan" = "#ffe066"
)
) +
labs(
title = "Distribusi Peluang Prediksi pada Data Testing",
x = "Peluang prediksi berlangganan deposito",
y = "Kepadatan",
fill = "Status aktual"
) +
theme_minimal(base_size = 12)
Interpretasi output
Grafik kepadatan memperlihatkan sebaran peluang prediksi pada dua kelompok aktual. Semakin terpisah dua kurva, semakin baik model membedakan nasabah yang berlangganan dan tidak berlangganan. Area tumpang tindih menunjukkan bagian data yang sulit diklasifikasikan.
5.11 Kesimpulan
Regresi logistik biner dapat digunakan untuk memprediksi peluang nasabah berlangganan deposito berjangka berdasarkan karakteristik nasabah dan informasi kampanye pemasaran. Model menghasilkan peluang prediksi yang kemudian dapat diubah menjadi kelas menggunakan threshold tertentu.
Pada data Bank Marketing, evaluasi tidak cukup hanya menggunakan akurasi karena kelas respon tidak seimbang. Oleh karena itu, metrik seperti sensitivity, specificity, F1-score, balanced accuracy, dan AUC perlu diperhatikan agar interpretasi performa model lebih lengkap.
Secara substantif, threshold klasifikasi dapat disesuaikan dengan tujuan bank. Jika bank ingin menangkap lebih banyak calon nasabah potensial, threshold dapat dibuat lebih rendah. Namun, konsekuensinya adalah lebih banyak nasabah yang sebenarnya tidak berlangganan dapat ikut diklasifikasikan sebagai calon pelanggan potensial.
6 Regresi Logistik Multinomial: Iris
6.1 Ringkasan Masalah
Regresi logistik multinomial digunakan ketika variabel respon bersifat kategorik nominal dengan lebih dari dua kategori. Pada bab ini digunakan dataset Iris dari UCI Machine Learning Repository untuk mengklasifikasikan spesies bunga iris berdasarkan ukuran sepal dan petal.
Variabel respon adalah species yang memiliki tiga kategori:
- Iris-setosa
- Iris-versicolor
- Iris-virginica
Catatan interpretasi: dataset Iris sebenarnya kurang kuat untuk interpretasi substantif seperti penelitian sosial atau ekonomi, karena variabel sepal dan petal lebih tepat disebut indikator morfologi untuk membedakan spesies, bukan prediktor kausal. Namun, dataset ini tetap dapat digunakan untuk mencoba regresi logistik multinomial karena responnya memiliki tiga kategori nominal dan hubungan antarvariabelnya masih dapat diinterpretasikan secara klasifikasi.
6.2 Formulasi Model
Misalkan variabel respon \(Y_i\) memiliki \(J\) kategori. Pada model baseline-category logit, salah satu kategori dipilih sebagai kategori referensi. Jika kategori referensi adalah kategori \(J\), maka untuk kategori \(j = 1,2,\ldots,J-1\):
\[ \log\left(\frac{P(Y_i=j)}{P(Y_i=J)}\right) = \beta_{0j} + \beta_{1j}X_{1i} + \beta_{2j}X_{2i} + \cdots + \beta_{kj}X_{ki}. \]
Koefisien pada model multinomial menunjukkan perubahan log-odds suatu kategori dibandingkan kategori referensi.
6.3 Sumber dan Persiapan Data
Dataset Iris berisi 150 observasi, 4 fitur numerik, dan 3 kelas spesies iris. Dataset asli dapat diakses melalui UCI Machine Learning Repository:
https://archive.ics.uci.edu/dataset/53/iris
required_packages_mpo <- c("dplyr", "ggplot2", "knitr", "scales", "nnet", "MASS", "tidyr")
missing_packages_mpo <- required_packages_mpo[
!vapply(required_packages_mpo, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages_mpo) > 0) {
stop(
"Paket berikut belum tersedia: ",
paste(missing_packages_mpo, collapse = ", "),
". Silakan install terlebih dahulu."
)
}
invisible(lapply(required_packages_mpo, library, character.only = TRUE))##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Warning: package 'tidyr' was built under R version 4.4.3iris_zip <- file.path(
"C:", "Users", "ASUS", "Documents", "CHEL", "Semester 4",
"Analisis Data Kategori", "Data Rpubs", "iris (3).zip"
)
if (!file.exists(iris_zip)) {
stop("File iris (3).zip tidak ditemukan pada path lokal.")
}
iris_raw <- read.csv(
unz(iris_zip, "iris.data"),
header = FALSE,
col.names = c(
"sepal_length", "sepal_width",
"petal_length", "petal_width",
"species"
)
)
iris_data <- iris_raw %>%
filter(species != "") %>%
mutate(species = factor(species))
iris_summary <- data.frame(
Keterangan = c("Jumlah observasi", "Jumlah prediktor", "Jumlah kategori respon"),
Nilai = c(nrow(iris_data), 4, nlevels(iris_data$species))
)
knitr::kable(
iris_summary,
caption = "Ringkasan dataset Iris"
)| Keterangan | Nilai |
|---|---|
| Jumlah observasi | 150 |
| Jumlah prediktor | 4 |
| Jumlah kategori respon | 3 |
Interpretasi output
Data Iris memiliki 150 observasi dengan empat prediktor numerik. Variabel respon memiliki 3 kategori, sehingga cocok digunakan untuk contoh regresi logistik multinomial.
6.4 Distribusi Kategori Respons
iris_class_summary <- iris_data %>%
count(species, name = "Jumlah") %>%
mutate(Proporsi = scales::percent(Jumlah / sum(Jumlah), accuracy = 0.1)) %>%
rename(`Spesies iris` = species)
knitr::kable(
iris_class_summary,
caption = "Distribusi spesies iris"
)| Spesies iris | Jumlah | Proporsi |
|---|---|---|
| Iris-setosa | 50 | 33.3% |
| Iris-versicolor | 50 | 33.3% |
| Iris-virginica | 50 | 33.3% |
Interpretasi output
Distribusi spesies pada dataset Iris seimbang. Setiap spesies memiliki jumlah observasi yang sama, sehingga evaluasi klasifikasi tidak terlalu dipengaruhi oleh ketidakseimbangan kelas.
ggplot(iris_data, aes(x = petal_length, y = petal_width, color = species)) +
geom_point(size = 2.4, alpha = 0.85) +
scale_color_manual(values = c("#d63384", "#f783ac", "#ffe066")) +
labs(
title = "Sebaran Panjang dan Lebar Petal Menurut Spesies",
x = "Panjang petal",
y = "Lebar petal",
color = "Spesies"
) +
theme_minimal(base_size = 12)
Interpretasi output
Plot memperlihatkan bahwa ukuran petal dapat membedakan spesies iris dengan cukup jelas. Hal ini menjelaskan mengapa dataset Iris sering digunakan sebagai contoh klasifikasi, meskipun interpretasi prediktornya lebih bersifat identifikasi biologis daripada hubungan kausal.
6.5 Estimasi Model Multinomial
set.seed(123)
iris_train_id <- unlist(
lapply(split(seq_len(nrow(iris_data)), iris_data$species), function(idx) {
sample(idx, size = floor(0.8 * length(idx)))
}),
use.names = FALSE
)
iris_train <- iris_data[iris_train_id, ]
iris_test <- iris_data[-iris_train_id, ]
iris_multinom_fit <- nnet::multinom(
species ~ sepal_length + sepal_width + petal_length + petal_width,
data = iris_train,
trace = FALSE
)
iris_model_summary <- data.frame(
Keterangan = c("Jumlah observasi training", "AIC", "Deviance residual"),
Nilai = c(
nrow(iris_train),
round(AIC(iris_multinom_fit), 3),
round(deviance(iris_multinom_fit), 3)
)
)
knitr::kable(
iris_model_summary,
caption = "Ringkasan model regresi logistik multinomial"
)| Keterangan | Nilai |
|---|---|
| Jumlah observasi training | 120.000 |
| AIC | 23.262 |
| Deviance residual | 3.262 |
Interpretasi output
Model multinomial mengestimasi log-odds setiap spesies relatif terhadap kategori referensi. Nilai AIC digunakan untuk melihat keseimbangan antara kecocokan model dan kompleksitas model.
6.6 Likelihood, Deviance, dan Goodness-of-Fit Multinomial
Model multinomial diestimasi dengan metode maximum likelihood. Untuk setiap observasi, model menghasilkan probabilitas prediksi pada seluruh kategori respons. Secara umum likelihood model multinomial dapat ditulis sebagai:
\[ L(\boldsymbol{\beta}) = \prod_{i=1}^{n}\prod_{j=1}^{J} \pi_{ij}^{y_{ij}}, \]
dengan \(\pi_{ij}=P(Y_i=j \mid X_i)\), dan \(y_{ij}=1\) jika observasi ke-\(i\) berada pada kategori \(j\).
Goodness-of-fit dapat dilihat dari log-likelihood, deviance, AIC, dan likelihood ratio test terhadap model null.
iris_multinom_null <- nnet::multinom(
species ~ 1,
data = iris_train,
trace = FALSE
)
ll_iris_full <- as.numeric(logLik(iris_multinom_fit))
ll_iris_null <- as.numeric(logLik(iris_multinom_null))
lr_iris <- 2 * (ll_iris_full - ll_iris_null)
df_iris <- attr(logLik(iris_multinom_fit), "df") -
attr(logLik(iris_multinom_null), "df")
iris_gof_table <- data.frame(
Ukuran = c(
"Log-likelihood model null",
"Log-likelihood model penuh",
"Likelihood ratio statistic",
"Derajat bebas LR",
"p-value LR",
"McFadden pseudo R2",
"AIC model penuh"
),
Nilai = c(
round(ll_iris_null, 3),
round(ll_iris_full, 3),
round(lr_iris, 3),
df_iris,
signif(pchisq(lr_iris, df = df_iris, lower.tail = FALSE), 3),
round(1 - ll_iris_full / ll_iris_null, 3),
round(AIC(iris_multinom_fit), 3)
)
)
knitr::kable(
iris_gof_table,
caption = "Goodness-of-fit model multinomial"
)| Ukuran | Nilai |
|---|---|
| Log-likelihood model null | -131.833 |
| Log-likelihood model penuh | -1.631 |
| Likelihood ratio statistic | 260.405 |
| Derajat bebas LR | 8.000 |
| p-value LR | 0.000 |
| McFadden pseudo R2 | 0.988 |
| AIC model penuh | 23.262 |
Interpretasi output
Likelihood ratio test membandingkan model penuh dengan model null yang hanya memiliki intercept. Jika p-value kecil, maka prediktor secara bersama-sama memberikan peningkatan kecocokan model dibandingkan model tanpa prediktor. McFadden pseudo-\(R^2\) bukan \(R^2\) regresi linear, tetapi memberi indikasi peningkatan log-likelihood relatif terhadap model null.
6.7 Koefisien dan Odds Ratio
iris_coef <- summary(iris_multinom_fit)$coefficients
iris_se <- summary(iris_multinom_fit)$standard.errors
iris_coef_table <- as.data.frame(as.table(iris_coef)) %>%
rename(Kategori = Var1, Variabel = Var2, Koefisien = Freq) %>%
mutate(
`Std. Error` = as.vector(iris_se),
`z-value` = Koefisien / `Std. Error`,
`p-value` = 2 * (1 - pnorm(abs(`z-value`))),
`Odds ratio` = exp(Koefisien)
) %>%
arrange(Kategori, `p-value`) %>%
mutate(
Koefisien = round(Koefisien, 3),
`Std. Error` = round(`Std. Error`, 3),
`z-value` = round(`z-value`, 3),
`Odds ratio` = round(`Odds ratio`, 3),
`p-value` = signif(`p-value`, 3)
)
knitr::kable(
iris_coef_table,
caption = "Koefisien model multinomial dan odds ratio"
)| Kategori | Variabel | Koefisien | Std. Error | z-value | p-value | Odds ratio |
|---|---|---|---|---|---|---|
| Iris-versicolor | sepal_width | -125.078 | 7.603 | -16.451 | 0.00e+00 | 0.000000e+00 |
| Iris-versicolor | petal_length | 191.716 | 17.403 | 11.016 | 0.00e+00 | 1.824894e+83 |
| Iris-versicolor | sepal_length | -50.146 | 8.376 | -5.987 | 0.00e+00 | 0.000000e+00 |
| Iris-versicolor | petal_width | 76.718 | 38.259 | 2.005 | 4.49e-02 | 2.081182e+33 |
| Iris-versicolor | (Intercept) | 44.886 | 88.179 | 0.509 | 6.11e-01 | 3.116534e+19 |
| Iris-virginica | sepal_width | -142.139 | 7.603 | -18.695 | 0.00e+00 | 0.000000e+00 |
| Iris-virginica | petal_length | 224.937 | 17.403 | 12.925 | 0.00e+00 | 4.884975e+97 |
| Iris-virginica | sepal_length | -62.248 | 8.376 | -7.432 | 0.00e+00 | 0.000000e+00 |
| Iris-virginica | petal_width | 155.843 | 38.259 | 4.073 | 4.63e-05 | 4.804688e+67 |
| Iris-virginica | (Intercept) | -132.545 | 88.179 | -1.503 | 1.33e-01 | 0.000000e+00 |
Interpretasi output
Odds ratio di atas 1 menunjukkan bahwa kenaikan prediktor meningkatkan odds suatu spesies dibandingkan kategori referensi. Odds ratio di bawah 1 menunjukkan penurunan odds relatif. Karena ini dataset morfologi, interpretasi yang paling tepat adalah bahwa ukuran sepal dan petal membantu membedakan pola karakteristik antarspesies.
6.8 Prediksi dan Confusion Matrix
iris_prob_test <- predict(iris_multinom_fit, newdata = iris_test, type = "probs")
iris_pred_test <- predict(iris_multinom_fit, newdata = iris_test, type = "class")
iris_confusion <- table(
Aktual = iris_test$species,
Prediksi = iris_pred_test
)
iris_accuracy <- sum(diag(iris_confusion)) / sum(iris_confusion)
knitr::kable(
addmargins(iris_confusion),
caption = "Confusion matrix model multinomial pada data testing"
)| Iris-setosa | Iris-versicolor | Iris-virginica | Sum | |
|---|---|---|---|---|
| Iris-setosa | 10 | 0 | 0 | 10 |
| Iris-versicolor | 0 | 10 | 0 | 10 |
| Iris-virginica | 0 | 2 | 8 | 10 |
| Sum | 10 | 12 | 8 | 30 |
iris_accuracy_table <- data.frame(
Metrik = "Akurasi testing",
Nilai = round(iris_accuracy, 3)
)
knitr::kable(
iris_accuracy_table,
caption = "Akurasi model multinomial"
)| Metrik | Nilai |
|---|---|
| Akurasi testing | 0.933 |
Interpretasi output
Confusion matrix menunjukkan jumlah spesies yang diklasifikasikan dengan benar dan salah. Jika akurasi tinggi, model berhasil memanfaatkan ukuran sepal dan petal untuk membedakan spesies iris pada data testing.
iris_prob_plot <- as.data.frame(iris_prob_test)
iris_prob_plot$Aktual <- iris_test$species
iris_prob_plot$Prediksi <- iris_pred_test
iris_prob_plot$Nomor <- seq_len(nrow(iris_prob_plot))
iris_prob_long <- iris_prob_plot %>%
tidyr::pivot_longer(
cols = all_of(levels(iris_data$species)),
names_to = "Spesies",
values_to = "Probabilitas"
)ggplot(iris_prob_long, aes(x = Nomor, y = Probabilitas, color = Spesies)) +
geom_line(linewidth = 0.9) +
geom_point(size = 1.8) +
scale_color_manual(values = c("#d63384", "#f783ac", "#ffe066")) +
labs(
title = "Probabilitas Prediksi Multinomial pada Data Testing",
x = "Observasi testing",
y = "Probabilitas prediksi",
color = "Spesies"
) +
theme_minimal(base_size = 12)
Interpretasi output
Grafik probabilitas menunjukkan tingkat keyakinan model untuk setiap spesies. Observasi yang memiliki satu probabilitas sangat dominan lebih mudah diklasifikasikan dibandingkan observasi dengan probabilitas yang saling berdekatan.
6.9 Metrik Evaluasi Multikelas
Akurasi saja tidak selalu cukup, terutama jika kelas respons tidak seimbang. Untuk klasifikasi multikelas, evaluasi dapat dilengkapi dengan recall per kelas dan balanced accuracy.
iris_class_metrics <- data.frame(
Kelas = rownames(iris_confusion),
Recall = diag(iris_confusion) / rowSums(iris_confusion),
Precision = diag(iris_confusion) / colSums(iris_confusion)
) %>%
mutate(
Recall = round(Recall, 3),
Precision = round(Precision, 3)
)
iris_balanced_accuracy <- mean(iris_class_metrics$Recall, na.rm = TRUE)
knitr::kable(
iris_class_metrics,
caption = "Recall dan precision per kelas pada model multinomial"
)| Kelas | Recall | Precision | |
|---|---|---|---|
| Iris-setosa | Iris-setosa | 1.0 | 1.000 |
| Iris-versicolor | Iris-versicolor | 1.0 | 0.833 |
| Iris-virginica | Iris-virginica | 0.8 | 1.000 |
knitr::kable(
data.frame(
Metrik = "Balanced accuracy",
Nilai = round(iris_balanced_accuracy, 3)
),
caption = "Balanced accuracy model multinomial"
)| Metrik | Nilai |
|---|---|
| Balanced accuracy | 0.933 |
Interpretasi output
Recall menunjukkan proporsi observasi aktual suatu kelas yang berhasil dikenali oleh model. Precision menunjukkan proporsi prediksi suatu kelas yang benar-benar berasal dari kelas tersebut. Balanced accuracy adalah rata-rata recall seluruh kelas.
6.10 Pemeriksaan Asumsi Model Multinomial
Pemeriksaan asumsi regresi logistik multinomial tidak hanya dilakukan setelah model jadi, tetapi juga dimulai dari bentuk variabel respon dan struktur data. Pada dataset Iris, pemeriksaan asumsi lebih bersifat teknis karena data ini memang dibuat untuk klasifikasi spesies berdasarkan indikator morfologi.
6.10.1 Respons Bersifat Nominal
Regresi logistik multinomial digunakan ketika kategori respons tidak memiliki urutan alami. Pada data Iris, kategori Iris-setosa, Iris-versicolor, dan Iris-virginica merupakan jenis spesies, bukan tingkatan dari rendah ke tinggi.
iris_response_check <- data.frame(
Aspek = c("Nama variabel respons", "Jumlah kategori", "Sifat kategori"),
Hasil = c(
"species",
nlevels(iris_data$species),
"Nominal, karena spesies tidak memiliki urutan alami"
)
)
knitr::kable(
iris_response_check,
caption = "Pemeriksaan skala respons multinomial"
)| Aspek | Hasil |
|---|---|
| Nama variabel respons | species |
| Jumlah kategori | 3 |
| Sifat kategori | Nominal, karena spesies tidak memiliki urutan alami |
Interpretasi output
Karena respons memiliki lebih dari dua kategori dan tidak berurutan, penggunaan regresi logistik multinomial sesuai secara konseptual. Jika kategori memiliki urutan alami, model ordinal lebih tepat.
6.10.2 Observasi Independen
Model mengasumsikan bahwa satu observasi tidak bergantung pada observasi lainnya. Dalam dataset Iris, setiap baris merepresentasikan satu sampel bunga. Dengan demikian, observasi diperlakukan independen.
Jika data berasal dari pengukuran berulang pada objek yang sama, lokasi yang sama, atau kelompok yang sama, maka asumsi independensi dapat terganggu dan model dengan efek acak atau struktur klaster perlu dipertimbangkan.
6.10.3 Kecukupan Frekuensi Setiap Kategori
Kategori respons sebaiknya memiliki jumlah observasi yang cukup. Kategori yang terlalu sedikit dapat menyebabkan estimasi tidak stabil.
iris_frequency_check <- iris_data %>%
count(species, name = "Jumlah") %>%
mutate(
Proporsi = scales::percent(Jumlah / sum(Jumlah), accuracy = 0.1),
`Status frekuensi` = ifelse(Jumlah >= 10, "Memadai", "Perlu perhatian")
) %>%
rename(`Spesies iris` = species)
knitr::kable(
iris_frequency_check,
caption = "Pemeriksaan frekuensi kategori respons"
)| Spesies iris | Jumlah | Proporsi | Status frekuensi |
|---|---|---|---|
| Iris-setosa | 50 | 33.3% | Memadai |
| Iris-versicolor | 50 | 33.3% | Memadai |
| Iris-virginica | 50 | 33.3% | Memadai |
Interpretasi output
Setiap kategori spesies memiliki frekuensi yang memadai. Hal ini mendukung kestabilan estimasi model multinomial.
6.10.4 Tidak Ada Multikolinearitas Berat
Prediktor numerik sebaiknya tidak saling berkorelasi terlalu tinggi. Korelasi yang sangat tinggi dapat membuat koefisien tidak stabil dan interpretasi individual menjadi sulit.
iris_numeric_predictors <- iris_data %>%
dplyr::select(sepal_length, sepal_width, petal_length, petal_width)
iris_cor_matrix <- cor(iris_numeric_predictors)
knitr::kable(
round(iris_cor_matrix, 3),
caption = "Matriks korelasi antar prediktor numerik"
)| sepal_length | sepal_width | petal_length | petal_width | |
|---|---|---|---|---|
| sepal_length | 1.000 | -0.109 | 0.872 | 0.818 |
| sepal_width | -0.109 | 1.000 | -0.421 | -0.357 |
| petal_length | 0.872 | -0.421 | 1.000 | 0.963 |
| petal_width | 0.818 | -0.357 | 0.963 | 1.000 |
iris_high_cor <- as.data.frame(as.table(abs(iris_cor_matrix))) %>%
filter(Var1 != Var2, Freq > 0.8) %>%
mutate(Korelasi = round(Freq, 3)) %>%
dplyr::select(`Prediktor 1` = Var1, `Prediktor 2` = Var2, Korelasi)
if (nrow(iris_high_cor) == 0) {
iris_high_cor <- data.frame(
`Prediktor 1` = "-",
`Prediktor 2` = "-",
Korelasi = "Tidak ada korelasi absolut di atas 0.8"
)
}
knitr::kable(
iris_high_cor,
caption = "Pasangan prediktor dengan korelasi tinggi"
)| Prediktor 1 | Prediktor 2 | Korelasi |
|---|---|---|
| petal_length | sepal_length | 0.872 |
| petal_width | sepal_length | 0.818 |
| sepal_length | petal_length | 0.872 |
| petal_width | petal_length | 0.963 |
| sepal_length | petal_width | 0.818 |
| petal_length | petal_width | 0.963 |
Interpretasi output
Jika terdapat korelasi absolut di atas 0.8, maka ada indikasi multikolinearitas. Pada data Iris, ukuran petal biasanya sangat berkorelasi karena keduanya mengukur bagian morfologi yang sama. Ini tidak membuat model otomatis salah, tetapi interpretasi koefisien per variabel harus dilakukan hati-hati.
6.10.5 Tidak Ada Perfect Separation
Perfect separation terjadi ketika kombinasi prediktor mampu memisahkan kategori respons secara sempurna. Jika terjadi, koefisien dapat menjadi sangat besar dan estimasi tidak stabil.
iris_train_prob_max <- apply(iris_prob_test, 1, max)
iris_separation_check <- data.frame(
Indikator = c(
"Probabilitas prediksi maksimum rata-rata",
"Jumlah observasi testing dengan probabilitas maksimum > 0.999"
),
Nilai = c(
round(mean(iris_train_prob_max), 4),
sum(iris_train_prob_max > 0.999)
)
)
knitr::kable(
iris_separation_check,
caption = "Indikasi separation melalui probabilitas prediksi ekstrem"
)| Indikator | Nilai |
|---|---|
| Probabilitas prediksi maksimum rata-rata | 0.9927 |
| Jumlah observasi testing dengan probabilitas maksimum > 0.999 | 28.0000 |
Interpretasi output
Probabilitas yang sangat dekat dengan 1 menunjukkan model sangat yakin pada kelas tertentu. Pada Iris, hal ini dapat terjadi karena spesies memang mudah dipisahkan oleh ukuran petal. Secara klasifikasi ini baik, tetapi untuk inferensi koefisien perlu diingat bahwa data terlalu bersih dapat menghasilkan estimasi yang sangat tajam.
6.10.6 Linearitas pada Skala Logit
Untuk prediktor kontinu, model mengasumsikan hubungan yang linear terhadap log-odds kategori target relatif terhadap kategori referensi. Pemeriksaan praktis dapat dilakukan dengan membandingkan model linear dengan model yang menambahkan suku kuadrat.
iris_multinom_quad <- nnet::multinom(
species ~ sepal_length + sepal_width + petal_length + petal_width +
I(sepal_length^2) + I(sepal_width^2) + I(petal_length^2) + I(petal_width^2),
data = iris_train,
trace = FALSE
)
iris_linearity_check <- data.frame(
Model = c("Model linear", "Model dengan suku kuadrat"),
AIC = c(
round(AIC(iris_multinom_fit), 3),
round(AIC(iris_multinom_quad), 3)
)
)
knitr::kable(
iris_linearity_check,
caption = "Perbandingan indikasi linearitas pada skala logit"
)| Model | AIC |
|---|---|
| Model linear | 23.262 |
| Model dengan suku kuadrat | 36.000 |
Interpretasi output
Jika AIC model dengan suku kuadrat jauh lebih kecil, terdapat indikasi hubungan non-linear pada skala logit. Jika perbedaannya kecil, model linear lebih sederhana dan masih dapat dipertahankan.
6.10.7 Independence of Irrelevant Alternatives (IIA)
Asumsi IIA menyatakan bahwa odds relatif antara dua kategori tidak berubah secara substansial hanya karena kategori lain ditambahkan atau dihilangkan. Dalam konteks Iris, asumsi ini dapat diperdebatkan karena spesies versicolor dan virginica secara morfologi lebih mirip dibandingkan dengan setosa.
iris_no_setosa <- iris_train %>%
filter(species != "Iris-setosa") %>%
mutate(species = droplevels(species))
iris_binary_like_fit <- nnet::multinom(
species ~ sepal_length + sepal_width + petal_length + petal_width,
data = iris_no_setosa,
trace = FALSE
)
iia_note_iris <- data.frame(
Aspek = c("Kategori yang dibandingkan ulang", "Tujuan pemeriksaan", "Catatan"),
Keterangan = c(
"Iris-versicolor dan Iris-virginica setelah Iris-setosa dikeluarkan",
"Melihat sensitivitas model ketika satu alternatif dihilangkan",
"Pemeriksaan ini bersifat diagnostik sederhana, bukan uji Hausman formal"
)
)
knitr::kable(
iia_note_iris,
caption = "Pemeriksaan konseptual asumsi IIA"
)| Aspek | Keterangan |
|---|---|
| Kategori yang dibandingkan ulang | Iris-versicolor dan Iris-virginica setelah Iris-setosa dikeluarkan |
| Tujuan pemeriksaan | Melihat sensitivitas model ketika satu alternatif dihilangkan |
| Catatan | Pemeriksaan ini bersifat diagnostik sederhana, bukan uji Hausman formal |
Interpretasi output
Jika kesimpulan antar kategori berubah drastis setelah satu kategori dikeluarkan, asumsi IIA perlu dicurigai. Pada praktik lanjutan, pemeriksaan IIA dapat dilakukan dengan uji Hausman-McFadden atau membandingkan model pilihan diskret yang lebih fleksibel.
6.10.8 Ringkasan Pemeriksaan Asumsi Multinomial
iris_assumption_summary <- data.frame(
Asumsi = c(
"Respons nominal",
"Observasi independen",
"Frekuensi kategori memadai",
"Tidak ada multikolinearitas berat",
"Tidak ada perfect separation ekstrem",
"Linearitas pada skala logit",
"IIA"
),
`Cara pemeriksaan` = c(
"Cek jenis kategori species",
"Cek desain data; satu baris satu bunga",
"Tabel frekuensi species",
"Matriks korelasi prediktor numerik",
"Probabilitas prediksi ekstrem",
"Bandingkan model linear dan kuadrat",
"Pemeriksaan sensitivitas kategori"
),
`Catatan pada data Iris` = c(
"Sesuai",
"Diasumsikan sesuai",
"Memadai dan seimbang",
"Perlu hati-hati karena ukuran petal berkorelasi",
"Mungkin muncul karena data sangat mudah dipisahkan",
"Dibaca melalui perbandingan AIC",
"Perlu catatan karena spesies tidak sepenuhnya substitusi bebas"
)
)
knitr::kable(
iris_assumption_summary,
caption = "Ringkasan asumsi regresi logistik multinomial"
)| Asumsi | Cara.pemeriksaan | Catatan.pada.data.Iris |
|---|---|---|
| Respons nominal | Cek jenis kategori species | Sesuai |
| Observasi independen | Cek desain data; satu baris satu bunga | Diasumsikan sesuai |
| Frekuensi kategori memadai | Tabel frekuensi species | Memadai dan seimbang |
| Tidak ada multikolinearitas berat | Matriks korelasi prediktor numerik | Perlu hati-hati karena ukuran petal berkorelasi |
| Tidak ada perfect separation ekstrem | Probabilitas prediksi ekstrem | Mungkin muncul karena data sangat mudah dipisahkan |
| Linearitas pada skala logit | Bandingkan model linear dan kuadrat | Dibaca melalui perbandingan AIC |
| IIA | Pemeriksaan sensitivitas kategori | Perlu catatan karena spesies tidak sepenuhnya substitusi bebas |
Jika asumsi tidak terpenuhi, beberapa alternatif yang dapat dipertimbangkan adalah menyederhanakan prediktor, menggabungkan kategori yang sangat kecil, menggunakan regularisasi, atau memakai model pilihan diskret yang lebih fleksibel seperti nested logit atau multinomial probit.
7 Regresi Logistik Ordinal: Car Evaluation
7.1 Ringkasan Masalah
Regresi logistik ordinal digunakan ketika variabel respon memiliki kategori yang berurutan. Pada bab ini digunakan dataset Car Evaluation dari UCI Machine Learning Repository. Variabel respon adalah tingkat evaluasi mobil:
\[ \text{unacc} < \text{acc} < \text{good} < \text{vgood} \]
Urutan tersebut menunjukkan tingkat penerimaan mobil dari paling rendah sampai paling tinggi.
7.2 Formulasi Model Ordinal
Model yang digunakan adalah proportional odds model atau cumulative logit model:
\[ \log\left(\frac{P(Y_i \le j)}{P(Y_i > j)}\right) = \alpha_j - \beta_1X_{1i} - \beta_2X_{2i} - \cdots - \beta_kX_{ki}. \]
Model ini memiliki beberapa titik potong \(\alpha_j\), tetapi koefisien \(\beta\) diasumsikan sama pada seluruh batas kategori. Asumsi ini disebut proportional odds assumption.
7.3 Sumber dan Persiapan Data
Dataset Car Evaluation memiliki 1.728 observasi, 6 prediktor kategorik, dan 4 tingkat kelas evaluasi. Dataset asli dapat diakses melalui UCI Machine Learning Repository:
https://archive.ics.uci.edu/dataset/19/car+evaluation
car_zip <- file.path(
"C:", "Users", "ASUS", "Documents", "CHEL", "Semester 4",
"Analisis Data Kategori", "Data Rpubs", "car+evaluation.zip"
)
if (!file.exists(car_zip)) {
stop("File car+evaluation.zip tidak ditemukan pada path lokal.")
}
car_raw <- read.csv(
unz(car_zip, "car.data"),
header = FALSE,
col.names = c("buying", "maint", "doors", "persons", "lug_boot", "safety", "class")
)
car_data <- car_raw %>%
mutate(
buying = factor(buying, levels = c("low", "med", "high", "vhigh")),
maint = factor(maint, levels = c("low", "med", "high", "vhigh")),
doors = factor(doors, levels = c("2", "3", "4", "5more")),
persons = factor(persons, levels = c("2", "4", "more")),
lug_boot = factor(lug_boot, levels = c("small", "med", "big")),
safety = factor(safety, levels = c("low", "med", "high")),
class_ord = factor(class, levels = c("unacc", "acc", "good", "vgood"), ordered = TRUE)
)
car_summary <- data.frame(
Keterangan = c("Jumlah observasi", "Jumlah prediktor", "Jumlah kategori respon"),
Nilai = c(nrow(car_data), 6, nlevels(car_data$class_ord))
)
knitr::kable(
car_summary,
caption = "Ringkasan dataset Car Evaluation"
)| Keterangan | Nilai |
|---|---|
| Jumlah observasi | 1728 |
| Jumlah prediktor | 6 |
| Jumlah kategori respon | 4 |
Interpretasi output
Dataset Car Evaluation cocok untuk regresi logistik ordinal karena variabel respon memiliki tingkatan yang jelas dari unacc sampai vgood.
7.4 Distribusi Respons Ordinal
car_class_summary <- car_data %>%
count(class_ord, name = "Jumlah") %>%
mutate(Proporsi = scales::percent(Jumlah / sum(Jumlah), accuracy = 0.1)) %>%
rename(`Evaluasi mobil` = class_ord)
knitr::kable(
car_class_summary,
caption = "Distribusi kelas evaluasi mobil"
)| Evaluasi mobil | Jumlah | Proporsi |
|---|---|---|
| unacc | 1210 | 70.0% |
| acc | 384 | 22.2% |
| good | 69 | 4.0% |
| vgood | 65 | 3.8% |
ggplot(car_data, aes(x = class_ord, fill = class_ord)) +
geom_bar(width = 0.65, color = "white", linewidth = 0.8) +
geom_text(
stat = "count",
aes(label = after_stat(count)),
vjust = -0.35,
fontface = "bold"
) +
scale_fill_manual(values = c("#d63384", "#f783ac", "#ffe066", "#2a9d8f")) +
labs(
title = "Distribusi Tingkat Evaluasi Mobil",
x = "Kelas evaluasi",
y = "Jumlah mobil"
) +
theme_minimal(base_size = 12) +
theme(legend.position = "none")
Interpretasi output
Distribusi kelas menunjukkan apakah kategori evaluasi mobil seimbang atau tidak. Jika kategori unacc mendominasi, maka evaluasi model perlu memperhatikan bahwa akurasi dapat dipengaruhi oleh kelas mayoritas.
7.5 Estimasi Model Ordinal
Dataset Car Evaluation dibentuk dari aturan evaluasi yang sangat
terstruktur. Variabel persons dan
safety dapat menimbulkan perfect separation pada model ordinal penuh,
misalnya kategori tertentu hampir selalu menghasilkan evaluasi unacc. Agar estimasi polr()
stabil saat dokumen di-knit, model utama di bawah menggunakan prediktor
buying, maint, doors, dan lug_boot.
set.seed(123)
car_train_id <- unlist(
lapply(split(seq_len(nrow(car_data)), car_data$class_ord), function(idx) {
sample(idx, size = floor(0.8 * length(idx)))
}),
use.names = FALSE
)
car_train <- car_data[car_train_id, ]
car_test <- car_data[-car_train_id, ]
car_polr_fit <- MASS::polr(
class_ord ~ buying + maint + doors + lug_boot,
data = car_train,
Hess = TRUE,
method = "logistic"
)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredcar_model_summary <- data.frame(
Keterangan = c("Jumlah observasi training", "AIC", "Deviance residual"),
Nilai = c(
nrow(car_train),
round(AIC(car_polr_fit), 3),
round(deviance(car_polr_fit), 3)
)
)
knitr::kable(
car_model_summary,
caption = "Ringkasan model regresi logistik ordinal"
)| Keterangan | Nilai |
|---|---|
| Jumlah observasi training | 1382.000 |
| AIC | 2125.041 |
| Deviance residual | 2097.041 |
Interpretasi output
Model ordinal mengestimasi pengaruh prediktor terhadap kecenderungan mobil berada pada kategori evaluasi yang lebih tinggi atau lebih rendah. Nilai AIC dapat digunakan untuk membandingkan beberapa spesifikasi model ordinal.
7.6 Likelihood, Cutpoint, dan Konvensi Tanda Model Ordinal
Model ordinal proportional odds juga diestimasi dengan maximum
likelihood. Untuk respons ordinal dengan \(J\) kategori, model membentuk \(J-1\) batas kumulatif. Pada
MASS::polr(), bentuk model yang digunakan adalah:
\[ \text{logit}\{P(Y \le j)\} = \zeta_j - \eta_i, \]
dengan:
\[ \eta_i = \beta_1X_{1i}+\beta_2X_{2i}+\cdots+\beta_kX_{ki}. \]
Parameter \(\zeta_j\) disebut cutpoint atau threshold. Cutpoint memisahkan kategori ordinal pada skala laten.
car_cutpoint_table <- coef(summary(car_polr_fit)) %>%
as.data.frame()
car_cutpoint_table$term <- rownames(car_cutpoint_table)
rownames(car_cutpoint_table) <- NULL
car_cutpoint_table <- car_cutpoint_table %>%
filter(grepl("\\|", term)) %>%
transmute(
Cutpoint = term,
Nilai = round(Value, 3),
`Std. Error` = round(`Std. Error`, 3),
`t-value` = round(`t value`, 3)
)
knitr::kable(
car_cutpoint_table,
caption = "Cutpoint pada model ordinal"
)| Cutpoint | Nilai | Std. Error | t-value |
|---|---|---|---|
| unacc|acc | 0.533 | 0.211 | 2.529 |
| acc|good | 2.329 | 0.224 | 10.376 |
| good|vgood | 3.172 | 0.245 | 12.971 |
Interpretasi output
Cutpoint bukan prediktor, melainkan batas antar kategori pada skala
laten. Pada output polr(), koefisien positif untuk
prediktor berarti nilai prediktor yang lebih tinggi cenderung menggeser
peluang ke kategori respons yang lebih tinggi, dengan asumsi variabel
lain konstan.
7.7 Goodness-of-Fit Model Ordinal
Goodness-of-fit model ordinal dapat dilihat dari AIC, log-likelihood, pseudo-\(R^2\), dan likelihood ratio test terhadap model null.
car_polr_null <- MASS::polr(
class_ord ~ 1,
data = car_train,
Hess = TRUE,
method = "logistic"
)
ll_car_full <- as.numeric(logLik(car_polr_fit))
ll_car_null <- as.numeric(logLik(car_polr_null))
lr_car <- 2 * (ll_car_full - ll_car_null)
df_car <- attr(logLik(car_polr_fit), "df") -
attr(logLik(car_polr_null), "df")
car_gof_table <- data.frame(
Ukuran = c(
"Log-likelihood model null",
"Log-likelihood model penuh",
"Likelihood ratio statistic",
"Derajat bebas LR",
"p-value LR",
"McFadden pseudo R2",
"AIC model penuh"
),
Nilai = c(
round(ll_car_null, 3),
round(ll_car_full, 3),
round(lr_car, 3),
df_car,
signif(pchisq(lr_car, df = df_car, lower.tail = FALSE), 3),
round(1 - ll_car_full / ll_car_null, 3),
round(AIC(car_polr_fit), 3)
)
)
knitr::kable(
car_gof_table,
caption = "Goodness-of-fit model ordinal"
)| Ukuran | Nilai |
|---|---|
| Log-likelihood model null | -1154.404 |
| Log-likelihood model penuh | -1048.520 |
| Likelihood ratio statistic | 211.767 |
| Derajat bebas LR | 11.000 |
| p-value LR | 0.000 |
| McFadden pseudo R2 | 0.092 |
| AIC model penuh | 2125.041 |
Interpretasi output
Jika p-value likelihood ratio test kecil, model dengan prediktor lebih baik daripada model null. AIC digunakan untuk membandingkan beberapa model, sedangkan pseudo-\(R^2\) memberi ukuran ringkas peningkatan kecocokan model.
7.8 Ringkasan Koefisien Ordinal
car_coef <- coef(summary(car_polr_fit))
car_coef_table <- as.data.frame(car_coef)
car_coef_table$term <- rownames(car_coef_table)
rownames(car_coef_table) <- NULL
car_coef_table <- car_coef_table %>%
filter(!grepl("\\|", term)) %>%
mutate(
`p-value` = 2 * (1 - pnorm(abs(`t value`))),
`Odds ratio` = exp(Value)
) %>%
transmute(
`Variabel/level` = term,
Koefisien = round(Value, 3),
`Std. Error` = round(`Std. Error`, 3),
`t-value` = round(`t value`, 3),
`Odds ratio` = round(`Odds ratio`, 3),
`p-value` = signif(`p-value`, 3)
)
knitr::kable(
car_coef_table,
caption = "Koefisien model ordinal dan odds ratio"
)| Variabel/level | Koefisien | Std. Error | t-value | Odds ratio | p-value |
|---|---|---|---|---|---|
| buyingmed | -0.117 | 0.161 | -0.727 | 0.889 | 4.67e-01 |
| buyinghigh | -1.049 | 0.172 | -6.103 | 0.350 | 0.00e+00 |
| buyingvhigh | -1.503 | 0.186 | -8.089 | 0.222 | 0.00e+00 |
| maintmed | -0.158 | 0.159 | -0.991 | 0.854 | 3.22e-01 |
| mainthigh | -0.853 | 0.172 | -4.959 | 0.426 | 7.00e-07 |
| maintvhigh | -1.550 | 0.189 | -8.211 | 0.212 | 0.00e+00 |
| doors3 | 0.300 | 0.178 | 1.689 | 1.350 | 9.12e-02 |
| doors4 | 0.498 | 0.177 | 2.814 | 1.645 | 4.89e-03 |
| doors5more | 0.424 | 0.178 | 2.382 | 1.528 | 1.72e-02 |
| lug_bootmed | 0.761 | 0.159 | 4.803 | 2.141 | 1.60e-06 |
| lug_bootbig | 0.930 | 0.158 | 5.892 | 2.534 | 0.00e+00 |
Interpretasi output
Pada model MASS::polr, koefisien positif menunjukkan
kecenderungan lebih besar menuju kategori evaluasi yang lebih tinggi.
Odds ratio di atas 1 menunjukkan peningkatan odds untuk berada pada
kategori evaluasi yang lebih baik, sedangkan odds ratio di bawah 1
menunjukkan kecenderungan ke kategori yang lebih rendah.
7.9 Prediksi Probabilitas dan Confusion Matrix
car_prob_test <- predict(car_polr_fit, newdata = car_test, type = "probs")
car_pred_test <- predict(car_polr_fit, newdata = car_test, type = "class")
car_confusion <- table(
Aktual = car_test$class_ord,
Prediksi = car_pred_test
)
car_accuracy <- sum(diag(car_confusion)) / sum(car_confusion)
knitr::kable(
addmargins(car_confusion),
caption = "Confusion matrix model ordinal pada data testing"
)| unacc | acc | good | vgood | Sum | |
|---|---|---|---|---|---|
| unacc | 227 | 15 | 0 | 0 | 242 |
| acc | 75 | 2 | 0 | 0 | 77 |
| good | 9 | 5 | 0 | 0 | 14 |
| vgood | 5 | 8 | 0 | 0 | 13 |
| Sum | 316 | 30 | 0 | 0 | 346 |
car_accuracy_table <- data.frame(
Metrik = "Akurasi testing",
Nilai = round(car_accuracy, 3)
)
knitr::kable(
car_accuracy_table,
caption = "Akurasi model ordinal"
)| Metrik | Nilai |
|---|---|
| Akurasi testing | 0.662 |
Interpretasi output
Confusion matrix memperlihatkan apakah model mampu memprediksi tingkat evaluasi mobil secara tepat. Kesalahan prediksi yang hanya bergeser satu tingkat, misalnya dari good ke acc, biasanya lebih ringan dibandingkan kesalahan ekstrem dari vgood ke unacc.
car_prob_plot <- as.data.frame(car_prob_test)
car_prob_plot$Nomor <- seq_len(nrow(car_prob_plot))
car_prob_plot$Aktual <- car_test$class_ord
car_prob_long <- car_prob_plot %>%
tidyr::pivot_longer(
cols = all_of(levels(car_data$class_ord)),
names_to = "Kelas",
values_to = "Probabilitas"
)ggplot(car_prob_long, aes(x = Nomor, y = Probabilitas, color = Kelas)) +
geom_line(linewidth = 0.85) +
scale_color_manual(values = c("#d63384", "#f783ac", "#ffe066", "#2a9d8f")) +
labs(
title = "Probabilitas Prediksi Ordinal pada Data Testing",
x = "Observasi testing",
y = "Probabilitas prediksi",
color = "Kelas"
) +
theme_minimal(base_size = 12)
Interpretasi output
Grafik memperlihatkan probabilitas prediksi untuk setiap tingkat evaluasi. Model ordinal tidak hanya memberikan satu kelas prediksi, tetapi juga peluang untuk masing-masing kategori berurutan.
7.10 Metrik Evaluasi Ordinal
Pada respons ordinal, kesalahan prediksi dapat dibaca dari jarak antar kategori. Salah prediksi satu tingkat biasanya lebih ringan daripada salah prediksi beberapa tingkat.
car_actual_score <- as.numeric(car_test$class_ord)
car_pred_score <- as.numeric(car_pred_test)
car_ordinal_metrics <- data.frame(
Metrik = c(
"Akurasi",
"Mean absolute ordinal error",
"Proporsi salah maksimal 1 tingkat"
),
Nilai = c(
round(car_accuracy, 3),
round(mean(abs(car_actual_score - car_pred_score)), 3),
round(mean(abs(car_actual_score - car_pred_score) <= 1), 3)
)
)
knitr::kable(
car_ordinal_metrics,
caption = "Metrik evaluasi tambahan untuk respons ordinal"
)| Metrik | Nilai |
|---|---|
| Akurasi | 0.662 |
| Mean absolute ordinal error | 0.416 |
| Proporsi salah maksimal 1 tingkat | 0.936 |
Interpretasi output
Mean absolute ordinal error menunjukkan rata-rata jarak kesalahan prediksi dalam satuan tingkat kategori. Nilai yang lebih kecil menunjukkan prediksi lebih dekat dengan kategori aktual.
7.11 Pemeriksaan Asumsi Model Ordinal
Regresi logistik ordinal memiliki asumsi yang lebih spesifik dibandingkan multinomial karena model memanfaatkan urutan kategori respons. Pemeriksaan utama adalah apakah kategori respons benar-benar ordinal dan apakah asumsi proportional odds masuk akal.
7.11.1 Respons Memiliki Urutan Alami
Pada data Car Evaluation, kelas respons disusun sebagai:
\[ \text{unacc} < \text{acc} < \text{good} < \text{vgood}. \]
car_response_check <- data.frame(
Aspek = c("Nama variabel respons", "Urutan kategori", "Sifat kategori"),
Hasil = c(
"class_ord",
paste(levels(car_data$class_ord), collapse = " < "),
"Ordinal, karena kategori menunjukkan tingkat evaluasi mobil"
)
)
knitr::kable(
car_response_check,
caption = "Pemeriksaan skala respons ordinal"
)| Aspek | Hasil |
|---|---|
| Nama variabel respons | class_ord |
| Urutan kategori | unacc < acc < good < vgood |
| Sifat kategori | Ordinal, karena kategori menunjukkan tingkat evaluasi mobil |
Interpretasi output
Karena kategori respons memiliki urutan yang jelas, regresi logistik ordinal lebih efisien daripada multinomial jika asumsi proportional odds masih dapat diterima.
7.11.2 Observasi Independen
Setiap baris pada dataset Car Evaluation adalah satu kombinasi karakteristik mobil. Model memperlakukan setiap kombinasi sebagai observasi independen. Jika data berasal dari penilaian berulang oleh evaluator yang sama, maka struktur dependensi perlu dimodelkan secara khusus.
7.11.3 Frekuensi Kategori Respons
car_frequency_check <- car_data %>%
count(class_ord, name = "Jumlah") %>%
mutate(
Proporsi = scales::percent(Jumlah / sum(Jumlah), accuracy = 0.1),
`Status frekuensi` = ifelse(Jumlah >= 10, "Memadai", "Perlu perhatian")
) %>%
rename(`Kelas evaluasi` = class_ord)
knitr::kable(
car_frequency_check,
caption = "Pemeriksaan frekuensi kategori ordinal"
)| Kelas evaluasi | Jumlah | Proporsi | Status frekuensi |
|---|---|---|---|
| unacc | 1210 | 70.0% | Memadai |
| acc | 384 | 22.2% | Memadai |
| good | 69 | 4.0% | Memadai |
| vgood | 65 | 3.8% | Memadai |
Interpretasi output
Kategori yang sangat kecil dapat membuat estimasi cutpoint dan koefisien tidak stabil. Pada data Car Evaluation, kelas tidak sepenuhnya seimbang sehingga interpretasi akurasi perlu dilakukan dengan hati-hati.
7.11.4 Tidak Ada Separation yang Mengganggu Estimasi
Dataset Car Evaluation dibentuk dari aturan evaluasi. Karena itu,
kombinasi tertentu dapat secara deterministik menghasilkan kelas
tertentu. Inilah alasan model utama tidak memasukkan persons dan safety, karena keduanya dapat menyebabkan
separation pada polr().
car_sparse_persons <- table(car_data$persons, car_data$class_ord)
car_sparse_safety <- table(car_data$safety, car_data$class_ord)
knitr::kable(
car_sparse_persons,
caption = "Tabulasi persons terhadap kelas evaluasi"
)| unacc | acc | good | vgood | |
|---|---|---|---|---|
| 2 | 576 | 0 | 0 | 0 |
| 4 | 312 | 198 | 36 | 30 |
| more | 322 | 186 | 33 | 35 |
| unacc | acc | good | vgood | |
|---|---|---|---|---|
| low | 576 | 0 | 0 | 0 |
| med | 357 | 180 | 39 | 0 |
| high | 277 | 204 | 30 | 65 |
Interpretasi output
Jika ada baris tabel yang hanya muncul pada satu kategori respons, maka separation dapat terjadi. Kondisi ini membuat estimasi maksimum likelihood sulit atau gagal konvergen.
7.11.5 Proportional Odds Assumption
Asumsi proportional odds menyatakan bahwa pengaruh prediktor sama pada setiap batas kumulatif kategori:
\[ \text{unacc} \;|\; \text{acc/good/vgood}, \]
\[ \text{unacc/acc} \;|\; \text{good/vgood}, \]
\[ \text{unacc/acc/good} \;|\; \text{vgood}. \]
Pemeriksaan praktis dapat dilakukan dengan membuat beberapa model logistik biner untuk setiap batas kumulatif, lalu membandingkan koefisiennya. Jika koefisien sangat berbeda antar batas, maka asumsi proportional odds perlu dicurigai.
car_cum_data <- car_train %>%
mutate(
class_score = as.numeric(class_ord),
y_le_unacc = as.integer(class_score <= 1),
y_le_acc = as.integer(class_score <= 2),
y_le_good = as.integer(class_score <= 3)
)
car_bin_formula <- ~ buying + maint + doors + lug_boot
car_cum_fits <- list(
`unacc vs lainnya` = glm(
y_le_unacc ~ buying + maint + doors + lug_boot,
data = car_cum_data,
family = binomial()
),
`unacc/acc vs good/vgood` = glm(
y_le_acc ~ buying + maint + doors + lug_boot,
data = car_cum_data,
family = binomial()
),
`unacc/acc/good vs vgood` = glm(
y_le_good ~ buying + maint + doors + lug_boot,
data = car_cum_data,
family = binomial()
)
)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredcar_cum_coef <- lapply(names(car_cum_fits), function(nm) {
out <- as.data.frame(coef(summary(car_cum_fits[[nm]])))
out$term <- rownames(out)
out$batas <- nm
rownames(out) <- NULL
out
}) %>%
bind_rows() %>%
filter(term != "(Intercept)") %>%
dplyr::select(batas, term, Estimate) %>%
tidyr::pivot_wider(names_from = batas, values_from = Estimate)
knitr::kable(
car_cum_coef %>% mutate(across(where(is.numeric), ~ round(.x, 3))),
caption = "Perbandingan koefisien pada beberapa batas kumulatif"
)| term | unacc vs lainnya | unacc/acc vs good/vgood | unacc/acc/good vs vgood |
|---|---|---|---|
| buyingmed | 0.002 | 0.516 | 0.092 |
| buyinghigh | 0.689 | 20.120 | 19.952 |
| buyingvhigh | 1.164 | 20.126 | 19.956 |
| maintmed | 0.034 | 0.692 | 0.114 |
| mainthigh | 0.587 | 2.237 | 0.897 |
| maintvhigh | 1.225 | 20.073 | 19.693 |
| doors3 | -0.273 | -0.414 | -0.265 |
| doors4 | -0.450 | -0.619 | -0.438 |
| doors5more | -0.342 | -0.627 | -0.726 |
| lug_bootmed | -0.654 | -1.197 | -19.261 |
| lug_bootbig | -0.766 | -1.541 | -19.716 |
Interpretasi output
Jika arah dan besar koefisien relatif mirip pada semua batas kumulatif, asumsi proportional odds lebih masuk akal. Jika terdapat perbedaan besar atau perubahan tanda, model ordinal proportional odds mungkin terlalu membatasi.
7.11.6 Alternatif Jika Asumsi Proportional Odds Tidak Terpenuhi
Jika asumsi proportional odds tidak terpenuhi, beberapa alternatif yang dapat digunakan adalah:
- Partial proportional odds model, yaitu sebagian prediktor boleh memiliki efek berbeda pada tiap batas.
- Generalized ordered logit, yaitu model ordinal yang lebih fleksibel.
- Regresi logistik multinomial, jika urutan kategori tidak ingin dipaksakan.
- Menggabungkan kategori respons jika ada kategori yang terlalu jarang atau sulit dibedakan.
Sebagai pembanding sederhana, model multinomial dapat diestimasi pada respons yang sama.
car_train_multi <- car_train %>%
mutate(class_nom = factor(class_ord, ordered = FALSE))
car_multinom_compare <- nnet::multinom(
class_nom ~ buying + maint + doors + lug_boot,
data = car_train_multi,
trace = FALSE
)
car_model_compare <- data.frame(
Model = c("Ordinal proportional odds", "Multinomial pembanding"),
AIC = c(
round(AIC(car_polr_fit), 3),
round(AIC(car_multinom_compare), 3)
)
)
knitr::kable(
car_model_compare,
caption = "Perbandingan model ordinal dan multinomial pembanding"
)| Model | AIC |
|---|---|
| Ordinal proportional odds | 2125.041 |
| Multinomial pembanding | 1978.866 |
Interpretasi output
Model multinomial memiliki parameter lebih banyak karena tidak memaksakan slope yang sama pada batas kumulatif. Jika AIC multinomial jauh lebih kecil, asumsi ordinal yang lebih ketat perlu dievaluasi kembali.
7.11.7 Ringkasan Pemeriksaan Asumsi Ordinal
car_assumption_summary <- data.frame(
Asumsi = c(
"Respons ordinal",
"Observasi independen",
"Frekuensi kategori memadai",
"Tidak ada separation berat",
"Proportional odds",
"Interpretasi tanda sesuai konvensi model"
),
`Cara pemeriksaan` = c(
"Cek urutan class_ord",
"Cek struktur data",
"Tabel frekuensi kelas",
"Tabulasi prediktor terhadap respons",
"Bandingkan koefisien logit kumulatif",
"Perhatikan konvensi MASS::polr"
),
`Catatan pada data Car Evaluation` = c(
"Sesuai",
"Diasumsikan sesuai",
"Tidak sepenuhnya seimbang",
"persons dan safety berpotensi separation",
"Perlu diperiksa dari tabel koefisien kumulatif",
"Koefisien polr positif mengarah ke kategori lebih tinggi"
)
)
knitr::kable(
car_assumption_summary,
caption = "Ringkasan asumsi regresi logistik ordinal"
)| Asumsi | Cara.pemeriksaan | Catatan.pada.data.Car.Evaluation |
|---|---|---|
| Respons ordinal | Cek urutan class_ord | Sesuai |
| Observasi independen | Cek struktur data | Diasumsikan sesuai |
| Frekuensi kategori memadai | Tabel frekuensi kelas | Tidak sepenuhnya seimbang |
| Tidak ada separation berat | Tabulasi prediktor terhadap respons | persons dan safety berpotensi separation |
| Proportional odds | Bandingkan koefisien logit kumulatif | Perlu diperiksa dari tabel koefisien kumulatif |
| Interpretasi tanda sesuai konvensi model | Perhatikan konvensi MASS::polr | Koefisien polr positif mengarah ke kategori lebih tinggi |
Pemeriksaan asumsi ordinal penting karena model ini berada di antara dua pilihan: lebih hemat parameter daripada multinomial, tetapi lebih ketat karena mengasumsikan efek prediktor yang sama pada seluruh batas kategori.
8 Regresi Poisson: Bike Sharing
8.1 Ringkasan Masalah
Regresi Poisson digunakan untuk memodelkan variabel respon berupa data hitung. Pada bab ini digunakan dataset Bike Sharing dari UCI Machine Learning Repository. Variabel respon adalah cnt, yaitu jumlah total peminjaman sepeda.
Karena cnt merupakan jumlah kejadian dalam satu hari, regresi Poisson menjadi model awal yang sesuai untuk mempelajari hubungan antara kondisi kalender, musim, cuaca, dan jumlah peminjaman sepeda.
8.2 Formulasi Model Poisson
Jika \(Y_i\) menyatakan jumlah kejadian pada observasi ke-\(i\), maka:
\[ Y_i \sim \text{Poisson}(\mu_i). \]
Model Poisson menggunakan link log:
\[ \log(\mu_i) = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \cdots + \beta_kX_{ki}. \]
Dengan transformasi balik:
\[ \hat{\mu}_i = \exp(\hat{\eta}_i). \]
Interpretasi koefisien dilakukan melalui incidence rate ratio:
\[ IRR_j = \exp(\beta_j). \]
Jika \(IRR_j > 1\), maka prediktor meningkatkan rata-rata jumlah kejadian. Jika \(IRR_j < 1\), maka prediktor menurunkan rata-rata jumlah kejadian.
8.3 Likelihood, Deviance, dan Offset pada Regresi Poisson
Regresi Poisson diestimasi dengan maximum likelihood. Untuk observasi ke-\(i\), fungsi peluang Poisson adalah:
\[ P(Y_i=y_i) = \frac{e^{-\mu_i}\mu_i^{y_i}}{y_i!}. \]
Dengan link log:
\[ \log(\mu_i)=\eta_i=X_i^\top\beta. \]
Jika setiap observasi memiliki exposure berbeda \(t_i\), model menjadi:
\[ \log(\mu_i)=\log(t_i)+X_i^\top\beta. \]
Komponen \(\log(t_i)\) disebut offset. Pada data Bike Sharing harian, semua observasi memiliki exposure satu hari, sehingga offset tidak diperlukan.
8.4 Sumber dan Persiapan Data
Dataset Bike Sharing memuat jumlah peminjaman sepeda dengan informasi kalender dan cuaca. Dataset asli dapat diakses melalui UCI Machine Learning Repository:
https://archive.ics.uci.edu/dataset/275/bike+sharing+dataset
bike_zip <- file.path(
"C:", "Users", "ASUS", "Documents", "CHEL", "Semester 4",
"Analisis Data Kategori", "Data Rpubs", "bike+sharing+dataset.zip"
)
if (!file.exists(bike_zip)) {
stop("File bike+sharing+dataset.zip tidak ditemukan pada path lokal.")
}
bike_day <- read.csv(
unz(bike_zip, "day.csv"),
stringsAsFactors = FALSE
)
bike_data <- bike_day %>%
mutate(
season = factor(
season,
levels = 1:4,
labels = c("winter", "spring", "summer", "fall")
),
yr = factor(yr, levels = c(0, 1), labels = c("2011", "2012")),
mnth = factor(mnth),
holiday = factor(holiday, levels = c(0, 1), labels = c("Tidak", "Ya")),
weekday = factor(weekday),
workingday = factor(workingday, levels = c(0, 1), labels = c("Tidak", "Ya")),
weathersit = factor(
weathersit,
levels = 1:3,
labels = c("Cerah/berawan", "Berkabut/mendung", "Hujan/salju ringan")
)
)
bike_summary <- data.frame(
Keterangan = c("Jumlah observasi harian", "Rata-rata peminjaman", "Median peminjaman"),
Nilai = c(
nrow(bike_data),
round(mean(bike_data$cnt), 2),
median(bike_data$cnt)
)
)
knitr::kable(
bike_summary,
caption = "Ringkasan dataset Bike Sharing harian"
)| Keterangan | Nilai |
|---|---|
| Jumlah observasi harian | 731.00 |
| Rata-rata peminjaman | 4504.35 |
| Median peminjaman | 4548.00 |
Interpretasi output
Data harian Bike Sharing memiliki variabel respon berupa jumlah peminjaman sepeda per hari. Karena nilai respon berupa hitungan non-negatif, regresi Poisson dapat digunakan sebagai model awal.
8.5 Distribusi Respons Hitungan
ggplot(bike_data, aes(x = cnt)) +
geom_histogram(bins = 30, fill = "#d63384", color = "white", alpha = 0.85) +
labs(
title = "Distribusi Jumlah Peminjaman Sepeda Harian",
x = "Jumlah peminjaman sepeda",
y = "Frekuensi"
) +
theme_minimal(base_size = 12)
Interpretasi output
Histogram menunjukkan sebaran jumlah peminjaman sepeda harian. Jika distribusi tidak simetris dan nilai respon berupa hitungan, model Poisson dengan link log menjadi pilihan awal yang lebih sesuai dibandingkan regresi linear biasa.
bike_season_summary <- bike_data %>%
group_by(season) %>%
summarise(
Jumlah_hari = n(),
`Rata-rata cnt` = mean(cnt),
`Median cnt` = median(cnt),
.groups = "drop"
) %>%
mutate(across(where(is.numeric), round, 2))
knitr::kable(
bike_season_summary,
caption = "Ringkasan jumlah peminjaman menurut musim"
)| season | Jumlah_hari | Rata-rata cnt | Median cnt |
|---|---|---|---|
| winter | 181 | 2604.13 | 2209.0 |
| spring | 184 | 4992.33 | 4941.5 |
| summer | 188 | 5644.30 | 5353.5 |
| fall | 178 | 4728.16 | 4634.5 |
Interpretasi output
Ringkasan menurut musim membantu melihat apakah jumlah peminjaman sepeda berbeda antar musim. Perbedaan ini masuk akal karena aktivitas bersepeda sangat dipengaruhi oleh kondisi cuaca dan waktu dalam setahun.
8.6 Estimasi Model Poisson
bike_poisson_fit <- glm(
cnt ~ season + yr + mnth + holiday + weekday + workingday +
weathersit + temp + hum + windspeed,
data = bike_data,
family = poisson(link = "log")
)
bike_model_summary <- data.frame(
Keterangan = c(
"Jumlah observasi",
"Null deviance",
"Residual deviance",
"Derajat bebas residual",
"AIC"
),
Nilai = c(
nobs(bike_poisson_fit),
round(bike_poisson_fit$null.deviance, 3),
round(bike_poisson_fit$deviance, 3),
bike_poisson_fit$df.residual,
round(AIC(bike_poisson_fit), 3)
)
)
knitr::kable(
bike_model_summary,
caption = "Ringkasan model regresi Poisson"
)| Keterangan | Nilai |
|---|---|
| Jumlah observasi | 731.0 |
| Null deviance | 668800.8 |
| Residual deviance | 116239.6 |
| Derajat bebas residual | 703.0 |
| AIC | 123694.2 |
Interpretasi output
Jika residual deviance lebih kecil daripada null deviance, maka prediktor membantu menjelaskan variasi jumlah peminjaman sepeda. Namun, pada regresi Poisson perlu diperiksa apakah terjadi overdispersion.
8.7 Goodness-of-Fit Model Poisson
Goodness-of-fit model Poisson dapat dilihat dari deviance, AIC, likelihood ratio test, dan pseudo-\(R^2\).
bike_poisson_null <- glm(
cnt ~ 1,
data = bike_data,
family = poisson(link = "log")
)
ll_bike_full <- as.numeric(logLik(bike_poisson_fit))
ll_bike_null <- as.numeric(logLik(bike_poisson_null))
lr_bike <- 2 * (ll_bike_full - ll_bike_null)
df_bike <- attr(logLik(bike_poisson_fit), "df") -
attr(logLik(bike_poisson_null), "df")
bike_gof_table <- data.frame(
Ukuran = c(
"Log-likelihood model null",
"Log-likelihood model penuh",
"Likelihood ratio statistic",
"Derajat bebas LR",
"p-value LR",
"McFadden pseudo R2",
"Residual deviance",
"p-value deviance GOF",
"AIC model penuh"
),
Nilai = c(
round(ll_bike_null, 3),
round(ll_bike_full, 3),
round(lr_bike, 3),
df_bike,
signif(pchisq(lr_bike, df = df_bike, lower.tail = FALSE), 3),
round(1 - ll_bike_full / ll_bike_null, 3),
round(bike_poisson_fit$deviance, 3),
signif(pchisq(bike_poisson_fit$deviance, df = bike_poisson_fit$df.residual, lower.tail = FALSE), 3),
round(AIC(bike_poisson_fit), 3)
)
)
knitr::kable(
bike_gof_table,
caption = "Goodness-of-fit model Poisson"
)| Ukuran | Nilai |
|---|---|
| Log-likelihood model null | -338099.722 |
| Log-likelihood model penuh | -61819.120 |
| Likelihood ratio statistic | 552561.203 |
| Derajat bebas LR | 27.000 |
| p-value LR | 0.000 |
| McFadden pseudo R2 | 0.817 |
| Residual deviance | 116239.580 |
| p-value deviance GOF | 0.000 |
| AIC model penuh | 123694.241 |
Interpretasi output
Likelihood ratio test menunjukkan apakah prediktor secara bersama-sama memperbaiki model dibandingkan model null. p-value deviance GOF yang sangat kecil dapat menjadi tanda bahwa model Poisson sederhana belum menangkap seluruh variasi data, terutama jika disertai overdispersion.
8.8 Incidence Rate Ratio
bike_coef <- as.data.frame(summary(bike_poisson_fit)$coefficients)
bike_coef$term <- rownames(bike_coef)
rownames(bike_coef) <- NULL
bike_irr_table <- bike_coef %>%
filter(term != "(Intercept)") %>%
mutate(
IRR = exp(Estimate),
ci_low = exp(Estimate - 1.96 * `Std. Error`),
ci_high = exp(Estimate + 1.96 * `Std. Error`)
) %>%
arrange(`Pr(>|z|)`) %>%
transmute(
`Variabel/level` = term,
IRR = round(IRR, 3),
`Interval kepercayaan 95 persen` = paste0(
round(ci_low, 3), " - ", round(ci_high, 3)
),
`p-value` = signif(`Pr(>|z|)`, 3)
)
knitr::kable(
head(bike_irr_table, 20),
caption = "Dua puluh prediktor paling menonjol berdasarkan p-value"
)| Variabel/level | IRR | Interval kepercayaan 95 persen | p-value |
|---|---|---|---|
| seasonspring | 1.314 | 1.304 - 1.323 | 0 |
| seasonsummer | 1.293 | 1.283 - 1.304 | 0 |
| seasonfall | 1.558 | 1.545 - 1.57 | 0 |
| yr2012 | 1.576 | 1.572 - 1.58 | 0 |
| mnth3 | 1.277 | 1.267 - 1.287 | 0 |
| mnth5 | 1.253 | 1.24 - 1.267 | 0 |
| mnth9 | 1.293 | 1.279 - 1.307 | 0 |
| holidayYa | 0.843 | 0.837 - 0.85 | 0 |
| weekday5 | 1.089 | 1.084 - 1.093 | 0 |
| weekday6 | 1.094 | 1.089 - 1.098 | 0 |
| weathersitBerkabut/mendung | 0.900 | 0.898 - 0.903 | 0 |
| weathersitHujan/salju ringan | 0.480 | 0.474 - 0.485 | 0 |
| temp | 3.162 | 3.11 - 3.214 | 0 |
| hum | 0.729 | 0.721 - 0.738 | 0 |
| windspeed | 0.527 | 0.519 - 0.536 | 0 |
| mnth4 | 1.217 | 1.205 - 1.23 | 0 |
| weekday4 | 1.079 | 1.074 - 1.083 | 0 |
| mnth10 | 1.207 | 1.194 - 1.22 | 0 |
| weekday2 | 1.072 | 1.068 - 1.077 | 0 |
| weekday3 | 1.069 | 1.064 - 1.073 | 0 |
Interpretasi output
IRR menunjukkan perubahan multiplikatif pada rata-rata jumlah peminjaman sepeda. Misalnya, IRR sebesar 1.10 berarti rata-rata jumlah peminjaman meningkat sekitar 10% dibandingkan kategori referensi atau untuk setiap kenaikan satu unit prediktor numerik.
8.9 Prediksi Rate Poisson
bike_data$pred_cnt_poisson <- predict(
bike_poisson_fit,
newdata = bike_data,
type = "response"
)
bike_prediction_preview <- bike_data %>%
dplyr::select(dteday, season, weathersit, temp, hum, windspeed, cnt, pred_cnt_poisson) %>%
head(10) %>%
mutate(pred_cnt_poisson = round(pred_cnt_poisson, 2))
knitr::kable(
bike_prediction_preview,
caption = "Contoh nilai aktual dan prediksi jumlah peminjaman"
)| dteday | season | weathersit | temp | hum | windspeed | cnt | pred_cnt_poisson |
|---|---|---|---|---|---|---|---|
| 2011-01-01 | winter | Berkabut/mendung | 0.344167 | 0.805833 | 0.1604460 | 985 | 1768.53 |
| 2011-01-02 | winter | Berkabut/mendung | 0.363478 | 0.696087 | 0.2485390 | 801 | 1617.84 |
| 2011-01-03 | winter | Cerah/berawan | 0.196364 | 0.437273 | 0.2483090 | 1349 | 1702.28 |
| 2011-01-04 | winter | Cerah/berawan | 0.200000 | 0.590435 | 0.1602960 | 1562 | 1746.10 |
| 2011-01-05 | winter | Cerah/berawan | 0.226957 | 0.436957 | 0.1869000 | 1600 | 1852.09 |
| 2011-01-06 | winter | Cerah/berawan | 0.204348 | 0.518261 | 0.0895652 | 1606 | 1889.98 |
| 2011-01-07 | winter | Berkabut/mendung | 0.196522 | 0.498696 | 0.1687260 | 1510 | 1627.77 |
| 2011-01-08 | winter | Berkabut/mendung | 0.165000 | 0.535833 | 0.2668040 | 959 | 1463.79 |
| 2011-01-09 | winter | Cerah/berawan | 0.138333 | 0.434167 | 0.3619500 | 822 | 1400.56 |
| 2011-01-10 | winter | Cerah/berawan | 0.150833 | 0.482917 | 0.2232670 | 1321 | 1617.98 |
ggplot(bike_data, aes(x = cnt, y = pred_cnt_poisson)) +
geom_point(alpha = 0.65, color = "#d63384") +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "#6c757d") +
labs(
title = "Perbandingan Jumlah Aktual dan Prediksi Poisson",
x = "Jumlah aktual",
y = "Jumlah prediksi"
) +
theme_minimal(base_size = 12)
Interpretasi output
Titik yang dekat dengan garis diagonal menunjukkan prediksi yang mendekati nilai aktual. Jika banyak titik jauh dari garis diagonal, berarti model masih belum menangkap seluruh variasi jumlah peminjaman sepeda.
8.10 Interval Prediksi Rata-Rata Poisson
Prediksi model Poisson dapat disertai interval kepercayaan pada skala rata-rata hitungan. Perhitungan dilakukan pada skala link, kemudian dikembalikan ke skala respons dengan fungsi eksponensial.
bike_pred_link <- predict(
bike_poisson_fit,
newdata = bike_data,
type = "link",
se.fit = TRUE
)
bike_pred_interval_preview <- bike_data %>%
mutate(
fit_link = bike_pred_link$fit,
se_link = bike_pred_link$se.fit,
pred_low = exp(fit_link - 1.96 * se_link),
pred_high = exp(fit_link + 1.96 * se_link)
) %>%
dplyr::select(dteday, cnt, pred_cnt_poisson, pred_low, pred_high) %>%
head(10) %>%
mutate(across(c(pred_cnt_poisson, pred_low, pred_high), ~ round(.x, 2)))
knitr::kable(
bike_pred_interval_preview,
caption = "Contoh interval kepercayaan prediksi rata-rata Poisson"
)| dteday | cnt | pred_cnt_poisson | pred_low | pred_high |
|---|---|---|---|---|
| 2011-01-01 | 985 | 1768.53 | 1756.85 | 1780.29 |
| 2011-01-02 | 801 | 1617.84 | 1606.94 | 1628.82 |
| 2011-01-03 | 1349 | 1702.28 | 1691.37 | 1713.26 |
| 2011-01-04 | 1562 | 1746.10 | 1735.11 | 1757.17 |
| 2011-01-05 | 1600 | 1852.09 | 1840.32 | 1863.94 |
| 2011-01-06 | 1606 | 1889.98 | 1877.67 | 1902.38 |
| 2011-01-07 | 1510 | 1627.77 | 1616.94 | 1638.67 |
| 2011-01-08 | 959 | 1463.79 | 1454.23 | 1473.41 |
| 2011-01-09 | 822 | 1400.56 | 1391.09 | 1410.09 |
| 2011-01-10 | 1321 | 1617.98 | 1607.56 | 1628.47 |
Interpretasi output
Interval ini adalah interval untuk rata-rata jumlah peminjaman yang diprediksi model, bukan interval prediksi observasi individual. Jika ingin interval prediksi observasi individual, variasi Poisson tambahan perlu diperhitungkan.
8.11 Pemeriksaan Asumsi Regresi Poisson
Regresi Poisson memiliki asumsi khusus untuk data hitung. Pada data Bike Sharing, pemeriksaan paling penting adalah memastikan respons berupa hitungan non-negatif, memeriksa equidispersion, melihat pola residual, dan menilai apakah ada kebutuhan offset/exposure.
8.11.1 Respons Berupa Hitungan Non-Negatif
Variabel cnt harus berupa bilangan cacah: 0, 1, 2, dan seterusnya.
bike_count_check <- data.frame(
Aspek = c(
"Nilai minimum cnt",
"Apakah semua cnt bilangan non-negatif?",
"Apakah semua cnt bilangan bulat?",
"Proporsi cnt = 0"
),
Hasil = c(
min(bike_data$cnt),
all(bike_data$cnt >= 0),
all(bike_data$cnt == floor(bike_data$cnt)),
scales::percent(mean(bike_data$cnt == 0), accuracy = 0.1)
)
)
knitr::kable(
bike_count_check,
caption = "Pemeriksaan karakteristik respons hitungan"
)| Aspek | Hasil |
|---|---|
| Nilai minimum cnt | 22 |
| Apakah semua cnt bilangan non-negatif? | TRUE |
| Apakah semua cnt bilangan bulat? | TRUE |
| Proporsi cnt = 0 | 0.0% |
Interpretasi output
Jika semua nilai cnt non-negatif dan bulat, maka syarat
dasar regresi Poisson terpenuhi. Jika terdapat nilai negatif atau
kontinu, model Poisson tidak sesuai.
8.11.2 Observasi Independen
Setiap baris pada day.csv adalah satu hari. Secara
praktik, data harian berpotensi memiliki dependensi waktu karena
peminjaman hari ini dapat berkaitan dengan hari sebelumnya. Model
Poisson sederhana mengabaikan struktur waktu ini, sehingga hasilnya
dibaca sebagai model awal.
bike_time_check <- data.frame(
Aspek = c("Unit observasi", "Potensi masalah", "Catatan model"),
Keterangan = c(
"Satu baris = satu hari",
"Data harian dapat memiliki autokorelasi",
"Model Poisson biasa dipakai sebagai model awal tanpa struktur time series"
)
)
knitr::kable(
bike_time_check,
caption = "Catatan independensi observasi pada data harian"
)| Aspek | Keterangan |
|---|---|
| Unit observasi | Satu baris = satu hari |
| Potensi masalah | Data harian dapat memiliki autokorelasi |
| Catatan model | Model Poisson biasa dipakai sebagai model awal tanpa struktur time series |
8.11.3 Mean dan Varians Respons
Pada distribusi Poisson murni, mean dan varians sama. Pemeriksaan
awal dapat dilakukan dengan membandingkan mean dan varians
cnt.
bike_mean_variance <- data.frame(
Statistik = c("Mean cnt", "Varians cnt", "Rasio varians/mean"),
Nilai = c(
round(mean(bike_data$cnt), 3),
round(var(bike_data$cnt), 3),
round(var(bike_data$cnt) / mean(bike_data$cnt), 3)
)
)
knitr::kable(
bike_mean_variance,
caption = "Pemeriksaan awal mean dan varians respons"
)| Statistik | Nilai |
|---|---|
| Mean cnt | 4504.349 |
| Varians cnt | 3752788.208 |
| Rasio varians/mean | 833.148 |
Interpretasi output
Jika varians jauh lebih besar daripada mean, terdapat indikasi awal overdispersion. Namun pemeriksaan yang lebih relevan dilakukan setelah model diestimasi menggunakan residual Pearson.
8.11.4 Linearitas pada Skala Log
Regresi Poisson mengasumsikan bahwa prediktor berhubungan linear dengan log rata-rata hitungan. Pemeriksaan praktis dapat dilakukan melalui plot residual Pearson terhadap nilai prediksi.
bike_diag <- bike_data %>%
mutate(
fitted_poisson = fitted(bike_poisson_fit),
pearson_resid = residuals(bike_poisson_fit, type = "pearson")
)
ggplot(bike_diag, aes(x = fitted_poisson, y = pearson_resid)) +
geom_point(alpha = 0.6, color = "#d63384") +
geom_smooth(method = "loess", se = FALSE, color = "#2a9d8f", linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "#6c757d") +
labs(
title = "Residual Pearson terhadap Nilai Prediksi Poisson",
x = "Nilai prediksi",
y = "Residual Pearson"
) +
theme_minimal(base_size = 12)## `geom_smooth()` using formula = 'y ~ x'
Interpretasi output
Jika pola residual membentuk lengkungan jelas, maka hubungan pada skala log mungkin belum cukup linear. Alternatifnya adalah menambahkan transformasi, suku kuadrat, spline, atau interaksi yang relevan.
8.11.5 Equidispersion dan Overdispersion
Asumsi penting model Poisson adalah:
\[ E(Y_i \mid X_i) = Var(Y_i \mid X_i). \]
Ukuran diagnostik yang umum digunakan adalah dispersion Pearson:
\[ \hat{\phi} = \frac{\sum r_{P,i}^2}{df}. \]
bike_dispersion_pearson <- sum(residuals(bike_poisson_fit, type = "pearson")^2) /
df.residual(bike_poisson_fit)
bike_dispersion_deviance <- bike_poisson_fit$deviance /
bike_poisson_fit$df.residual
bike_dispersion_table <- data.frame(
Ukuran = c("Dispersion Pearson", "Residual deviance / df residual"),
Nilai = c(
round(bike_dispersion_pearson, 3),
round(bike_dispersion_deviance, 3)
),
Interpretasi = c(
ifelse(
bike_dispersion_pearson < 1.5,
"Tidak ada indikasi overdispersion berat",
ifelse(bike_dispersion_pearson < 2.5, "Indikasi sedang", "Indikasi kuat")
),
ifelse(
bike_dispersion_deviance < 1.5,
"Tidak ada indikasi overdispersion berat",
ifelse(bike_dispersion_deviance < 2.5, "Indikasi sedang", "Indikasi kuat")
)
)
)
knitr::kable(
bike_dispersion_table,
caption = "Pemeriksaan overdispersion model Poisson"
)| Ukuran | Nilai | Interpretasi |
|---|---|---|
| Dispersion Pearson | 152.752 | Indikasi kuat |
| Residual deviance / df residual | 165.348 | Indikasi kuat |
Interpretasi output
Nilai dispersion sekitar 1 mendukung asumsi Poisson. Nilai yang jauh lebih besar dari 1 menunjukkan overdispersion, sehingga standard error model Poisson dapat terlalu kecil.
8.11.6 Pemeriksaan Zero Inflation
Zero inflation terjadi jika jumlah nol aktual jauh lebih banyak daripada yang diprediksi model Poisson. Pada data Bike Sharing harian, nilai nol biasanya tidak menjadi masalah karena layanan tetap memiliki peminjaman hampir setiap hari.
bike_observed_zero <- mean(bike_data$cnt == 0)
bike_expected_zero <- mean(exp(-fitted(bike_poisson_fit)))
bike_zero_table <- data.frame(
Ukuran = c("Proporsi nol aktual", "Proporsi nol harapan menurut Poisson"),
Nilai = scales::percent(c(bike_observed_zero, bike_expected_zero), accuracy = 0.01)
)
knitr::kable(
bike_zero_table,
caption = "Pemeriksaan indikasi zero inflation"
)| Ukuran | Nilai |
|---|---|
| Proporsi nol aktual | 0.00% |
| Proporsi nol harapan menurut Poisson | 0.00% |
Interpretasi output
Jika proporsi nol aktual jauh lebih besar daripada proporsi nol harapan, model zero-inflated Poisson dapat dipertimbangkan. Jika nol aktual sangat sedikit atau tidak ada, zero inflation bukan perhatian utama.
8.11.7 Exposure dan Offset
Offset diperlukan ketika setiap observasi memiliki kesempatan kejadian yang berbeda. Misalnya jumlah kecelakaan dibandingkan antar wilayah dengan jumlah penduduk berbeda, maka perlu offset \(\log(\text{populasi})\).
Pada data Bike Sharing harian, setiap baris mewakili satu hari. Karena exposure waktunya sama, yaitu satu hari, model ini tidak memakai offset. Jika data dibandingkan antar durasi yang berbeda, maka offset harus ditambahkan.
bike_offset_note <- data.frame(
Aspek = c("Unit exposure", "Offset yang digunakan", "Alasan"),
Keterangan = c(
"Satu hari",
"Tidak ada offset",
"Setiap observasi memiliki durasi pengamatan yang sama"
)
)
knitr::kable(
bike_offset_note,
caption = "Pemeriksaan kebutuhan offset"
)| Aspek | Keterangan |
|---|---|
| Unit exposure | Satu hari |
| Offset yang digunakan | Tidak ada offset |
| Alasan | Setiap observasi memiliki durasi pengamatan yang sama |
8.11.8 Pengamatan Berpengaruh
Pengamatan dengan nilai Cook’s distance tinggi dapat memiliki pengaruh besar terhadap estimasi model.
bike_cooks <- cooks.distance(bike_poisson_fit)
bike_influence_table <- data.frame(
Statistik = c(
"Cook's distance maksimum",
"Jumlah observasi dengan Cook's distance > 4/n"
),
Nilai = c(
round(max(bike_cooks), 4),
sum(bike_cooks > 4 / nrow(bike_data))
)
)
knitr::kable(
bike_influence_table,
caption = "Pemeriksaan pengamatan berpengaruh"
)| Statistik | Nilai |
|---|---|
| Cook’s distance maksimum | 9.1685 |
| Jumlah observasi dengan Cook’s distance > 4/n | 633.0000 |
Interpretasi output
Jika banyak observasi memiliki Cook’s distance tinggi, maka model perlu diperiksa lebih lanjut. Pada data deret waktu harian, hari-hari ekstrem seperti cuaca buruk, libur besar, atau anomali sistem dapat menjadi pengamatan berpengaruh.
8.11.9 Alternatif Jika Asumsi Poisson Bermasalah
Jika asumsi Poisson tidak terpenuhi, beberapa alternatif yang dapat digunakan adalah:
- Quasi-Poisson, jika masalah utama adalah overdispersion dan fokus pada standard error.
- Negative binomial, jika varians jauh lebih besar daripada mean.
- Zero-inflated Poisson, jika terlalu banyak nilai nol.
- Model dengan offset, jika exposure berbeda antar observasi.
- Model time series atau generalized additive model, jika terdapat pola waktu atau non-linearitas kuat.
bike_quasi_fit <- glm(
cnt ~ season + yr + mnth + holiday + weekday + workingday +
weathersit + temp + hum + windspeed,
data = bike_data,
family = quasipoisson(link = "log")
)
bike_quasi_table <- data.frame(
Model = c("Poisson", "Quasi-Poisson"),
`Parameter dispersi` = c(
round(summary(bike_poisson_fit)$dispersion, 3),
round(summary(bike_quasi_fit)$dispersion, 3)
)
)
knitr::kable(
bike_quasi_table,
caption = "Perbandingan parameter dispersi"
)| Model | Parameter.dispersi |
|---|---|
| Poisson | 1.000 |
| Quasi-Poisson | 152.753 |
Interpretasi output
Quasi-Poisson tidak mengubah bentuk rata-rata prediksi, tetapi menyesuaikan standard error ketika terdapat overdispersion. Oleh karena itu, quasi-Poisson sering digunakan sebagai pembanding ketika data hitung memiliki variasi lebih besar daripada yang diasumsikan model Poisson.
8.11.10 Ringkasan Pemeriksaan Asumsi Poisson
bike_assumption_summary <- data.frame(
Asumsi = c(
"Respons hitungan non-negatif",
"Observasi independen",
"Linearitas pada skala log",
"Equidispersion",
"Tidak ada zero inflation berat",
"Exposure sesuai",
"Tidak didominasi pengamatan berpengaruh"
),
`Cara pemeriksaan` = c(
"Cek minimum, integer, dan proporsi nol",
"Cek unit observasi dan potensi autokorelasi",
"Plot residual Pearson vs fitted",
"Dispersion Pearson dan deviance/df",
"Bandingkan nol aktual dan nol harapan",
"Cek apakah semua observasi memiliki durasi sama",
"Cook's distance"
),
`Catatan pada data Bike Sharing` = c(
"Sesuai",
"Perlu catatan karena data harian",
"Dilihat dari pola residual",
"Perlu diperiksa dari nilai dispersion",
"Biasanya bukan masalah utama pada data harian ini",
"Tidak memakai offset karena unitnya satu hari",
"Hari ekstrem perlu diperhatikan"
)
)
knitr::kable(
bike_assumption_summary,
caption = "Ringkasan asumsi regresi Poisson"
)| Asumsi | Cara.pemeriksaan | Catatan.pada.data.Bike.Sharing |
|---|---|---|
| Respons hitungan non-negatif | Cek minimum, integer, dan proporsi nol | Sesuai |
| Observasi independen | Cek unit observasi dan potensi autokorelasi | Perlu catatan karena data harian |
| Linearitas pada skala log | Plot residual Pearson vs fitted | Dilihat dari pola residual |
| Equidispersion | Dispersion Pearson dan deviance/df | Perlu diperiksa dari nilai dispersion |
| Tidak ada zero inflation berat | Bandingkan nol aktual dan nol harapan | Biasanya bukan masalah utama pada data harian ini |
| Exposure sesuai | Cek apakah semua observasi memiliki durasi sama | Tidak memakai offset karena unitnya satu hari |
| Tidak didominasi pengamatan berpengaruh | Cook’s distance | Hari ekstrem perlu diperhatikan |
8.12 Perbandingan Poisson dan OLS Log(Y)
Pendekatan OLS terhadap \(\log(Y)\) kadang digunakan untuk data positif yang sangat miring. Namun, pendekatan ini memiliki beberapa keterbatasan:
- Tidak cocok jika terdapat nilai \(Y = 0\), karena \(\log(0)\) tidak terdefinisi.
- Model tidak secara langsung memodelkan rata-rata hitungan.
- Prediksi pada skala log perlu ditransformasikan kembali dan dapat menimbulkan bias.
- Poisson lebih sesuai ketika respon memang berupa jumlah kejadian.
bike_lm_log <- lm(
log(cnt) ~ season + yr + mnth + holiday + weekday + workingday +
weathersit + temp + hum + windspeed,
data = bike_data
)
model_compare_count <- data.frame(
Model = c("Poisson", "OLS log(cnt)"),
Keterangan = c(
"Memodelkan rata-rata hitungan dengan link log",
"Memodelkan log jumlah peminjaman dengan regresi linear"
),
`Ukuran ringkas` = c(
paste0("AIC = ", round(AIC(bike_poisson_fit), 3)),
paste0("R-squared = ", round(summary(bike_lm_log)$r.squared, 3))
)
)
knitr::kable(
model_compare_count,
caption = "Perbandingan ringkas Poisson dan OLS log(cnt)"
)| Model | Keterangan | Ukuran.ringkas |
|---|---|---|
| Poisson | Memodelkan rata-rata hitungan dengan link log | AIC = 123694.241 |
| OLS log(cnt) | Memodelkan log jumlah peminjaman dengan regresi linear | R-squared = 0.763 |
Interpretasi output
Poisson lebih alami untuk data hitung karena menghasilkan prediksi rata-rata yang selalu positif dan langsung berada pada skala jumlah kejadian. OLS pada \(\log(Y)\) dapat dipakai sebagai pendekatan eksploratif ketika semua nilai respon positif, tetapi bukan pengganti utama untuk model data hitung.
8.13 Kesimpulan
Regresi Poisson dapat digunakan untuk memodelkan jumlah peminjaman sepeda berdasarkan musim, kalender, dan cuaca. Koefisien lebih mudah dibaca dalam bentuk IRR, yaitu perubahan multiplikatif pada rata-rata jumlah peminjaman.
Namun, data Bike Sharing kemungkinan memiliki overdispersion karena jumlah peminjaman harian sangat bervariasi. Oleh karena itu, setelah model Poisson diestimasi, pemeriksaan dispersion perlu dilakukan dan hasilnya dapat dibandingkan dengan quasi-Poisson.
8.14 Perbandingan Multinomial, Ordinal, dan Poisson
Ketiga model pada bab ini digunakan untuk respons non-kontinu, tetapi jenis responsnya berbeda. Pemilihan model harus dimulai dari skala variabel respon.
comparison_mop <- data.frame(
Aspek = c(
"Jenis respons",
"Contoh pada dokumen",
"Apakah kategori berurutan?",
"Fungsi link utama",
"Interpretasi utama",
"Asumsi khas",
"Ukuran efek"
),
`Logistik Multinomial` = c(
"Kategori nominal > 2",
"Spesies Iris",
"Tidak",
"Baseline-category logit",
"Log-odds kategori target dibanding referensi",
"IIA",
"Relative risk ratio / odds ratio relatif"
),
`Logistik Ordinal` = c(
"Kategori ordinal > 2",
"Evaluasi mobil",
"Ya",
"Cumulative logit",
"Pergeseran peluang menuju kategori lebih tinggi/rendah",
"Proportional odds",
"Odds ratio kumulatif"
),
Poisson = c(
"Hitungan non-negatif",
"Jumlah peminjaman sepeda",
"Tidak relevan",
"Log link",
"Perubahan rata-rata jumlah kejadian",
"Equidispersion",
"Incidence rate ratio"
)
)
knitr::kable(
comparison_mop,
caption = "Perbandingan model multinomial, ordinal, dan Poisson"
)| Aspek | Logistik.Multinomial | Logistik.Ordinal | Poisson |
|---|---|---|---|
| Jenis respons | Kategori nominal > 2 | Kategori ordinal > 2 | Hitungan non-negatif |
| Contoh pada dokumen | Spesies Iris | Evaluasi mobil | Jumlah peminjaman sepeda |
| Apakah kategori berurutan? | Tidak | Ya | Tidak relevan |
| Fungsi link utama | Baseline-category logit | Cumulative logit | Log link |
| Interpretasi utama | Log-odds kategori target dibanding referensi | Pergeseran peluang menuju kategori lebih tinggi/rendah | Perubahan rata-rata jumlah kejadian |
| Asumsi khas | IIA | Proportional odds | Equidispersion |
| Ukuran efek | Relative risk ratio / odds ratio relatif | Odds ratio kumulatif | Incidence rate ratio |
Interpretasi output
Multinomial dan ordinal sama-sama menangani respons kategorik dengan lebih dari dua kategori. Perbedaannya terletak pada ada atau tidaknya urutan kategori. Poisson berbeda karena responsnya berupa jumlah kejadian, bukan kategori.
8.15 Pedoman Pemilihan Model
Gunakan pertanyaan berikut untuk memilih model:
- Jika respons hanya memiliki dua kategori, gunakan regresi logistik biner.
- Jika respons memiliki lebih dari dua kategori dan tidak berurutan, gunakan regresi logistik multinomial.
- Jika respons memiliki lebih dari dua kategori dan berurutan, gunakan regresi logistik ordinal.
- Jika respons berupa jumlah kejadian atau data hitung, gunakan regresi Poisson.
- Jika model Poisson mengalami overdispersion kuat, pertimbangkan quasi-Poisson atau negative binomial.
- Jika asumsi proportional odds ordinal tidak masuk akal, pertimbangkan partial proportional odds atau model multinomial.
- Jika asumsi IIA multinomial tidak masuk akal, pertimbangkan model pilihan diskret yang lebih fleksibel.
8.16 Template Pelaporan Hasil
8.16.1 Template Model Multinomial
Dengan kategori referensi [kategori referensi], koefisien variabel [nama variabel] pada kategori [kategori target] sebesar [nilai koefisien]. Nilai \(\exp(\beta)\) sebesar [nilai] menunjukkan bahwa setiap kenaikan satu unit prediktor meningkatkan atau menurunkan odds relatif memilih kategori target dibanding kategori referensi sebesar [persentase], dengan asumsi variabel lain konstan.
8.16.2 Template Model Ordinal
Pada model ordinal, koefisien variabel [nama variabel] sebesar [nilai koefisien] dengan odds ratio
sebesar [nilai OR]. Jika menggunakan
MASS::polr(), koefisien positif menunjukkan kecenderungan
menuju kategori respons yang lebih tinggi, dengan asumsi variabel lain
konstan. Interpretasi ini tetap harus mempertimbangkan asumsi
proportional odds.
8.16.3 Template Model Poisson
Koefisien variabel [nama variabel] sebesar [nilai koefisien] dengan \(IRR=\exp(\beta)\) sebesar [nilai IRR]. Hal ini menunjukkan bahwa setiap kenaikan satu unit prediktor meningkatkan atau menurunkan rata-rata jumlah kejadian sebesar [persentase], dengan asumsi variabel lain dan exposure konstan.
8.17 Kesalahan Umum dalam Interpretasi
Beberapa kesalahan yang perlu dihindari adalah:
- Menganggap koefisien logistik sama seperti koefisien regresi linear.
- Mengabaikan kategori referensi pada model multinomial.
- Memakai model multinomial untuk respons ordinal tanpa alasan.
- Memakai model ordinal tanpa memeriksa proportional odds.
- Mengabaikan overdispersion pada regresi Poisson.
- Menggunakan OLS pada \(\log(Y)\) tanpa memperhatikan nilai nol dan target model.
- Menyimpulkan kausalitas dari data observasional tanpa desain kausal.
- Hanya melaporkan p-value tanpa ukuran efek, interval, diagnostik, dan visualisasi prediksi.
9 Referensi
Agresti, A. (2013). Introduction to Categorical Data Analysis. Wiley.
Dobson, A. J. (2002). An Introduction to Generalized Linear Models. Chapman & Hall/CRC.
McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models. Chapman & Hall.