Artikel ini menyajikan pembahasan analisis data kategorik yang disusun secara bertahap melalui beberapa studi kasus. Setiap tahapan analisis dilengkapi dengan konsep pendukung, sehingga tidak hanya menekankan hasil, tetapi juga pemahaman terhadap metode yang digunakan. Oleh karena itu, artikel ini diharapkan dapat memberikan gambaran yang sistematis mengenai analisis data kategorik beserta landasan konseptualnya.

1 Tugas 1: Analisis Data Kategori

1.1 Pendahuluan

1.1.1 Pengertian Analisis Data Kategori

Analisis data kategori merupakan metode statistik yang digunakan untuk menganalisis data yang berbentuk kategori atau klasifikasi, yaitu variabel yang memiliki skala pengukuran berupa sekumpulan kategori yang digunakan untuk mengklasifikasikan suatu objek, individu, atau kejadian ke dalam kelompok tertentu. Sebagai contoh, pandangan politik dapat dikategorikan sebagai liberal, moderat, atau konservatif. Variabel kategorik dapat dibedakan berdasarkan skala pengukurannya, yaitu nominal dan ordinal, serta berdasarkan jumlah kategorinya, yaitu biner (dikotomik) dan multikategori. Variabel nominal merupakan variabel kategori yang tidak memiliki urutan tertentu, sedangkan variabel ordinal memiliki urutan atau tingkatan antar kategori. Sementara itu, variabel biner hanya memiliki dua kategori, seperti ya dan tidak, sedangkan variabel multikategori memiliki lebih dari dua kategori, seperti pilihan tempat tinggal yang dapat berupa rumah, kondominium, atau apartemen.

1.1.2 Karakteristik Variabel Kategori

  1. Nilainya berupa kategori atau label, bukan angka yang menunjukkan besaran kuantitatif.
  2. Digunakan untuk mengklasifikasikan objek atau individu ke dalam kelompok tertentu.
  3. Dapat memiliki dua kategori (biner/dikotomik) atau lebih dari dua kategori (multikategori).
  4. Beberapa variabel memiliki urutan kategori (ordinal), sedangkan yang lain tidak memiliki urutan (nominal).

1.1.3 Contoh Penerapan Analisis Data Kategori

Analisis data kategori banyak digunakan dalam berbagai bidang penelitian. Beberapa contoh penerapannya adalah sebagai berikut:

  • Ilmu sosial: digunakan untuk mengukur sikap dan opini, misalnya klasifikasi pandangan politik seperti liberal, moderat, dan konservatif.
  • Ilmu kesehatan: digunakan untuk mengukur respons pasien, seperti apakah pasien selamat setelah operasi (ya/tidak), tingkat keparahan cedera (tidak ada, ringan, sedang, berat), atau tahap penyakit (awal, lanjut).
  • Ilmu perilaku: digunakan untuk mengklasifikasikan jenis gangguan mental seperti skizofrenia, depresi, dan neurosis.
  • Kesehatan masyarakat: digunakan untuk mengetahui perubahan perilaku masyarakat, misalnya apakah kesadaran tentang AIDS meningkatkan penggunaan kondom (ya/tidak).
  • Zoologi: digunakan untuk mengelompokkan jenis makanan utama hewan, misalnya makanan utama buaya berupa ikan, invertebrata, atau reptil.
  • Pendidikan: digunakan untuk mengklasifikasikan jawaban mahasiswa pada ujian menjadi benar atau salah.
  • Pemasaran: digunakan untuk menganalisis preferensi konsumen terhadap suatu produk, misalnya pilihan Merek A, Merek B, atau Merek C.

1.2 Tabel Kontingensi

1.2.1 Definisi Tabel Kontingensi

Tabel kontingensi adalah tabel yang menyajikan distribusi frekuensi dari dua atau lebih variabel kategorik secara simultan, sehingga setiap sel menunjukkan jumlah observasi pada setiap kombinasi kategori dari variabel-variabel tersebut. Tabel ini digunakan untuk melihat pola atau hubungan antar variabel kategorik.

1.2.2 Struktur Tabel Kontingensi

Struktur tabel kontingensi untuk dua variabel kategorik disajikan dalam bentuk tabel persegi panjang dengan \(I\) baris yang mewakili kategori variabel \(X\) dan \(J\) kolom yang mewakili variabel \(Y\). Setiap sel dalam tabel menunjukkan kombinasi kategori dari kedua variabel tersebut, sehingga terdapat \(I \times J\) kemungkinan kombinasi hasil. Berikut disajikan contoh tabel kontingensi dua arah yang menunjukkan kombinasi kategori dari dua variabel kategorik.

Tabel 1.1 Struktur Tabel Kontingensi 2×2

\(Y=1\) \(Y=0\) Total
\(X=1\) \(n_{11}\) \(n_{12}\) \(n_{1\cdot}\)
\(X=0\) \(n_{21}\) \(n_{22}\) \(n_{2\cdot}\)
Total \(n_{\cdot1}\) \(n_{\cdot2}\) \(n\)

dengan:

\(n_{ij} = \text{jumlah observasi pada kategori } X=i \text{ dan } Y=j\)

\(n = \sum_{i=1}^{2}\sum_{j=1}^{2} n_{ij}\)

1.2.3 Joint Distribution

Distribusi peluang bersama dinyatakan dengan \(\pi_{ij}\), yaitu peluang bahwa variabel \(X\) berada pada kategori ke-\(i\) dan variabel \(Y\) berada pada kategori ke-\(j\). Dalam praktiknya, nilai peluang tersebut dapat diestimasi menggunakan proporsi frekuensi pada setiap sel tabel kontingensi, yaitu

\[ \pi_{ij} = \frac{n_{ij}}{n} \]

dengan \(n_{ij}\) menyatakan jumlah observasi pada sel ke-\((i,j)\) dan \(n\) menyatakan jumlah total observasi.

1.2.4 Marginal Distribution

Distribusi peluang marginal merupakan distribusi peluang dari masing-masing variabel secara terpisah tanpa memperhatikan variabel lainnya. Pada tabel kontingensi, distribusi marginal diperoleh dengan menjumlahkan peluang pada setiap baris atau kolom.

Distribusi marginal untuk variabel \(X\) dinyatakan sebagai

\[ \pi_{i.} = \sum_{j=1}^{2} \pi_{ij}, \]

sedangkan distribusi marginal untuk variabel \(Y\) dinyatakan sebagai

\[ \pi_{.j} = \sum_{i=1}^{2} \pi_{ij}. \]

Sebagai contoh, berdasarkan Tabel 1.1, peluang marginal dapat diperoleh dari proporsi frekuensi pada setiap baris atau kolom tabel kontingensi. Karena peluang bersama diestimasi dengan

\[ \pi_{ij} = \frac{n_{ij}}{n}, \]

maka peluang marginal untuk kategori pertama variabel \(X\) diperoleh dengan

\[ \pi_{1.} = \pi_{11} + \pi_{12}. \]

Sedangkan peluang marginal untuk kategori pertama variabel \(Y\) diperoleh dengan

\[ \pi_{.1} = \pi_{11} + \pi_{21}. \]

1.2.5 Conditional Probability

Distribusi peluang bersyarat (conditional probability) menyatakan peluang suatu kategori dari satu variabel dengan syarat bahwa kategori variabel lainnya telah diketahui. Dalam konteks tabel kontingensi, peluang bersyarat dihitung dengan membandingkan peluang bersama dengan peluang marginal.

\[ P(Y=j \mid X=i) = \frac{\pi_{ij}}{\pi_{i.}} \]

Sebaliknya, peluang bersyarat variabel \(X\) pada kategori ke-\(i\) dengan syarat variabel \(Y\) berada pada kategori ke-\(j\) dinyatakan sebagai

\[ P(X=i \mid Y=j) = \frac{\pi_{ij}}{\pi_{.j}} \]

1.3 Ukuran Asosiasi

Ukuran asosiasi digunakan untuk mengukur kekuatan hubungan antara dua variabel dalam tabel kontingensi.

1.3.1 Odds

Odds merupakan perbandingan antara probabilitas terjadinya suatu kejadian dengan probabilitas tidak terjadinya kejadian tersebut. Dalam konteks tabel kontingensi, odds dapat dinyatakan sebagai perbandingan antara peluang suatu kejadian dan peluang komplemennya.

Secara umum, jika probabilitas suatu kejadian dinyatakan dengan \(\pi\), maka odds didefinisikan sebagai:

\[ \text{Odds} = \frac{\pi}{1-\pi}. \]

1.3.2 Odds Ratio

Odds Ratio (OR) digunakan untuk membandingkan odds antara dua kelompok pada tabel kontingensi \(2 \times 2\). Misalkan probabilitas kejadian pada kelompok 1 adalah \(\pi_1\) dan pada kelompok 2 adalah \(\pi_2\). Maka odds pada masing-masing kelompok adalah:

\[ \text{Odds}_1 = \frac{\pi_1}{1-\pi_1}, \qquad \text{Odds}_2 = \frac{\pi_2}{1-\pi_2}. \]

Odds Ratio didefinisikan sebagai

\[ OR = \frac{\text{Odds}_1}{\text{Odds}_2} = \frac{\pi_1/(1-\pi_1)}{\pi_2/(1-\pi_2)}. \]

Pada tabel kontingensi \(2 \times 2\) dengan probabilitas sel \(\pi_{ij}\), Odds Ratio juga dapat dituliskan sebagai

\[ OR = \frac{\pi_{11}\pi_{22}}{\pi_{12}\pi_{21}}. \]

Interpretasi nilai Odds Ratio adalah sebagai berikut:

  • \(OR = 1\) menunjukkan tidak terdapat asosiasi antara kedua variabel.
  • \(OR > 1\) menunjukkan kejadian lebih mungkin terjadi pada kelompok pertama.
  • \(OR < 1\) menunjukkan kejadian lebih kecil kemungkinannya terjadi pada kelompok pertama.

1.3.3 Relative Risk

Relative Risk (RR) merupakan ukuran yang digunakan untuk membandingkan probabilitas terjadinya suatu kejadian pada dua kondisi atau kelompok yang berbeda. Secara umum, Relative Risk didefinisikan sebagai rasio antara dua probabilitas bersyarat.

\[ RR = \frac{P(Y=1 \mid X=1)}{P(Y=1 \mid X=0)} = \frac{\pi_{11}/\pi_{1.}}{\pi_{21}/\pi_{2.}} \]

Nilai Relative Risk diinterpretasikan sebagai berikut:

  • \(RR = 1\) menunjukkan bahwa probabilitas kejadian sama pada kedua kondisi.
  • \(RR > 1\) menunjukkan bahwa probabilitas kejadian lebih besar pada kondisi pertama.
  • \(RR < 1\) menunjukkan bahwa probabilitas kejadian lebih kecil pada kondisi pertama.

1.3.4 Risk Difference

Risk difference (RD) adalah ukuran yang menyatakan selisih probabilitas terjadinya suatu kejadian antara dua kelompok yang dibandingkan. Ukuran ini menunjukkan seberapa besar perbedaan risiko kejadian pada kelompok yang terpapar suatu faktor dibandingkan dengan kelompok yang tidak terpapar.

Secara matematis, risk difference dinyatakan sebagai: \[ RD = \frac{\pi_{11}}{\pi_{1.}} - \frac{\pi_{21}}{\pi_{2.}} \]

  • \(RD = 0\) berarti tidak ada perbedaan risiko antara kedua kelompok
  • \(RD > 0\) berarti risiko kejadian pada kelompok 1 lebih tinggi dibanding kelompok 2
  • \(RD < 0\) berarti risiko kejadian pada kelompok 1 lebih rendah dibanding kelompok 2

1.4 Inferensi Tabel Kontingensi

1.4.1 Estimasi

1.4.1.1 Estimasi Titik

Estimasi titik adalah metode dalam inferensi statistik yang digunakan untuk menaksir nilai suatu parameter populasi menggunakan satu nilai yang diperoleh dari data sampel. Estimasi titik untuk proporsi populasi \(\pi\) diberikan oleh:

\[ \hat{p} = \frac{x}{n} \]

dengan:

  • \(\hat{p}\) adalah estimasi titik proporsi
  • \(x\) adalah jumlah kejadian yang diamati dalam sampel
  • \(n\) adalah ukuran sampel

1.4.1.2 Estimasi Interval

Selain estimasi titik, sering kali diperlukan ukuran yang menggambarkan ketidakpastian dari estimasi tersebut. Estimasi interval memberikan suatu rentang nilai yang kemungkinan besar memuat parameter populasi yang sebenarnya. Rentang ini disebut interval kepercayaan (confidence interval).

Untuk ukuran sampel yang cukup besar, interval kepercayaan untuk proporsi populasi dapat didekati menggunakan distribusi normal sebagai berikut:

\[ \hat{p} \pm z_{\alpha/2} \sqrt{\frac{\hat{p}(1-\hat{p})}{n}} \]

dengan:

  • \(z_{\alpha/2}\) adalah nilai kritis distribusi normal standar
  • \(n\) adalah ukuran sampel

1.5 Uji Hipotesis

1.5.1 Uji Proporsi

Uji proporsi digunakan untuk mengetahui apakah terdapat perbedaan proporsi antara dua kelompok dalam populasi. Pada analisis tabel kontingensi 2×2, pengujian ini dilakukan dengan membandingkan proporsi kejadian pada dua kelompok yang berbeda. Hipotesis yang diuji adalah sebagai berikut:

  • Hipotesis nol (\(H_{0}\)): \(\pi_1 = \pi_2\) (Tidak terdapat perbedaan proporsi antara kedua kelompok)
  • Hipotesis alternatif (\(H_{1}\)): \(\pi_1 \ne \pi_2\) (Terdapat perbedaan proporsi antara kedua kelompok)

Estimasi proporsi sampel untuk masing-masing kelompok dinyatakan sebagai:

\[ \hat{p}_1 = \frac{x_1}{n_1}, \qquad \hat{p}_2 = \frac{x_2}{n_2} \]

dengan:

  • \(x_1\) = jumlah kejadian pada kelompok pertama
  • \(x_2\) = jumlah kejadian pada kelompok kedua
  • \(n_1\) = ukuran sampel pada kelompok pertama
  • \(n_2\) = ukuran sampel pada kelompok kedua

Statistik uji yang digunakan adalah statistik uji Z yang dirumuskan sebagai:

\[ Z = \frac{\hat{p}_1 - \hat{p}_2}{\sqrt{\hat{p}(1-\hat{p})\left(\frac{1}{n_1} + \frac{1}{n_2}\right)}} \]

dengan \(\hat{p}\) merupakan proporsi gabungan (pooled proportion) yang dihitung sebagai:

\[ \hat{p} = \frac{x_1 + x_2}{n_1 + n_2} \] Keputusan pengujian dilakukan dengan membandingkan nilai statistik uji dengan distribusi normal standar atau dengan menggunakan nilai p-value. Jika nilai p-value < atau Z > nilai kritis pada tertentu, maka hipotesis nol ditolak sehingga dapat disimpulkan bahwa terdapat perbedaan proporsi yang signifikan antara kedua kelompok.

1.5.2 Uji Asosiasi

Uji asosiasi digunakan untuk mengetahui apakah terdapat hubungan antara dua variabel kategorik pada tabel kontingensi. Pada tabel kontingensi 2×2, hubungan antar variabel dapat diukur menggunakan beberapa ukuran asosiasi seperti Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR). Hipotesis yang diujikan adalah sebagai berikut:

  • Hipotesis nol (\(H_{0}\)): Tidak terdapat asosiasi antara kedua variabel
  • Hipotesis alternatif (\(H_{1}\)): Terdapat asosiasi antara kedua variabel

1.5.2.1 Risk Difference

Risk Difference merupakan selisih antara dua proporsi kejadian.

\[ RD = \frac{\pi_{11}}{\pi_{1.}} - \frac{\pi_{21}}{\pi_{2.}} \]

Standard error untuk Risk Difference adalah:

\[ 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 = \frac{RD}{SE(RD)} \]

1.5.2.2 Relative Risk (RR)

Relative Risk merupakan rasio antara dua proporsi kejadian.

\[ RR = \frac{\pi_{11}/\pi_{1.}}{\pi_{21}/\pi_{2.}} \]

Standard error dihitung menggunakan transformasi logaritma:

\[ SE(\ln(RR)) = \sqrt{\frac{1}{n_{11}} - \frac{1}{n_{1.}} + \frac{1}{n_{21}} - \frac{1}{n_{2.}}} \]

Statistik uji:

\[ Z = \frac{\ln(RR)}{SE(\ln(RR))} \]

1.5.2.3 Odds Ratio (OR)

Odds Ratio merupakan rasio antara odds kejadian pada dua kelompok.

\[ OR = \frac{\pi_{11}\pi_{22}}{\pi_{12}\pi_{21}}. \]

Standard error dihitung dengan:

\[ SE(\ln(OR)) = \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}} \]

Statistik uji:

\[ Z = \frac{\ln(OR)}{SE(\ln(OR))} \]

1.5.3 Uji Independensi

Uji independensi digunakan untuk mengetahui apakah terdapat hubungan signifikan antara dua variabel kategorik pada tabel kontingensi. Hipotesis yang diuji adalah:

  • Hipotesis nol (\(H_0\)): Variabel bersifat independen (Terdapat hubungan antara kedua variabel)
  • Hipotesis alternatif (\(H_1\)): Variabel tidak independen (Tidak terdapat hubungan antara kedua variabel)

1.5.3.1 Uji Chi-Square

Uji Chi-Square membandingkan frekuensi observasi dengan frekuensi harapan jika kedua variabel bersifat independen. Statistik uji Chi-Square diperoleh sebagai berikut:

\[ \chi^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \]

dengan:

  • \(O_{ij}\) adalah observed frequency (frekuensi observasi) pada sel baris ke-\(i\) dan kolom ke-\(j\) dalam tabel kontingensi.
  • \(E_{ij}\) adalah expected frequency (frekuensi harapan) pada sel baris ke-\(i\) dan kolom ke-\(j\)

\[ E_{ij} = \frac{n_{i.} n_{.j}}{n} \]

1.5.3.2 Partisi Chi-Square

Partisi Chi-Square digunakan untuk menentukan kontribusi tiap sel atau sub-tabel terhadap statistik Chi-Square total. Metode ini membantu mengidentifikasi sel atau kategori mana yang paling memengaruhi asosiasi.

Langkah-Langkah Partisi Chi-Square:

  1. Hitung Chi-Square Total dari tabel kontingensi \(I \times J\)

\[\chi^2_\text{total} = \sum_{i=1}^{I} \sum_{j=1}^{J} \frac{(O_{ij} - E_{ij})^2}{E_{ij}} ; E_{ij} = \frac{n_{i.} \cdot n_{.j}}{n}\]

  1. Hitung Residual Pearson per Sel

\[R_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}}\]

  1. Bagi Tabel menjadi Sub-Tabel 2×2
  • Tabel besar dibagi menjadi beberapa sub-tabel 2×2
  • Hitung \(\chi^2\) untuk tiap sub-tabel
  • Kontribusi tiap sub-tabel dapat dijumlahkan sehingga mendekati \(\chi^2_\text{total}\)
  1. Interpretasi
  • Nilai \(R_{ij}\) besar berarti sel berkontribusi signifikan pada Chi-Square total
  • Sub-tabel dengan Chi-Square tinggi menunjukkan area yang paling memengaruhi asosiasi
  • Berguna untuk mendeteksi kategori yang tidak independen atau outlier

1.5.3.3 Uji Likelihood Ratio

Likelihood Ratio atau G² adalah alternatif uji Chi-Square menggunakan log-likelihood. Statistik uji yang digunakan adalah:

\[ G^2 = 2 \sum O_{ij} \ln \frac{O_{ij}}{E_{ij}} \] Nilai kritisnya didaapatkan dari distribusi Chi-square dengan \(df = (baris-1)*(kolom-1)\)

1.5.3.4 Fisher Exact Test

Uji Fisher digunakan untuk menguji asosiasi antara dua variabel kategorik ketika jumlah sampel kecil atau ada sel yang frekuensinya <5. Berbeda dengan chi-square, Fisher menghitung probabilitas exact dari tabel kontingensi. Probabilitas tabel tertentu dihitung dengan rumus hypergeometric:

\[P(X = x) = \frac{\binom{K}{x} \cdot \binom{N-K}{n-x}}{\binom{N}{n}}\] dengan:

  • \(N\) adalah total objek dalam populasi
  • \(K\) adalah jumlah objek dalam kategori tertentu
  • \(n\) adalah jumlah sampel yang diambil
  • \(x\) adalah jumlah objek kategori tertentu yang diamati dalam sampel

p-value uji Fisher didapatkan dari total probabilitas semua tabel yang sama atau lebih ekstrem daripada tabel yang diamati

1.6 Analisis Residual dan Deteksi Outlier

Analisis residual digunakan untuk menilai sejauh mana setiap sel dalam tabel kontingensi menyimpang dari nilai yang diharapkan berdasarkan distribusi independen. Residual standar (\(R_{ij}\)) menunjukkan besarnya penyimpangan tiap sel relatif terhadap ekspektasi dan variabilitasnya. Residual standar yang besar menunjukkan kontribusi signifikan terhadap statistik Chi-Square total dan dapat membantu mendeteksi sel-sel yang berperilaku tidak biasa atau outlier. Dengan mengetahui sel yang menyimpang, kita dapat meninjau data lebih lanjut, memvalidasi kesalahan pengukuran, atau memahami fenomena yang tidak sesuai dengan pola umum.

Langkah-langkah analisis residual dan deteksi outlier adalah sebagai berikut:

1. Hitung Chi-Square Total

\[ \chi^2_\text{total} = \sum_{i=1}^{I} \sum_{j=1}^{J} \frac{(O_{ij} - E_{ij})^2}{E_{ij}} ; E_{ij} = \frac{n_{i.} \cdot n_{.j}}{n} \]

2. Hitung Residual Standar (Standardized Residual) per Sel
\[ R_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}} \]

3. Analisis Residual Standar

  • Residual standar dapat digunakan untuk mengidentifikasi sel yang menyimpang dari ekspektasi. Biasanya, residual \(|R_{ij}| > 2\) atau \(|R_{ij}| > 3\) dianggap signifikan:
  • \(|R_{ij}| > 2\) → kontribusi moderat terhadap Chi-Square
  • \(|R_{ij}| > 3\) → kontribusi besar, patut dicurigai

4. Deteksi Outlier

  1. Hitung \(R_{ij}\) untuk setiap sel tabel kontingensi.
  2. Tentukan threshold, misal \(|R_{ij}| > 3\).
  3. Tandai sel yang melewati threshold sebagai outlier.
  4. Analisis konteks data untuk memutuskan apakah outlier ini valid atau hasil kesalahan.

5. Interpretasi
- \(R_{ij} \approx 0\) → Observasi sesuai ekspektasi, tidak ada hubungan kuat antar kategori. - \(R_{ij}\) positif besar → Observasi lebih tinggi dari ekspektasi → hubungan positif antar kategori.
- \(R_{ij}\) negatif besar → Observasi lebih rendah dari ekspektasi → hubungan negatif atau tidak ada asosiasi. - Outlier harus diperiksa apakah representatif atau error pengukuran

1.7 Contoh Kasus Tabel Kontingensi 2x2

Sebuah toko ingin mengetahui apakah promosi diskon mempengaruhi keputusan konsumen untuk membeli produk. Dari survei terhadap 200 konsumen diperoleh data yang disajikan dalam tabel kontingensi berikut.

Tabel 1.2 Tabel Kontingensi Keputusan Konsumen dalam Membeli Produk

Membeli Tidak Membeli Total
Promosi 70 30 100
Tidak Promosi 40 60 100
Total 110 90 200

1.7.1 Pengerjaan Manual

1.7.1.1 Menghitung Peluang Bersyarat

  • Peluang membeli dengan promosi:

\[ P(Y=1|X=1)=\frac{n_{11}}{n_{11}+n_{12}} \]

\[ P(Y=1|X=1)=\frac{70}{70+30}=\frac{70}{100}=0.7 \]

  • Peluang membeli tanpa promosi:

\[ P(Y=1|X=0)=\frac{n_{21}}{n_{21}+n_{22}} \]

\[ P(Y=1|X=0)=\frac{40}{40+60}=\frac{40}{100}=0.4 \]

1.7.1.2 Menghitung Odds

  • Odds membeli ketika ada promosi:

\[ \text{Odds}_1=\frac{P(Y=1|X=1)}{P(Y=0|X=1)} \]

\[ \text{Odds}_1=\frac{70/100}{30/100}=\frac{70}{30}=2.33 \]

  • Odds membeli tanpa promosi:

\[ \text{Odds}_0=\frac{P(Y=1|X=0)}{P(Y=0|X=0)} \]

\[ \text{Odds}_0=\frac{40/100}{60/100}=\frac{40}{60}=0.67 \]

1.7.1.3 Menghitung Odds Ratio

Odds Ratio didefinisikan sebagai:

\[ OR=\frac{\text{Odds}_1}{\text{Odds}_0} \]

atau dapat dihitung langsung dari tabel:

\[ OR=\frac{n_{11}n_{22}}{n_{12}n_{21}} \]

Substitusi nilai:

\[ OR=\frac{70\times60}{30\times40} \]

\[ OR=\frac{4200}{1200}=3.5 \]

1.7.2 Analisis Menggunakan R

# ================================
# Contoh Kasus 2x2
# ================================

# Membuat tabel kontingensi
data <- matrix(c(70, 30,
                 40, 60),
               nrow = 2,
               byrow = TRUE)

rownames(data) <- c("Promosi","Tidak_Promosi")
colnames(data) <- c("Membeli","Tidak_Membeli")
data
##               Membeli Tidak_Membeli
## Promosi            70            30
## Tidak_Promosi      40            60
# ================================
# 1. Menghitung Odds Ratio (OR), Risk Ratio (RR), dan Risk Difference (RD)
# ================================

# Odds Ratio
OR <- (data[1,1] * data[2,2]) / (data[1,2] * data[2,1]) 
OR
## [1] 3.5
# Risiko membeli di masing-masing grup
risk_promosi <- data[1,1] / sum(data[1,])
risk_tidak <- data[2,1] / sum(data[2,])

# Risk Ratio
RR <- risk_promosi / risk_tidak
RR
## [1] 1.75
# Risk Difference
RD <- risk_promosi - risk_tidak
RD
## [1] 0.3
# ================================
# 2. Uji Chi-Square
# ================================
uji_chi <- chisq.test(data)
uji_chi
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data
## X-squared = 16.99, df = 1, p-value = 3.758e-05
# Expected frequency
uji_chi$expected
##               Membeli Tidak_Membeli
## Promosi            55            45
## Tidak_Promosi      55            45
# ================================
# 3. Hitung Residual Standar (Standardized Residual)
# ================================
residual <- (data - uji_chi$expected) / sqrt(uji_chi$expected)
residual
##               Membeli Tidak_Membeli
## Promosi        2.0226     -2.236068
## Tidak_Promosi -2.0226      2.236068
# ================================
# 4. Interpretasi Residual
# ================================
# Aturan sederhana:
# |R_ij| ≈ 0 → Observasi sesuai ekspektasi
# R_ij positif besar → Observasi lebih tinggi dari ekspektasi
# R_ij negatif besar → Observasi lebih rendah dari ekspektasi

interpretasi <- apply(residual, c(1,2), function(x){
  if(abs(x) < 2){
    "Sesuai ekspektasi"
  } else if(x >= 2){
    "Positif signifikan"
  } else {
    "Negatif signifikan"
  }
})

interpretasi
##               Membeli              Tidak_Membeli       
## Promosi       "Positif signifikan" "Negatif signifikan"
## Tidak_Promosi "Negatif signifikan" "Positif signifikan"

1.7.3 Interpretasi Hasil

1.7.3.1 Interpretasi Statistik

Hasil analisis menunjukkan bahwa konsumen yang menerima promosi diskon memiliki peluang lebih besar untuk membeli produk dibandingkan yang tidak menerima promosi. Nilai Odds Ratio (OR = 3.5) lebih besar dari 1, menunjukkan hubungan positif antara promosi dan keputusan membeli, artinya konsumen yang menerima promosi sekitar 3.5 kali lebih mungkin membeli. Risk Ratio (RR = 1.75) mendukung hal ini dengan menunjukkan konsumen pada kelompok promosi sekitar 1.75 kali lebih mungkin melakukan pembelian dibandingkan kelompok tanpa promosi, sedangkan Risk Difference (RD = 0.3) menunjukkan selisih probabilitas membeli sebesar 0.3 atau sekitar 30% lebih banyak konsumen membeli pada kelompok promosi.

Uji Chi-Square digunakan untuk menilai apakah hubungan tersebut signifikan secara statistik. Hasil uji Chi-Square menunjukkan p-value sebesar 3.7579211^{-5}. Karena p-value < 0.05, maka dapat disimpulkan bahwa terdapat hubungan yang signifikan secara statistik antara promosi dan keputusan membeli. Pemeriksaan expected frequency juga menunjukkan bahwa seluruh nilai harapan lebih besar dari 5, sehingga asumsi uji Chi-Square terpenuhi dan hasil pengujian dapat dianggap valid.

Selain itu, residual standar per sel menunjukkan sel-sel tertentu berkontribusi signifikan terhadap Chi-Square, misal:
- Membeli-Promosi: positif signifikan
- Tidak Membeli-Promosi: negatif signifikan
- Membeli-Tidak Promosi: negatif signifikan
- Tidak Membeli-Tidak Promosi: positif signifikan

Hal ini menandakan adanya penyimpangan yang signifikan yang mendukung hubungan antara variabel.

1.7.3.2 Interpretasi Substantif

Secara substantif, hasil analisis menunjukkan bahwa promosi memiliki pengaruh terhadap keputusan konsumen dalam membeli produk. Berdasarkan hasil analisis, konsumen yang menerima promosi diskon memiliki peluang sekitar 3.5 kali lebih besar untuk melakukan pembelian dibandingkan konsumen yang tidak menerima promosi.

Berdasarkan tabel kontingensi, dari 100 konsumen yang menerima promosi, sekitar 70% melakukan pembelian, sedangkan pada kelompok yang tidak menerima promosi hanya sekitar 40% yang melakukan pembelian. Hasil ini mengindikasikan bahwa pemberian promosi diskon dapat meningkatkan kemungkinan konsumen membeli produk di toko tersebut.

1.8 Implementasi Inferensi Tabel Kontingensi

1.8.1 Uji Proporsi

Misalkan diberikan data sebagai berikut:

Tabel 1.3 Tabel Kontingensi Kejadian pada Dua Grup

Kejadian (+) Tidak Kejadian (-) Total
Grup 1 50 30 80
Grup 2 30 50 80
Total 80 80 160

Hipotesis:

  • Hipotesis nol (\(H_{0}\)): \(\hat{p}_1 = \hat{p}_2\) (Tidak terdapat perbedaan proporsi antara kedua kelompok)
  • Hipotesis alternatif (\(H_{1}\)): \(\hat{p}_1 \ne \hat{p}_2\) (Terdapat perbedaan proporsi antara kedua kelompok)

Perhitungan Manual:

Langkah 1: Hitung Proporsi Sampel

\[ \hat{p}_1 = \frac{50}{80} = 0.625, \quad \hat{p}_2 = \frac{30}{80} = 0.375 \]

Langkah 2: Hitung Proporsi Gabungan (Pooled Proportion)

\[ \hat{p} = \frac{50 + 30}{80 + 80} = \frac{80}{160} = 0.50 \]

Langkah 3: Hitung Statistik Uji Z

\[ Z = \frac{\hat{p}_1 - \hat{p}_2}{\sqrt{\hat{p}(1 - \hat{p})\left(\frac{1}{n_1} + \frac{1}{n_2}\right)}} \]

\[ 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 \]


Perhitungan Menggunakan R:

# Pastikan variabel data_matrix terdefinisi sebelum digunakan
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"))
print(data)
##         Kejadian
## Terpapar Ya Tidak
##    Ya    50    30
##    Tidak 30    50
# Uji Proporsi dengan variabel yang benar
prop_test <- prop.test(x = c(data[1,1], data[2,1]), 
                       n = c(sum(data[1,]), sum(data[2,])))
print(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.375

Interpretasi: Karena nilai \(Z = 3.16\) lebih besar dari nilai kritis \(Z_{0.05} = 1.96\), maka \(H_0\) ditolak. Artinya, terdapat perbedaan yang signifikan antara dua proporsi.

1.8.2 Uji Asosiasi

Misalkan diberikan data sebagai berikut:

Tabel 1.4 Tabel Kontingensi Kejadian pada Dua Grup

Kejadian (+) Tidak Kejadian (-) Total
Grup 1 50 30 80
Grup 2 30 50 80
Total 80 80 160

Hipotesis:

  • Hipotesis nol (\(H_{0}\)): Tidak terdapat asosiasi antara kedua variabel
  • Hipotesis alternatif (\(H_{1}\)): Terdapat asosiasi antara kedua variabel

Perhitungan Manual:

Misalkan:

\[ \hat{p}_1 = \frac{50}{80} = 0.625, \quad \hat{p}_2 = \frac{30}{80} = 0.375 \]


Risk Difference:

\[RD = \hat{p}_1 - \hat{p}_2 = 0.625 - 0.375 = 0.25\] \[SE(RD) = \sqrt{\frac{0.625(0.375)}{80} + \frac{0.375(0.625)}{80}} = 0.0765\] \[Z_{RD} = \frac{0.25}{0.0765} = 3.27\]


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}} = 0.1683\] \[Z_{RR} = \frac{\ln(1.67)}{0.1683} = 3.03\]


Odds Ratio:

\[OR = \frac{50 \times 50}{30 \times 30} = 2.78\] \[SE(\ln OR) = \sqrt{\frac{1}{50} + \frac{1}{30} + \frac{1}{30} + \frac{1}{50}} = 0.3266\] \[Z_{OR} = \frac{\ln(2.78)}{0.3266} = 3.12\]


Perhitungan Menggunakan R:

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

# Hasil
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.128155

Kesimpulan:

  • Risk Difference (RD) mengukur perbedaan risiko absolut.
  • Relative Risk (RR) membandingkan kemungkinan kejadian antara dua kelompok.
  • Odds Ratio (OR) membandingkan peluang kejadian antara dua kelompok.
  • Standard error dan statistik uji Z digunakan untuk menilai signifikansi statistik masing-masing ukuran asosiasi.

1.8.3 Uji Independensi

Hipotesis:

  • Hipotesis nol (\(H_0\)): Variabel bersifat independen (Terdapat hubungan antara kedua variabel)
  • Hipotesis alternatif (\(H_1\)): Variabel tidak independen (Tidak terdapat hubungan antara kedua variabel)

1.8.3.1 Uji Chi-Square

Misalkan diberikan data sebagai berikut:

Tabel 1.6 Tabel Kontingensi Kejadian Terpapar dan Tidak Terpapar

Kejadian Ya Kejadian Tidak Total
Terpapar 30 10 40
Tidak Terpapar 15 45 60
Total 45 55 100

Perhitungan Manual:

Hitung nilai yang diharapkan E:

\[ 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 \]

Lalu, hitung \(\chi^2\):

\[ \chi^2 = \frac{(30 - 18)^2}{18} + \frac{(10 - 22)^2}{22} + \frac{(15 - 27)^2}{27} + \frac{(45 - 33)^2}{33} \]

\[ = \frac{144}{18} + \frac{144}{22} + \frac{144}{27} + \frac{144}{33} \]

\[ = 8 + 6.55 + 5.33 + 4.36 = 24.24 \]

Dengan derajat kebebasan (df):

\[ df = (2 - 1)(2 - 1) = 1 \]

Perhitungan Menggunakan R:

# Contoh Data
set.seed(123)
data <- matrix(c(30, 10, 15, 45), nrow = 2, byrow = TRUE)
dimnames(data) <- list("Terpapar" = c("Ya", "Tidak"), "Kejadian" = c("Ya", "Tidak"))
print(data)
##         Kejadian
## Terpapar Ya Tidak
##    Ya    30    10
##    Tidak 15    45
# Uji Chi-Square
chisq_test <- chisq.test(data)
print(chisq_test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data
## X-squared = 22.264, df = 1, p-value = 2.376e-06

Interpretasi: Karena nilai \(\chi^2\) hitung = 24.24 lebih besar dari \(\chi^2\) tabel = 3.841 (df = 1, α = 0.05), maka \(H_0\) ditolak. Artinya, terdapat hubungan yang signifikan antara kedua variabel dalam tabel kontingensi.

1.8.3.2 Partisi Chi-Square

Misalkan kita memiliki data berikut:

Tabel 1.6 Tabel Kontingensi Gender dan Preferensi Partai Politik

Democrat Republican Independent Total
Female 495 272 590 1357
Male 330 265 498 1093
Total 825 537 1088 2450

Hipotesis:

  • Hipotesis Nol (\(H_0\)): Tidak ada hubungan antara variabel Gender dan Identifikasi Partai Politik.
  • Hipotesis Alternatif (\(H_1\)): Ada hubungan antara variabel Gender dan Identifikasi Partai Politik.

Perhitungan Manual:

Langkah 1: Lakukan Chi-Square Secara Keseluruhan

  1. Hitung Frekuensi Ekspektasi
  • Democrat, Female: \(E_{11} = \frac{1357 \times 825}{2450} = 456.95\)
  • Republican, Female: \(E_{12} = \frac{1357 \times 537}{2450} = 297.43\)
  • Independent, Female: \(E_{13} = \frac{1357 \times 1088}{2450} = 602.62\)
  • Democrat, Male: \(E_{21} = \frac{1093 \times 825}{2450} = 368.05\)
  • Republican, Male: \(E_{22} = \frac{1093 \times 537}{2450} = 239.57\)
  • Independent, Male: \(E_{23} = \frac{1093 \times 1088}{2450} = 485.38\)
  1. Hitung Statistik Uji Chi-Square \[ \chi^2 = \frac{(495 - 456.95)^2}{456.95} + \frac{(272 - 297.43)^2}{297.43} + \frac{(590 - 602.62)^2}{602.62} + \frac{(330 - 368.05)^2}{368.05} + \frac{(265 - 239.57)^2}{239.57} + \frac{(498 - 485.38)^2}{485.38} \]

\[ \chi^2 = 12.57 \]

  1. Bandingkan dengan Derajat bebas

\[ df = (I - 1)(J - 1) = (2 - 1)(3 - 1) = 2 \]

Langkah 2: Lakukan Partisi Chi-Square

Partisi 1: Democrat vs. Republican

Tabel 1.7 Partisi 1: Democrat vs. Republican

Democrat Republican
Female 495 272
Male 330 265

Chi-Square: \(\chi^2 = 11.536, \quad p < 0.001\)

Partisi 2: Democrat + Republican vs. Independent

Tabel 1.8 Partisi 2: Democrat + Republican vs. Independent

Democrat + Republican Independent
Female 767 590
Male 595 498

Chi-Square: \(\chi^2 = 1.065, \quad p = 0.698\)

Perhitungan Menggunakan R:

# Data Observasi
data_matrix <- matrix(c(495, 272, 590, 330, 265, 498), nrow = 2, byrow = TRUE)
colnames(data_matrix) <- c("Democrat", "Republican", "Independent")
rownames(data_matrix) <- c("Female", "Male")

# Uji Chi-Square
chi_test <- chisq.test(data_matrix)

# Hasil
list(Chi_Square = chi_test$statistic, P_Value = chi_test$p.value, Decision = ifelse(chi_test$p.value < 0.05, "Reject H0", "Fail to Reject H0"))
## $Chi_Square
## X-squared 
##  12.56926 
## 
## $P_Value
## [1] 0.00186475
## 
## $Decision
## [1] "Reject H0"
# Uji Chi-Square Partisi 1
chi_test1 <- chisq.test(data_matrix)

# Data Partisi 2
data_matrix2 <- matrix(c(767, 590, 595, 498), nrow = 2, byrow = TRUE)
colnames(data_matrix2) <- c("Dem+Rep", "Independent")
rownames(data_matrix2) <- c("Female", "Male")

# Uji Chi-Square Partisi 2
chi_test2 <- chisq.test(data_matrix2)

# Hasil
list(Chi_Square_Partisi1 = chi_test1, Chi_Square_Partisi2 = chi_test2)
## $Chi_Square_Partisi1
## 
##  Pearson's Chi-squared test
## 
## data:  data_matrix
## X-squared = 12.569, df = 2, p-value = 0.001865
## 
## 
## $Chi_Square_Partisi2
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data_matrix2
## X-squared = 0.98267, df = 1, p-value = 0.3215

Interpretasi:

  • Pada uji Chi-Square keseluruhan didapat \(\chi^2 = 12.57 > 5.99\), maka \(H_0\) ditolak. Terdapat hubungan yang signifikan antara Gender dan Identifikasi Partai Politik.
  • Partisi pertama menunjukkan perbedaan Democrat vs Republican signifikan (p < 0.001).
  • Partisi kedua menunjukkan (Democrat + Republican) vs Independent tidak signifikan (p = 0.698).
  • Kesimpulan: Gender berpengaruh signifikan pada pilihan antara Democrat dan Republican, tetapi tidak berpengaruh signifikan terhadap pilihan Independent.

1.8.3.3 Uji Likelihood Ratio

Diperoleh data sebagai berikut:

Tabel 1.9 Tabel Kontingensi Status Merokok dan Kanker Paru

Cancer (+) Control (-) Total
Smoker 688 650 1338
Non-Smoker 21 59 80
Total 709 709 1418

Perhitungan Manual:

Langkah 1: Hitung Frekuensi Ekspektasi

\[ \hat{\mu}_{11} = \frac{1338 \times 709}{1418} = 669 \]

\[ \hat{\mu}_{12} = \frac{1338 \times 709}{1418} = 669 \]

\[ \hat{\mu}_{21} = \frac{80 \times 709}{1418} = 40 \]

\[ \hat{\mu}_{22} = \frac{80 \times 709}{1418} = 40 \]

Langkah 2: Hitung Statistik Uji \(G^2\)

\[ G^2 = 2 \Bigg[ 688 \ln\left(\frac{688}{669}\right) + 650 \ln\left(\frac{650}{669}\right) + 21 \ln\left(\frac{21}{40}\right) + 59 \ln\left(\frac{59}{40}\right) \Bigg] \]

\[ G^2 = 2 \times (19.27 - 18.73 - 13.53 + 22.93) \]

\[ G^2 = 2 \times 9.94 = 19.88 \]

Langkah 3: Derajat bebas

\[ df = (2 - 1)(2 - 1) = 1 \]

Nilai kritis \(\chi^2\) pada \(\alpha = 0.05\) adalah 3.841.

Perhitunngan Menggunakan R:

# Data Observasi
data_matrix <- matrix(c(688, 650, 21, 59), nrow = 2, byrow = TRUE)
colnames(data_matrix) <- c("Cancer (+)", "Control (-)")
rownames(data_matrix) <- c("Smoker", "Non-Smoker")

# Hitung Frekuensi Ekspektasi
data_expected <- chisq.test(data_matrix)$expected

# Hitung Statistik G²
G2 <- 2 * sum(data_matrix * log(data_matrix / data_expected))

# Nilai kritis chi-square untuk df = 1 dan alpha = 0.05
critical_value <- qchisq(0.95, df = 1)

# Hasil
list(G2 = G2, Critical_Value = critical_value, Decision = ifelse(G2 > critical_value, "Reject H0", "Fail to Reject H0"))
## $G2
## [1] 19.87802
## 
## $Critical_Value
## [1] 3.841459
## 
## $Decision
## [1] "Reject H0"

Interpretasi: Karena \(G^2 = 9.71 > 3.841\), maka \(H_0\) ditolak, sehingga terdapat hubungan yang signifikan antara merokok dan kanker paru-paru.

1.8.3.4 Uji Exact Fisher

Misalkan diberikan contoh sebagai berikut:

  • \(N = 40 \quad (\text{total populasi})\)
  • \(K = 29 \quad (\text{jumlah bola putih dalam populasi})\)
  • \(n = 20 \quad (\text{jumlah sampel yang diambil})\)
  • \(x = 18 \quad (\text{jumlah bola putih dalam sampel})\)

Perhitungan Menggunakan R:

# Definisi parameter
N <- 40   # Total populasi
K <- 29   # Jumlah kategori sukses (bola putih)
n <- 20   # Jumlah sampel diambil
x <- 18   # Jumlah sukses dalam sampel

# Hitung probabilitas P(X = 18)
dhyper(x, m = K, n = N - K, k = n)
## [1] 0.01380413
choose(29, 18) * choose(11, 2) / choose(40, 20)
## [1] 0.01380413
choose(29, 20) * choose(11, 0) / choose(40, 20)
## [1] 7.26533e-05
choose(29, 19) * choose(11, 1) / choose(40, 20)
## [1] 0.001598373
choose(29, 18) * choose(11, 2) / choose(40, 20)
## [1] 0.01380413
choose(29, 17) * choose(11, 3) / choose(40, 20)
## [1] 0.06211857
choose(29, 16) * choose(11, 4) / choose(40, 20)
## [1] 0.162464
choose(29, 15) * choose(11, 5) / choose(40, 20)
## [1] 0.2599423
choose(29, 14) * choose(11, 6) / choose(40, 20)
## [1] 0.2599423
choose(29, 13) * choose(11, 7) / choose(40, 20)
## [1] 0.162464
choose(29, 12) * choose(11, 8) / choose(40, 20)
## [1] 0.06211857
choose(29, 11) * choose(11, 9) / choose(40, 20)
## [1] 0.01380413
choose(29, 10) * choose(11, 10) / choose(40, 20)
## [1] 0.001598373
choose(29, 9) * choose(11, 11) / choose(40, 20)
## [1] 7.26533e-05
choose(29, 18) * choose(11, 2) / choose(40, 20)
## [1] 0.01380413
p.value <- 0.00007+0.00160+0.01380+0.01380+0.00160+0.00007
p.value
## [1] 0.03094
# Menggunakan Package
data <- matrix(c(18, 2, 11, 9), nrow = 2, byrow = TRUE)
fisher.test(data)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  data
## p-value = 0.03095
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   1.147793 78.183838
## sample estimates:
## odds ratio 
##   6.994073

Interpretasi: Karena p-value = 0.03095 < 0.05, maka \(H_0\) ditolak. Artinya, terdapat hubungan yang signifikan antara kedua variabel yang diuji.

1.8.4 Analisis Residual

Misalkan terdapat tabel kontingensi sebagai berikut:

Tabel 1.10 Tabel Kontingensi Kategori A dan Kategori B pada Grup

Kategori A Kategori B Total
Grup 1 20 10 30
Grup 2 30 20 50
Total 50 30 80

Perhitungan Manual:

  1. Hitung Ekspektasi Frekuensi
  • \(E_{11} = \frac{30 \times 50}{80} = 18.75\)
  • \(E_{12} = \frac{30 \times 30}{80} = 11.25\)
  • \(E_{21} = \frac{50 \times 50}{80} = 31.25\)
  • \(E_{22} = \frac{50 \times 30}{80} = 18.75\)
  1. Hitung Pearson Residual
  • \(r_{11} = \frac{20 - 18.75}{\sqrt{18.75}} = 0.29\)
  • \(r_{12} = \frac{10 - 11.25}{\sqrt{11.25}} = -0.37\)
  • \(r_{21} = \frac{30 - 31.25}{\sqrt{31.25}} = -0.22\)
  • \(r_{22} = \frac{20 - 18.75}{\sqrt{18.75}} = 0.29\)

Perhitungan Menggunakan R:

# Data Observasi
observed <- matrix(c(20, 10, 30, 20), nrow = 2, byrow = TRUE)

# Hitung nilai ekspektasi
expected <- chisq.test(observed)$expected

# Pearson Residual
pearson_residual <- (observed - expected) / sqrt(expected)

# Standardized Residual
row_sum <- rowSums(observed)
col_sum <- colSums(observed)
total_sum <- sum(observed)

standardized_residual <- (observed - expected) / sqrt(expected * (1 - row_sum / total_sum) * (1 - col_sum / total_sum))

# Menampilkan hasil
list(
  Pearson_Residual = pearson_residual,
  Standardized_Residual = standardized_residual
)
## $Pearson_Residual
##            [,1]       [,2]
## [1,]  0.2886751 -0.3726780
## [2,] -0.2236068  0.2886751
## 
## $Standardized_Residual
##            [,1]       [,2]
## [1,]  0.5962848 -0.7698004
## [2,] -0.4618802  0.5962848

Interpretasi: Semua nilai residual < 2, sehingga tidak ada sel yang memberikan kontribusi besar terhadap nilai Chi-Square. Hal ini menunjukkan bahwa kontribusi tiap sel relatif kecil dan pola asosiasi tidak didominasi oleh satu sel tertentu.

1.8.5 Deteksi Outlier

Misalkan diberikan tabel kontingensi sebagai berikut:

Tabel 1.11 Tabel Kontingensi Sukses dan Gagal pada Dua Grup

Sukses Gagal Total
Grup A 50 20 70
Grup B 30 50 80
Total 80 70 150

Setelah menghitung residual, didapatkan:

Tabel 1.12 Pearson Residual

Pearson Residual Standardized Residual
Grup A, Sukses 2.06 2.45
Grup A, Gagal -2.20 -2.60
Grup B, Sukses -1.96 -2.10
Grup B, Gagal 2.06 2.45
  • Sel (A, Gagal) memiliki residual negatif terbesar (-2.20, -2.60) → Ini menunjukkan bahwa Grup A mengalami lebih sedikit kegagalan dari yang diprediksi oleh model independensi. Jika |r| melebihi 3, maka sel ini dapat dianggap sebagai outlier yang signifikan.
  • Sel (A, Sukses) dan (B, Gagal) memiliki residual positif cukup besar → Menunjukkan bahwa ada ketidakseimbangan dalam hubungan antara kategori.

Perhitungan Menggunakan R:

# Data Observasi
observed <- matrix(c(50, 20, 30, 50), nrow = 2, byrow = TRUE)
colnames(observed) <- c("Sukses", "Gagal")
rownames(observed) <- c("Grup A", "Grup B")

# Hitung nilai ekspektasi
expected <- chisq.test(observed)$expected

# Pearson Residual
pearson_residual <- (observed - expected) / sqrt(expected)

# Standardized Residual
row_sum <- rowSums(observed)
col_sum <- colSums(observed)
total_sum <- sum(observed)
standardized_residual <- (observed - expected) / sqrt(expected * (1 - row_sum / total_sum) * (1 - col_sum / total_sum))

# Menampilkan hasil
list(
  Observed = observed,
  Expected = expected,
  Pearson_Residual = pearson_residual,
  Standardized_Residual = standardized_residual
)
## $Observed
##        Sukses Gagal
## Grup A     50    20
## Grup B     30    50
## 
## $Expected
##          Sukses    Gagal
## Grup A 37.33333 32.66667
## Grup B 42.66667 37.33333
## 
## $Pearson_Residual
##           Sukses     Gagal
## Grup A  2.073070 -2.216205
## Grup B -1.939179  2.073070
## 
## $Standardized_Residual
##           Sukses     Gagal
## Grup A  4.155384 -4.442293
## Grup B -3.887006  4.155384

1.8.6 Tugas Mahasiswa

1.8.6.1 Studi Kasus 1

Buatlah fungsi untuk menghitung dan melakukan pegujian hipotesis untuk RD, RR, dan OR. Gunakan data berikut Dataset dari Agresti (2019, hlm. 35, Tabel 2.3):

Tabel 1.13 Tabel Kontingensi Merokok dan Kanker Paru

Smoker Lung Cancer (Cases) Control
Yes 688 650
No 21 59

Struktur tabel untuk pembuatan function

Tabel 1.14 Struktur Tabel Kontingensi Studi Kasus 1

Exposure Cases Control Total
Yes a c a+c
No b d b+d
Total a+b c+d a+b+c+d
# Uji Proporsi
## Membuat Fungsi
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))
}
## Input data
hasil <- prop_diff(a = 688, b = 21, c = 650, d = 59)

## Menampilkan hasil
print(hasil)
## $estimate
## [1] 0.2517003
## 
## $ci
## [1] 0.1516343 0.3517663
# Relative Risk
## Membuat Fungsi
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))
}
## Input data
hasil <- relative_risk(a = 688, b = 21, c = 650, d = 59)

## Menampilkan hasil
print(hasil)
## $estimate
## [1] 1.958858
## 
## $ci
## [1] 1.351735 2.838667
# Odds Ratio
## Membuat Fungsi
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))
}

## Input data
hasil <- odds_ratio(a = 688, b = 21, c = 650, d = 59)

## Menampilkan hasil
print(hasil)
## $estimate
## [1] 2.973773
## 
## $ci
## [1] 1.786737 4.949427
# 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.949474
# Perbandingan dengan Output R
library(epiR)
## Warning: package 'epiR' was built under R version 4.4.3
## Loading required package: survival
## Package epiR 2.0.91 is loaded
## Type help(epi.about) for summary information
## Type browseVignettes(package = 'epiR') to learn how to use epiR for applied epidemiological analyses
## 
table_data <- matrix(c(a, c, b, d), nrow = 2, byrow = TRUE)
colnames(table_data) <- c("Lung Cancer", "Control")
rownames(table_data) <- c("Yes", "No")
res <- epi.2by2(table_data)
print(res)
##              Outcome+    Outcome-      Total                 Inc risk *
## Exposure+         688         650       1338     51.42 (48.70 to 54.13)
## Exposure-          21          59         80     26.25 (17.04 to 37.29)
## Total             709         709       1418     50.00 (47.36 to 52.64)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 1.96 (1.35, 2.84)
## Inc odds ratio                                 2.97 (1.79, 4.95)
## Attrib risk in the exposed *                   25.17 (15.16, 35.18)
## Attrib fraction in the exposed (%)            48.95 (28.08, 65.39)
## Attrib risk in the population *                23.75 (13.76, 33.74)
## Attrib fraction in the population (%)         47.50 (29.16, 64.02)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 19.129 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

Interpretasi:

  • Risk Difference (RD) = 0.2517 (atau 25.17%) → Risiko kanker paru pada perokok lebih tinggi 25.17% secara absolut dibandingkan non-perokok.
  • Relative Risk (RR) = 1.96 → Perokok memiliki 1.96 kali lipat lebih tinggi risiko terkena kanker paru dibandingkan non-perokok.
  • Odds Ratio (OR) = 2.97 → Odds perokok terkena kanker paru 2.97 kali lebih besar dibandingkan non-perokok.

1.8.6.2 Studi Kasus 2

Seorang peneliti ingin mengetahui apakah ada hubungan antara kebiasaan merokok dan kejadian kanker paru-paru. Data yang diperoleh dari sebuah penelitian medis ditampilkan dalam tabel kontingensi berikut:

Tabel 1.15 Tabel Kontingensi Kebiasaan Merokok dan Kanker Paru-Paru

Kanker Paru (+) Kanker Paru (-) Total
Perokok 450 200 650
Bukan Perokok 50 300 350
Total 500 500 1000
  1. Hitung frekuensi ekspektasi untuk masing-masing sel berdasarkan asumsi independensi.
  2. Hitung Pearson Residual dan Standardized Residual untuk setiap sel.
  3. Interpretasikan hasil residual yang diperoleh.
  4. Tentukan apakah ada kategori yang bisa dianggap sebagai outlier dalam hubungan antara merokok dan kanker paru-paru.
# Data Observasi
observed <- matrix(c(450, 200, 50, 300), nrow = 2, byrow = TRUE)
colnames(observed) <- c("Kanker Paru (+)", "Kanker Paru (-)")
rownames(observed) <- c("Perokok", "Bukan Perokok")

# Hitung nilai ekspektasi
expected <- chisq.test(observed)$expected

# Pearson Residual
pearson_residual <- (observed - expected) / sqrt(expected)

# Standardized Residual
row_sum <- rowSums(observed)
col_sum <- colSums(observed)
total_sum <- sum(observed)
standardized_residual <- (observed - expected) / sqrt(expected * (1 - row_sum / total_sum) * (1 - col_sum / total_sum))

# Menampilkan hasil
list(
  Observed = observed,
  Expected = expected,
  Pearson_Residual = pearson_residual,
  Standardized_Residual = standardized_residual
)
## $Observed
##               Kanker Paru (+) Kanker Paru (-)
## Perokok                   450             200
## Bukan Perokok              50             300
## 
## $Expected
##               Kanker Paru (+) Kanker Paru (-)
## Perokok                   325             325
## Bukan Perokok             175             175
## 
## $Pearson_Residual
##               Kanker Paru (+) Kanker Paru (-)
## Perokok              6.933752       -6.933752
## Bukan Perokok       -9.449112        9.449112
## 
## $Standardized_Residual
##               Kanker Paru (+) Kanker Paru (-)
## Perokok              16.57484       -16.57484
## Bukan Perokok       -16.57484        16.57484

Interpretasi:

  • Perokok → Kanker Paru (+): residual positif besar → kasus lebih banyak dari yang diharapkan
  • Perokok → Kanker Paru (-): residual negatif besar → kasus lebih sedikit dari yang diharapkan
  • Bukan Perokok → Kanker Paru (+): residual negatif besar → kasus lebih sedikit dari yang diharapkan
  • Bukan Perokok → Kanker Paru (-): residual positif besar → kasus lebih banyak dari yang diharapkan
  • Kesimpulan: terdapat hubungan yang kuat antara status merokok dan kejadian kanker paru-paru.

2 Tugas 2: Inferensi Tabel Kontingensi Dua Arah

2.1 Pendahuluan

Inferensi tabel kontingensi dua arah merupakan salah satu metode dalam statistika yang digunakan untuk menganalisis hubungan antara dua variabel kategorik. Melalui penyajian data dalam bentuk tabel, metode ini memungkinkan peneliti untuk mengidentifikasi ada atau tidaknya keterkaitan antara kedua variabel, serta mengukur kekuatan asosiasi yang terjadi. Analisis ini sering digunakan dalam berbagai bidang untuk memahami pola hubungan dalam data kategorik secara lebih sistematis.

Load package yang diperlukan:

# Load paket yang diperlukan
library(epitools)     # untuk RR, OR, dan uji terkait
## 
## Attaching package: 'epitools'
## The following object is masked from 'package:survival':
## 
##     ratetable
library(vcd)          # untuk mosaic plot dan ukuran asosiasi
## Warning: package 'vcd' was built under R version 4.4.3
## Loading required package: grid
## 
## Attaching package: 'vcd'
## The following object is masked from 'package:epitools':
## 
##     oddsratio
library(ggplot2)      # untuk visualisasi
## Warning: package 'ggplot2' was built under R version 4.4.3
library(knitr)        # untuk tabel rapi
library(kableExtra)   # untuk styling tabel
## Warning: package 'kableExtra' was built under R version 4.4.3
library(MASS)         # untuk likelihood ratio
## Warning: package 'MASS' was built under R version 4.4.3
library(tidyverse)    # untuk manipulasi data
## Warning: package 'tidyverse' was built under R version 4.4.3
## Warning: package 'tidyr' was built under R version 4.4.3
## Warning: package 'readr' was built under R version 4.4.3
## Warning: package 'purrr' was built under R version 4.4.3
## Warning: package 'dplyr' was built under R version 4.4.3
## Warning: package 'forcats' was built under R version 4.4.3
## Warning: package 'lubridate' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.1.6
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.4     ✔ tidyr     1.3.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::group_rows() masks kableExtra::group_rows()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ dplyr::select()     masks MASS::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(scales)       # untuk format angka
## Warning: package 'scales' was built under R version 4.4.3
## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
library(corrplot)     # untuk visualisasi residual
## Warning: package 'corrplot' was built under R version 4.4.3
## corrplot 0.95 loaded

2.2 Kasus 1: Tabel Kontingensi 2×2 (Merokok dan Kanker Paru)

2.2.1 Penyusunan Tabel Kontingensi

Berikut merupakan data hubungan antara kebiasaan merokok (smoking status) dan kejadian kanker paru (lung cancer):

# Data
tabel1 <- matrix(
  c(688, 650, 21, 59),
  nrow = 2,
  byrow = TRUE
)

# Hitung total
row_total <- rowSums(tabel1)
col_total <- colSums(tabel1)
grand_total <- sum(tabel1)

# Tampilkan sebagai HTML
cat('
<p style="text-align:center; font-weight:bold;">
Tabel 2.1 Tabel Kontingensi 2x2 Kebiasaan Merokok dan Kejadian Kanker Paru
</p>
<table class="my-table">
  <tr>
    <th></th>
    <th>Cancer (+)</th>
    <th>Control (-)</th>
    <th>Total</th>
  </tr>
  <tr>
    <td><b>Smoker</b></td>
    <td>', tabel1[1,1], '</td>
    <td>', tabel1[1,2], '</td>
    <td><b>', row_total[1], '</b></td>
  </tr>
  <tr>
    <td><b>Non-Smoker</b></td>
    <td>', tabel1[2,1], '</td>
    <td>', tabel1[2,2], '</td>
    <td><b>', row_total[2], '</b></td>
  </tr>
  <tr>
    <td><b>Total</b></td>
    <td><b>', col_total[1], '</b></td>
    <td><b>', col_total[2], '</b></td>
    <td><b>', grand_total, '</b></td>
  </tr>
</table>
')

Tabel 2.1 Tabel Kontingensi 2x2 Kebiasaan Merokok dan Kejadian Kanker Paru

Cancer (+) Control (-) Total
Smoker 688 650 1338
Non-Smoker 21 59 80
Total 709 709 1418

Keterangan notasi:

  • \(a = 688\) (Smoker, Cancer+)
  • \(b = 650\) (Smoker, Control-)
  • \(c = 21\) (Non-Smoker, Cancer+)
  • \(d = 59\) (Non-Smoker, Control-)
  • \(n_1 = 1338\) (total Smoker)
  • \(n_2 = 80\) (total Non-Smoker)
  • \(N = 1418\)

2.2.2 Estimasi

2.2.2.1 Estimasi Titik Proporsi

Estimasi titik proporsi digunakan untuk menggambarkan peluang kejadian pada masing-masing kelompok.

Perhitungan Manual: Secara umum, estimasi proporsi sampel untuk masing-masing kelompok dinyatakan sebagai:

\[ \hat{p}_1 = \frac{x_1}{n_1}, \qquad \hat{p}_2 = \frac{x_2}{n_2} \]

  • Estimasi titik proporsi pada kelompok Smoker:

\[\hat{p}_1 = \frac{a}{n_1} = \frac{688}{1338}\]

  • Estimasi titik proporsi pada kelompok Non-Smoker:

\[\hat{p}_2 = \frac{c}{n_2} = \frac{21}{80}\] Perhitungan Menggunakan R:

a <- 688; b <- 650; c_val <- 21; d <- 59
n1 <- a + b   # total Smoker
n2 <- c_val + d  # total Non-Smoker
N  <- n1 + n2

p1_hat <- a / n1
p2_hat <- c_val / n2

cat("=== Estimasi Titik Proporsi ===\n")
## === Estimasi Titik Proporsi ===
cat(sprintf("Proporsi kanker paru pada Smoker     : p1 = %d/%d = %.4f (%.2f%%)\n",
            a, n1, p1_hat, p1_hat * 100))
## Proporsi kanker paru pada Smoker     : p1 = 688/1338 = 0.5142 (51.42%)
cat(sprintf("Proporsi kanker paru pada Non-Smoker : p2 = %d/%d = %.4f (%.2f%%)\n",
            c_val, n2, p2_hat, p2_hat * 100))
## Proporsi kanker paru pada Non-Smoker : p2 = 21/80 = 0.2625 (26.25%)

2.2.2.2 Estimasi Interval

Untuk memperoleh estimasi yang lebih akurat, digunakan interval kepercayaan terhadap proporsi masing-masing kelompok serta ukuran asosiasi yang meliputi Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR).

2.2.2.2.1 Proporsi Masing-Masing Kelompok

Salah satu metode yang direkomendasikan adalah metode Wilson, karena lebih stabil dibandingkan metode Wald. Secara umum, interval Wilson dirumuskan sebagai berikut:

\[ CI_{95\%}(\hat{p}) = \frac{\hat{p} + \frac{z^2}{2n} \pm z\sqrt{\frac{\hat{p}(1-\hat{p})}{n} + \frac{z^2}{4n^2}}}{1 + \frac{z^2}{n}} \]

Interval kepercayaan proporsi dihitung secara manual dan menggunakan fungsi prop.test() di R, dengan hasil sebagai berikut:

# Perhitungan Menggunakan Fungsi `prop.test()`
ci_p1_prop <- prop.test(a, n1, conf.level = 0.95)$conf.int
ci_p2_prop <- prop.test(c_val, n2, conf.level = 0.95)$conf.int

ci_p1_prop

[1] 0.4870445 0.5412736 attr(,“conf.level”) [1] 0.95

ci_p2_prop

[1] 0.1733064 0.3748263 attr(,“conf.level”) [1] 0.95

# Perhitungan Manual
z <- qnorm(0.975)
RD <- p1_hat - p2_hat
SE_RD <- sqrt(p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2)
CI_RD_lower <- RD - z * SE_RD
CI_RD_upper <- RD + z * SE_RD

## Fungsi Wilson CI
wilson_ci <- function(x, n, conf = 0.95) {
  z_val <- qnorm(1 - (1 - conf) / 2)
  p_hat <- x / n
  center <- (p_hat + z_val^2 / (2 * n)) / (1 + z_val^2 / n)
  margin <- z_val * sqrt(p_hat * (1 - p_hat) / n + z_val^2 / (4 * n^2)) / (1 + z_val^2 / n)
  c(lower = center - margin, upper = center + margin)
}

ci_p1 <- wilson_ci(a, n1)
ci_p2 <- wilson_ci(c_val, n2)

bt1 <- binom.test(a, n1)
bt2 <- binom.test(c_val, n2)
  
# HTML table
cat('
<p style="text-align:center; font-weight:bold;">
Tabel 2.2 Estimasi Proporsi dan Interval Kepercayaan 95%
</p>

<table class="my-table">
  <tr>
    <th>Kelompok</th>
    <th>n</th>
    <th>Kejadian</th>
    <th>p̂</th>
    <th>CI Lower (Wilson)</th>
    <th>CI Upper (Wilson)</th>
    <th>CI Lower (Exact)</th>
    <th>CI Upper (Exact)</th>
  </tr>
  <tr>
    <td><b>Smoker</b></td>
    <td>', n1, '</td>
    <td>', a, '</td>
    <td>', round(p1_hat,4), '</td>
    <td>', round(ci_p1["lower"],4), '</td>
    <td>', round(ci_p1["upper"],4), '</td>
    <td>', round(bt1$conf.int[1],4), '</td>
    <td>', round(bt1$conf.int[2],4), '</td>
  </tr>
  <tr>
    <td><b>Non-Smoker</b></td>
    <td>', n2, '</td>
    <td>', c_val, '</td>
    <td>', round(p2_hat,4), '</td>
    <td>', round(ci_p2["lower"],4), '</td>
    <td>', round(ci_p2["upper"],4), '</td>
    <td>', round(bt2$conf.int[1],4), '</td>
    <td>', round(bt2$conf.int[2],4), '</td>
  </tr>
</table>
')

Tabel 2.2 Estimasi Proporsi dan Interval Kepercayaan 95%

Kelompok n Kejadian CI Lower (Wilson) CI Upper (Wilson) CI Lower (Exact) CI Upper (Exact)
Smoker 1338 688 0.5142 0.4874 0.5409 0.487 0.5413
Non-Smoker 80 21 0.2625 0.1786 0.3682 0.1704 0.3729

Interpretasi: Dengan tingkat kepercayaan 95%, dapat dinyatakan bahwa proporsi sebenarnya pada populasi kelompok perokok berada dalam interval (0.4874, 0.5409), sedangkan pada kelompok non-perokok berada dalam interval (0.1786, 0.3682).

2.2.2.2.2 Risk Difference (RD)

Interval kepercayaan 95% untuk Risk Difference dihitung sebagai berikut:

\[ RD = \hat{p}_1 - \hat{p}_2, \quad CI_{95\%}(RD) = RD \pm z \cdot SE(RD) \]

RD <- p1_hat - p2_hat
SE_RD <- sqrt(p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2)
CI_RD_lower <- RD - z * SE_RD
CI_RD_upper <- RD + z * SE_RD

cat("=== Risk Difference (RD) ===\n")
## === Risk Difference (RD) ===
cat(sprintf("RD = %.4f - %.4f = %.4f\n", p1_hat, p2_hat, RD))
## RD = 0.5142 - 0.2625 = 0.2517
cat(sprintf("SE(RD) = %.4f\n", SE_RD))
## SE(RD) = 0.0511
cat(sprintf("95%% CI RD : (%.4f, %.4f)\n", CI_RD_lower, CI_RD_upper))
## 95% CI RD : (0.1516, 0.3518)

Interpretasi: Dengan tingkat kepercayaan 95%, dapat dinyatakan bahwa perbedaan proporsi antara kelompok perokok dan non-perokok berada dalam interval (0.1516, 0.3518).

2.2.2.2.3 Relative Risk (RR)

Interval kepercayaan 95% untuk Relative Risk dihitung sebagai berikut:

\[RR = \frac{\hat{p}_1}{\hat{p}_2}, \quad SE(\ln RR) = \sqrt{\frac{1-\hat{p}_1}{n_1 \hat{p}_1} + \frac{1-\hat{p}_2}{n_2 \hat{p}_2}}\]

RR <- p1_hat / p2_hat
SE_lnRR <- sqrt((1 - p1_hat) / (n1 * p1_hat) + (1 - p2_hat) / (n2 * p2_hat))
CI_RR_lower <- exp(log(RR) - z * SE_lnRR)
CI_RR_upper <- exp(log(RR) + z * SE_lnRR)

cat("=== Relative Risk (RR) ===\n")
## === Relative Risk (RR) ===
cat(sprintf("RR = %.4f / %.4f = %.4f\n", p1_hat, p2_hat, RR))
## RR = 0.5142 / 0.2625 = 1.9589
cat(sprintf("SE(ln RR) = %.4f\n", SE_lnRR))
## SE(ln RR) = 0.1893
cat(sprintf("95%% CI RR : (%.4f, %.4f)\n", CI_RR_lower, CI_RR_upper))
## 95% CI RR : (1.3517, 2.8387)

Interpretasi: Berdasarkan hasil perhitungan dengan tingkat kepercayaan 95%, diperoleh nilai relative risk (RR) sebesar 1.9589 dengan interval kepercayaan (1.3517 , 2.8387). Hal ini menunjukkan bahwa risiko kejadian pada kelompok perokok diperkirakan sebesar 1.9589 kali lebih besar dibandingkan dengan kelompok non-perokok.

2.2.2.2.4 Odds Ratio (OR)

Interval kepercayaan 95% untuk Odds Ratio dihitung sebagai berikut:

\[OR = \frac{a \cdot d}{b \cdot c}, \quad SE(\ln OR) = \sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}}\]

# Data
a <- 688
b <- 650
c_val <- 21
d <- 59

n1 <- a + b
n2 <- c_val + d

# Proporsi
p1_hat <- a / n1
p2_hat <- c_val / n2

# Z
z <- qnorm(0.975)

OR <- (a * d) / (b * c_val)
SE_lnOR <- sqrt(1/a + 1/b + 1/c_val + 1/d)
CI_OR_lower <- exp(log(OR) - z * SE_lnOR)
CI_OR_upper <- exp(log(OR) + z * SE_lnOR)

cat("=== Odds Ratio (OR) ===\n")
## === Odds Ratio (OR) ===
cat(sprintf("RR = %.4f / %.4f = %.4f\n", p1_hat, p2_hat, OR))
## RR = 0.5142 / 0.2625 = 2.9738
cat(sprintf("SE(ln OR) = %.4f\n", SE_lnOR))
## SE(ln OR) = 0.2599
cat(sprintf("95%% CI RR : (%.4f, %.4f)\n", CI_OR_lower, CI_OR_upper))
## 95% CI RR : (1.7867, 4.9494)

Interpretasi: Dengan tingkat kepercayaan 95%, diperoleh interval kepercayaan untuk Odds Ratio antara kelompok perokok dan non-perokok berada dalam interval (1.7867, 4.9494), dengan nilai OR sebesar 2.9738. Hal ini menunjukkan bahwa peluang terjadinya kejadian pada kelompok perokok adalah sebesar 2.9738 kali lebih besar dibandingkan dengan kelompok non-perokok.

# ===== Proporsi (Wilson test) =====
ci_p1 <- prop.test(a, n1, conf.level = 0.95)$conf.int
ci_p2 <- prop.test(c_val, n2, conf.level = 0.95)$conf.int

# ===== RD =====
RD <- p1_hat - p2_hat
SE_RD <- sqrt(p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2)
CI_RD_lower <- RD - z * SE_RD
CI_RD_upper <- RD + z * SE_RD

# ===== RR =====
RR <- p1_hat / p2_hat
SE_lnRR <- sqrt((1 - p1_hat)/(a) + (1 - p2_hat)/(c_val))
CI_RR_lower <- exp(log(RR) - z * SE_lnRR)
CI_RR_upper <- exp(log(RR) + z * SE_lnRR)

# ===== OR =====
OR <- (a * d) / (b * c_val)
SE_lnOR <- sqrt(1/a + 1/b + 1/c_val + 1/d)
CI_OR_lower <- exp(log(OR) - z * SE_lnOR)
CI_OR_upper <- exp(log(OR) + z * SE_lnOR)

OR <- (a * d) / (b * c_val)
SE_lnOR <- sqrt(1/a + 1/b + 1/c_val + 1/d)
CI_OR_lower <- exp(log(OR) - z * SE_lnOR)
CI_OR_upper <- exp(log(OR) + z * SE_lnOR)

cat('
<p style="text-align:center; font-weight:bold; margin-bottom:5px;">
Tabel 2.3 Ringkasan Estimasi dan Interval Kepercayaan 95%
</p>

<table class="my-table">
  <tr>
    <th>Parameter</th>
    <th>Estimasi</th>
    <th>CI 95% Lower</th>
    <th>CI 95% Upper</th>
    <th>Keterangan</th>
  </tr>
  <tr>
    <td><b>Proporsi (Perokok)</b></td>
    <td>', round(p1_hat,4), '</td>
    <td>', round(ci_p1[1],4), '</td>
    <td>', round(ci_p1[2],4), '</td>
    <td>Proporsi kejadian pada kelompok perokok</td>
  </tr>

  <tr>
    <td><b>Proporsi (Non-Perokok)</b></td>
    <td>', round(p2_hat,4), '</td>
    <td>', round(ci_p2[1],4), '</td>
    <td>', round(ci_p2[2],4), '</td>
    <td>Proporsi kejadian pada kelompok non-perokok</td>
  </tr>
  <tr>
    <td><b>Risk Difference (RD)</b></td>
    <td>', round(RD,4), '</td>
    <td>', round(CI_RD_lower,4), '</td>
    <td>', round(CI_RD_upper,4), '</td>
    <td>Perbedaan risiko absolut</td>
  </tr>
  <tr>
    <td><b>Relative Risk (RR)</b></td>
    <td>', round(RR,4), '</td>
    <td>', round(CI_RR_lower,4), '</td>
    <td>', round(CI_RR_upper,4), '</td>
    <td>Rasio risiko relatif</td>
  </tr>
  <tr>
    <td><b>Odds Ratio (OR)</b></td>
    <td>', round(OR,4), '</td>
    <td>', round(CI_OR_lower,4), '</td>
    <td>', round(CI_OR_upper,4), '</td>
    <td>Rasio odds (cocok studi kasus-kontrol)</td>
  </tr>
</table>
')

Tabel 2.3 Ringkasan Estimasi dan Interval Kepercayaan 95%

Parameter Estimasi CI 95% Lower CI 95% Upper Keterangan
Proporsi (Perokok) 0.5142 0.487 0.5413 Proporsi kejadian pada kelompok perokok
Proporsi (Non-Perokok) 0.2625 0.1733 0.3748 Proporsi kejadian pada kelompok non-perokok
Risk Difference (RD) 0.2517 0.1516 0.3518 Perbedaan risiko absolut
Relative Risk (RR) 1.9589 1.3517 2.8387 Rasio risiko relatif
Odds Ratio (OR) 2.9738 1.7867 4.9494 Rasio odds (cocok studi kasus-kontrol)

2.2.3 Uji Hipotesis

2.2.3.1 Uji Proporsi

Hipotesis:

  • Hipotesis Nol \(H_0\): Tidak ada perbedaan proporsi antara kelompok perokok dan non-perokok \((p_1 = p_2)\)
  • Hipotesis Alternatif \(H_1\): Terdapat perbedaan proporsi antara kelompok perokok dan non-perokok \((p_1 \neq p_2)\)

Statistik uji:

\[Z = \frac{\hat{p}_1 - \hat{p}_2}{\sqrt{\hat{p}(1-\hat{p})\left(\frac{1}{n_1} + \frac{1}{n_2}\right)}}\]

di mana \(\hat{p} = \frac{a+c}{N}\) adalah proporsi gabungan.

p_pool <- (a + c_val) / N
SE_pool <- sqrt(p_pool * (1 - p_pool) * (1/n1 + 1/n2))
Z_stat <- (p1_hat - p2_hat) / SE_pool
p_value_Z <- 2 * pnorm(-abs(Z_stat))

cat("=== Uji Dua Proporsi (Z-test) ===\n")
## === Uji Dua Proporsi (Z-test) ===
cat(sprintf("Proporsi gabungan p̄ = (%d + %d) / %d = %.4f\n", a, c_val, N, p_pool))
## Proporsi gabungan p̄ = (688 + 21) / 1418 = 0.5000
cat(sprintf("SE gabungan         = %.4f\n", SE_pool))
## SE gabungan         = 0.0575
cat(sprintf("Statistik Z         = %.4f\n", Z_stat))
## Statistik Z         = 4.3737
cat(sprintf("p-value (2-sisi)    = %.6f\n", p_value_Z))
## p-value (2-sisi)    = 0.000012
cat(sprintf("Keputusan           : %s H0 pada α = 0.05\n",
            ifelse(p_value_Z < 0.05, "Tolak", "Gagal Tolak")))
## Keputusan           : Tolak H0 pada α = 0.05
# Verifikasi dengan prop.test
prop_test <- prop.test(c(a, c_val), c(n1, n2), correct = FALSE)
cat("\n--- Verifikasi dengan prop.test() ---\n")
## 
## --- Verifikasi dengan prop.test() ---
print(prop_test)
## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  c(a, c_val) out of c(n1, n2)
## 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.2625000

Interpretasi: Karena p-value < 0,05, maka keputusan yang diambil adalah menolak hipotesis nol (\(H_0\)). Artinya, terdapat perbedaan proporsi yang signifikan antara kedua kelompok.

2.2.3.2 Uji Chi-Square

Hipotesis:

  • Hipotesis Nol \(H_0\): Merokok dan Kanker Paru independen
  • Hipotesis Alternatif \(H_1\): Terdapat hubungan antara variabel Merokok dan Kanker Paru

Statistik uji:

\[\chi^2 = \sum_{i,j} \frac{(O_{ij} - E_{ij})^2}{E_{ij}}, \quad E_{ij} = \frac{n_i \cdot n_j}{N}\]

chi_test <- chisq.test(tabel1, correct = FALSE)

cat("=== Uji Chi-Square Independensi ===\n")
## === Uji Chi-Square Independensi ===
cat("\nFrekuensi Observasi:\n")
## 
## Frekuensi Observasi:
print(tabel1)
##      [,1] [,2]
## [1,]  688  650
## [2,]   21   59
cat("\nFrekuensi Harapan (Expected):\n")
## 
## Frekuensi Harapan (Expected):
print(round(chi_test$expected, 2))
##      [,1] [,2]
## [1,]  669  669
## [2,]   40   40
cat(sprintf("\nStatistik χ² = %.4f\n", chi_test$statistic))
## 
## Statistik χ² = 19.1292
cat(sprintf("df           = %d\n", chi_test$parameter))
## df           = 1
cat(sprintf("p-value      = %.6f\n", chi_test$p.value))
## p-value      = 0.000012
cat(sprintf("Keputusan    : %s H0 pada α = 0.05\n",
            ifelse(chi_test$p.value < 0.05, "Tolak", "Gagal Tolak")))
## Keputusan    : Tolak H0 pada α = 0.05

Interpretasi: Karena p-value < 0,05, maka keputusan yang diambil adalah menolak hipotesis nol (\(H_0\)). Artinya, terdapat hubungan yang signifikan antara kedua variabel.

2.2.3.3 Uji Likelihood Ratio (G²)

Uji Likelihood Ratio merupakan alternatif dari uji Chi-Square untuk menguji asosiasi antar variabel kategorik.

Hipotesis:

  • Hipotesis Nol \(H_0\): Merokok dan Kanker Paru independen
  • Hipotesis Alternatif \(H_1\): Terdapat hubungan antara variabel Merokok dan Kanker Paru

Statistik uji:

\[G^2 = 2 \sum_{i,j} O_{ij} \ln\left(\frac{O_{ij}}{E_{ij}}\right)\]

O <- as.vector(tabel1)
E <- as.vector(chi_test$expected)
G2 <- 2 * sum(O * log(O / E))
df_G2 <- (nrow(tabel1) - 1) * (ncol(tabel1) - 1)
p_value_G2 <- pchisq(G2, df = df_G2, lower.tail = FALSE)

cat("=== Uji Likelihood Ratio (G²) ===\n")
## === Uji Likelihood Ratio (G²) ===
cat("\nRincian perhitungan G²:\n")
## 
## Rincian perhitungan G²:
lr_detail <- data.frame(
  Sel = c("Smoker-Cancer+", "Smoker-Control-", "NonSmoker-Cancer+", "NonSmoker-Control-"),
  O = O,
  E = round(E, 4),
  O_lnO_E = round(O * log(O / E), 4)
)
print(lr_detail)
##                  Sel   O   E  O_lnO_E
## 1     Smoker-Cancer+ 688 669  19.2673
## 2    Smoker-Control-  21  40 -13.5315
## 3  NonSmoker-Cancer+ 650 669 -18.7276
## 4 NonSmoker-Control-  59  40  22.9308
cat(sprintf("\nG² = 2 × Σ(O × ln(O/E)) = 2 × %.4f = %.4f\n", sum(O * log(O / E)), G2))
## 
## G² = 2 × Σ(O × ln(O/E)) = 2 × 9.9390 = 19.8780
cat(sprintf("df     = %d\n", df_G2))
## df     = 1
cat(sprintf("p-value = %.6f\n", p_value_G2))
## p-value = 0.000008
cat(sprintf("Keputusan : %s H0 pada α = 0.05\n",
            ifelse(p_value_G2 < 0.05, "Tolak", "Gagal Tolak")))
## Keputusan : Tolak H0 pada α = 0.05

Interpretasi: Karena p-value < 0,05, maka keputusan yang diambil adalah menolak hipotesis nol (\(H_0\)). Artinya, terdapat hubungan yang signifikan antara kedua variabel.

2.2.3.4 Fisher Exact Test

Fisher exact test digunakan ketika asumsi chi-square tidak terpenuhi (frekuensi harapan < 5), dengan menghitung probabilitas eksak berdasarkan distribusi hipergeometrik.

Hipotesis:

  • Hipotesis Nol \(H_0\): Merokok dan Kanker Paru independen
  • Hipotesis Alternatif \(H_1\): Terdapat hubungan antara variabel Merokok dan Kanker Paru

Statistik Uji:

\[p = \frac{\binom{n_1}{a}\binom{n_2}{c}}{\binom{N}{a+c}} = \frac{n_1! \cdot n_2! \cdot (a+c)! \cdot (b+d)!}{N! \cdot a! \cdot b! \cdot c! \cdot d!}\]

fisher_test <- fisher.test(tabel1)

cat("=== Fisher Exact Test ===\n")
## === Fisher Exact Test ===
cat(sprintf("Odds Ratio (MLE) = %.4f\n", fisher_test$estimate))
## Odds Ratio (MLE) = 2.9716
cat(sprintf("95%% CI OR        = (%.4f, %.4f)\n",
            fisher_test$conf.int[1], fisher_test$conf.int[2]))
## 95% CI OR        = (1.7556, 5.2107)
cat(sprintf("p-value          = %.6f\n", fisher_test$p.value))
## p-value          = 0.000015
cat(sprintf("Keputusan        : %s H0 pada α = 0.05\n",
            ifelse(fisher_test$p.value < 0.05, "Tolak", "Gagal Tolak")))
## Keputusan        : Tolak H0 pada α = 0.05

Interpretasi: Karena p-value < 0,05, maka keputusan yang diambil adalah menolak hipotesis nol (\(H_0\)). Artinya, terdapat hubungan yang signifikan antara kedua variabel.

2.2.3.5 Perbandingan Hasil Uji

comparison <- data.frame(
  Uji = c("Uji Dua Proporsi (Z)", "Chi-Square (χ²)", "Likelihood Ratio (G²)", "Fisher Exact Test"),
  Hipotesis_Nol = rep("p₁ = p₂ (independensi)", 4),
  Statistik_Uji = c(
    paste0("Z = ", round(Z_stat, 4)),
    paste0("χ² = ", round(chi_test$statistic, 4)),
    paste0("G² = ", round(G2, 4)),
    "— (probabilitas eksak)"
  ),
  df = c("—", "1", "1", "—"),
  p_value = c(
    formatC(p_value_Z, format = "f", digits = 6),
    formatC(chi_test$p.value, format = "f", digits = 6),
    formatC(p_value_G2, format = "f", digits = 6),
    formatC(fisher_test$p.value, format = "f", digits = 6)
  ),
  Keputusan = rep("Tolak H₀", 4),
  Interpretasi = c(
    "Ada perbedaan proporsi signifikan",
    "Ada asosiasi signifikan (asimtotik)",
    "Ada asosiasi signifikan (asimtotik)",
    "Ada asosiasi signifikan (eksak)"
  )
)

cat('
<p style="text-align:center; font-weight:bold; margin-bottom:5px;">
Tabel 2.4 Perbandingan Hasil Keempat Metode Pengujian
</p>

<table class="my-table">
  <tr>
    <th>Metode Uji</th>
    <th>Hipotesis Nol</th>
    <th>Statistik Uji</th>
    <th>df</th>
    <th>p-value</th>
    <th>Keputusan</th>
    <th>Interpretasi</th>
  </tr>
  <tr>
    <td><b>Uji Dua Proporsi (Z)</b></td>
    <td>p₁ = p₂ (independensi)</td>
    <td>Z = ', round(Z_stat,4), '</td>
    <td>—</td>
    <td>', formatC(p_value_Z, format="f", digits=6), '</td>
    <td>Tolak H₀</td>
    <td>Ada perbedaan proporsi signifikan</td>
  </tr>
  <tr>
    <td><b>Chi-Square (χ²)</b></td>
    <td>p₁ = p₂ (independensi)</td>
    <td>χ² = ', round(chi_test$statistic,4), '</td>
    <td>1</td>
    <td>', formatC(chi_test$p.value, format="f", digits=6), '</td>
    <td>Tolak H₀</td>
    <td>Ada asosiasi signifikan (asimtotik)</td>
  </tr>
  <tr>
    <td><b>Likelihood Ratio (G²)</b></td>
    <td>p₁ = p₂ (independensi)</td>
    <td>G² = ', round(G2,4), '</td>
    <td>1</td>
    <td>', formatC(p_value_G2, format="f", digits=6), '</td>
    <td>Tolak H₀</td>
    <td>Ada asosiasi signifikan (asimtotik)</td>
  </tr>
  <tr>
    <td><b>Fisher Exact Test</b></td>
    <td>p₁ = p₂ (independensi)</td>
    <td>—</td>
    <td>—</td>
    <td>', formatC(fisher_test$p.value, format="f", digits=6), '</td>
    <td>Tolak H₀</td>
    <td>Ada asosiasi signifikan (eksak)</td>
  </tr>
</table>
')

Tabel 2.4 Perbandingan Hasil Keempat Metode Pengujian

Metode Uji Hipotesis Nol Statistik Uji df p-value Keputusan Interpretasi
Uji Dua Proporsi (Z) p₁ = p₂ (independensi) Z = 4.3737 0.000012 Tolak H₀ Ada perbedaan proporsi signifikan
Chi-Square (χ²) p₁ = p₂ (independensi) χ² = 19.1292 1 0.000012 Tolak H₀ Ada asosiasi signifikan (asimtotik)
Likelihood Ratio (G²) p₁ = p₂ (independensi) G² = 19.878 1 0.000008 Tolak H₀ Ada asosiasi signifikan (asimtotik)
Fisher Exact Test p₁ = p₂ (independensi) 0.000015 Tolak H₀ Ada asosiasi signifikan (eksak)

Interpretasi:

  • Uji Z dan Chi-Square ekuivalen secara matematis (\(Z^2 = \chi^2\)); keduanya bersifat asimtotik.
  • Likelihood Ratio (G²) juga asimtotik tetapi berbasis informasi; umumnya lebih akurat untuk sampel kecil.
  • Fisher Exact Test memberikan probabilitas eksak dari distribusi hipergeometrik; tidak bergantung pada asumsi asimtotik.
  • Keempat uji menghasilkan keputusan yang konsisten: tolak \(H_0\).

2.2.4 Visualisasi

par(mfrow = c(1, 2))

# Mosaic plot
mosaicplot(tabel1,
           main = "Mosaic Plot: Merokok vs Kanker Paru",
           color = c("#e74c3c", "#3498db"),
           xlab = "Status Merokok",
           ylab = "Status Kanker",
           cex.axis = 1.1)

# Bar plot proporsi
prop_data <- data.frame(
  Kelompok = c("Smoker", "Non-Smoker"),
  Proporsi = c(p1_hat, p2_hat),
  Lower = c(ci_p1["lower"], ci_p2["lower"]),
  Upper = c(ci_p1["upper"], ci_p2["upper"])
)

barplot(
  height = c(p1_hat, p2_hat),
  names.arg = c("Smoker", "Non-Smoker"),
  col = c("#e74c3c", "#3498db"),
  main = "Proporsi Kanker Paru per Kelompok",
  ylab = "Proporsi",
  ylim = c(0, 0.8),
  border = NA
)
arrows(x0 = c(0.7, 1.9), y0 = c(ci_p1["lower"], ci_p2["lower"]),
       x1 = c(0.7, 1.9), y1 = c(ci_p1["upper"], ci_p2["upper"]),
       angle = 90, code = 3, length = 0.1, lwd = 2)
abline(h = p1_hat, lty = 2, col = "gray")

par(mfrow = c(1, 1))

2.2.5 Kesimpulan

Berdasarkan seluruh analisis pada Kasus 1, diperoleh kesimpulan sebagai berikut:

  1. Estimasi proporsi: Proporsi kanker paru pada kelompok Smoker (\(\hat{p}_1 = 0.5142\)) jauh lebih tinggi dibandingkan Non-Smoker (\(\hat{p}_2 = 0.2625\)).

  2. Ukuran asosiasi:

    • \(RD = 0.2517\): Perokok memiliki risiko kanker paru 25.2% lebih tinggi secara absolut dibanding non-perokok.
    • \(RR = 1.9589\): Risiko kanker paru perokok 1.96 kali lipat dibandingkan non-perokok.
    • \(OR = 2.9738\): Odds kejadian kanker paru perokok 2.97 kali lebih tinggi dibandingkan non-perokok.
  3. Pengujian hipotesis: Keempat metode uji (dua proporsi, chi-square, likelihood ratio, Fisher exact) secara konsisten menolak \(H_0\) dengan p-value yang sangat kecil (< 0.0001).

  4. Kesimpulan substantif: Terdapat hubungan yang sangat signifikan dan kuat antara kebiasaan merokok dan kejadian kanker paru. Bukti statistik ini konsisten dengan berbagai penelitian epidemiologi yang telah mapan.


2.3 Kasus 2: Tabel Kontingensi 2×3 (Gender dan Identifikasi Partai Politik)

2.3.1 Penyusunan Tabel Kontingensi

Berikut merupakan data gender dan identifikasi partai politik:

tabel2 <- matrix(
  c(495, 272, 590,
    330, 265, 498),
  nrow = 2,
  byrow = TRUE,
  dimnames = list(
    c("Female", "Male"),
    c("Democrat", "Republican", "Independent")
  )
)

# Hitung total
row_total <- rowSums(tabel2)
col_total <- colSums(tabel2)
grand_total <- sum(tabel2)

cat('
<p style="text-align:center; font-weight:bold; margin-bottom:5px;">
Tabel 2.5 Kontingensi 2×3: Gender dan Identifikasi Partai Politik
</p>

<table class="my-table">
  <tr>
    <th>Gender</th>
    <th>Democrat</th>
    <th>Republican</th>
    <th>Independent</th>
    <th>Total</th>
  </tr>
  <tr>
    <td><b>Female</b></td>
    <td>', tabel2[1,1], '</td>
    <td>', tabel2[1,2], '</td>
    <td>', tabel2[1,3], '</td>
    <td><b>', row_total[1], '</b></td>
  </tr>
  <tr>
    <td><b>Male</b></td>
    <td>', tabel2[2,1], '</td>
    <td>', tabel2[2,2], '</td>
    <td>', tabel2[2,3], '</td>
    <td><b>', row_total[2], '</b></td>
  </tr>
  <tr>
    <td><b>Total</b></td>
    <td><b>', col_total[1], '</b></td>
    <td><b>', col_total[2], '</b></td>
    <td><b>', col_total[3], '</b></td>
    <td><b>', grand_total, '</b></td>
  </tr>
</table>
')

Tabel 2.5 Kontingensi 2×3: Gender dan Identifikasi Partai Politik

Gender Democrat Republican Independent Total
Female 495 272 590 1357
Male 330 265 498 1093
Total 825 537 1088 2450

2.3.2 Frekuensi Harapan

Frekuensi harapan dihitung dengan rumus:

\[E_{ij} = \frac{n_{i\cdot} \times n_{\cdot j}}{N}\]

chi2_full <- chisq.test(tabel2, correct = FALSE)

E_mat <- chi2_full$expected

# Hitung total
row_total <- rowSums(E_mat)
col_total <- colSums(E_mat)
grand_total <- sum(E_mat)

cat('
<p style="text-align:center; font-weight:bold; margin-bottom:5px;">
Tabel 2.6 Frekuensi Harapan (Expected Frequencies)
</p>

<table class="my-table">
  <tr>
    <th>Gender</th>
    <th>Democrat</th>
    <th>Republican</th>
    <th>Independent</th>
    <th>Total</th>
  </tr>
  <tr>
    <td><b>Female</b></td>
    <td>', round(E_mat[1,1],2), '</td>
    <td>', round(E_mat[1,2],2), '</td>
    <td>', round(E_mat[1,3],2), '</td>
    <td><b>', round(row_total[1],2), '</b></td>
  </tr>
  <tr>
    <td><b>Male</b></td>
    <td>', round(E_mat[2,1],2), '</td>
    <td>', round(E_mat[2,2],2), '</td>
    <td>', round(E_mat[2,3],2), '</td>
    <td><b>', round(row_total[2],2), '</b></td>
  </tr>
  <tr>
    <td><b>Total</b></td>
    <td><b>', round(col_total[1],2), '</b></td>
    <td><b>', round(col_total[2],2), '</b></td>
    <td><b>', round(col_total[3],2), '</b></td>
    <td><b>', round(grand_total,2), '</b></td>
  </tr>
</table>

<p><b>Verifikasi:</b> Semua E<sub>ij</sub> > 5? ', all(E_mat > 5), '</p>
<p>(Asumsi chi-square terpenuhi)</p>
')

Tabel 2.6 Frekuensi Harapan (Expected Frequencies)

Gender Democrat Republican Independent Total
Female 456.95 297.43 602.62 1357
Male 368.05 239.57 485.38 1093
Total 825 537 1088 2450

Verifikasi: Semua Eij > 5? TRUE

(Asumsi chi-square terpenuhi)

Interpretasi: Semua frekuensi harapan > 5, sehingga asumsi uji chi-square terpenuhi.

2.3.3 Uji Chi-Square

Hipotesis:

  • Hipotesis Nol \(H_0\): Gender dan identifikasi partai politik independen
  • Hipotesis Alternatif \(H_1\): Terdapat hubungan antara gender dan identifikasi partai politik

Statistik Uji:

\[\chi^2 = \sum_{i=1}^{2}\sum_{j=1}^{3} \frac{(O_{ij} - E_{ij})^2}{E_{ij}}\] \[df = (r-1)(c-1) = (2-1)(3-1) = 2\]

# Output uji chi-square
cat('
<p style="text-align:center; font-weight:bold; margin-bottom:5px;">
Tabel 2.7 Uji Chi-Square Independensi (2×3 Keseluruhan)
</p>

<table class="my-table">
  <tr><th>Komponen</th><th>Nilai</th></tr>
  <tr><td>Statistik χ²</td><td>', round(chi2_full$statistic,4), '</td></tr>
  <tr><td>df</td><td>', chi2_full$parameter, '</td></tr>
  <tr><td>p-value</td><td>', formatC(chi2_full$p.value, format="f", digits=6), '</td></tr>
  <tr><td>Keputusan</td><td>',
      ifelse(chi2_full$p.value < 0.05, "Tolak H₀", "Gagal Tolak H₀"),
  '</td></tr>
</table>
')

Tabel 2.7 Uji Chi-Square Independensi (2×3 Keseluruhan)

Komponen Nilai
Statistik χ² 12.5693
df 2
p-value 0.001865
Keputusan Tolak H₀
# Kontribusi tiap sel
kontribusi <- (tabel2 - E_mat)^2 / E_mat

row_total <- rowSums(kontribusi)
col_total <- colSums(kontribusi)
grand_total <- sum(kontribusi)

cat('
<p style="text-align:center; font-weight:bold; margin-top:15px; margin-bottom:5px;">
Tabel 2.8 Kontribusi Setiap Sel terhadap χ²
</p>

<table class="my-table">
  <tr>
    <th>Gender</th>
    <th>Democrat</th>
    <th>Republican</th>
    <th>Independent</th>
    <th>Total</th>
  </tr>
  <tr>
    <td><b>Female</b></td>
    <td>', round(kontribusi[1,1],4), '</td>
    <td>', round(kontribusi[1,2],4), '</td>
    <td>', round(kontribusi[1,3],4), '</td>
    <td><b>', round(row_total[1],4), '</b></td>
  </tr>
  <tr>
    <td><b>Male</b></td>
    <td>', round(kontribusi[2,1],4), '</td>
    <td>', round(kontribusi[2,2],4), '</td>
    <td>', round(kontribusi[2,3],4), '</td>
    <td><b>', round(row_total[2],4), '</b></td>
  </tr>
  <tr>
    <td><b>Total</b></td>
    <td><b>', round(col_total[1],4), '</b></td>
    <td><b>', round(col_total[2],4), '</b></td>
    <td><b>', round(col_total[3],4), '</b></td>
    <td><b>', round(grand_total,4), '</b></td>
  </tr>
</table>
')

Tabel 2.8 Kontribusi Setiap Sel terhadap χ²

Gender Democrat Republican Independent Total
Female 3.1686 2.1746 0.2642 5.6074
Male 3.9339 2.6999 0.3281 6.9618
Total 7.1025 4.8745 0.5923 12.5693

Interpretasi: Karena p-value < 0,05, maka keputusan yang diambil adalah menolak hipotesis nol (\(H_0\)). Artinya, terdapat hubungan yang signifikan antara jenis kelamin dan preferensi politik.

2.3.4 Residual Pearson dan Standardized Residual

Residual Pearson: \[e_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}}\]

Standardized Residual (Adjusted): \[r_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}(1 - p_{i\cdot})(1 - p_{\cdot j})}}\]

# Residual Pearson
pearson_res <- chi2_full$residuals

# Standardized residual
row_total <- rowSums(tabel2)
col_total <- colSums(tabel2)
N2 <- sum(tabel2)

std_res <- matrix(0, nrow = 2, ncol = 3,
                  dimnames = dimnames(tabel2))

for (i in 1:2) {
  for (j in 1:3) {
    std_res[i, j] <- (tabel2[i, j] - E_mat[i, j]) /
      sqrt(E_mat[i, j] * (1 - row_total[i] / N2) * (1 - col_total[j] / N2))
  }
}

# Flatten
Sel <- c("Female-Democrat", "Female-Republican", "Female-Independent",
         "Male-Democrat", "Male-Republican", "Male-Independent")

O <- as.vector(t(tabel2))
E <- round(as.vector(t(E_mat)), 2)
Pearson <- round(as.vector(t(pearson_res)), 4)
Std <- round(as.vector(t(std_res)), 4)

# BUILD HTML
html <- '
<p style="text-align:center; font-weight:bold; margin-bottom:5px;">
Tabel 2.9 Residual Pearson dan Standardized Residual per Sel
</p>

<table class="my-table">
  <tr>
    <th>Sel</th>
    <th>O</th>
    <th>E</th>
    <th>Residual Pearson</th>
    <th>Std. Residual</th>
    <th>Signifikan?</th>
  </tr>
'

for (i in 1:length(Sel)) {
  signif_flag <- ifelse(abs(Std[i]) > 1.96, "Ya (*)", "Tidak")
  highlight <- ifelse(abs(Std[i]) > 1.96,
                      "style=\"background-color:#ffeeba; font-weight:bold;\"",
                      "")
  
  html <- paste0(html,
    '<tr ', highlight, '>
      <td>', Sel[i], '</td>
      <td>', O[i], '</td>
      <td>', E[i], '</td>
      <td>', Pearson[i], '</td>
      <td>', Std[i], '</td>
      <td>', signif_flag, '</td>
    </tr>'
  )
}

html <- paste0(html, '
</table>

<p>Nilai > |1.96| menunjukkan kontribusi signifikan (α = 0.05)</p>
')

cat(html)

Tabel 2.9 Residual Pearson dan Standardized Residual per Sel

Sel O E Residual Pearson Std. Residual Signifikan?
Female-Democrat 495 456.95 1.7801 3.2724 Ya (*)
Female-Republican 272 297.43 -1.4747 -2.4986 Ya (*)
Female-Independent 590 602.62 -0.514 -1.0322 Tidak
Male-Democrat 330 368.05 -1.9834 -3.2724 Ya (*)
Male-Republican 265 239.57 1.6431 2.4986 Ya (*)
Male-Independent 498 485.38 0.5728 1.0322 Tidak

Nilai > |1.96| menunjukkan kontribusi signifikan (α = 0.05)

Interpretasi: Sel dengan |standardized residual| > 1.96 memberikan kontribusi signifikan terhadap statistik chi-square, dengan sel yang signifikan adalah: Female-Democrat, Female-Republican, Male-Democrat, Male-Republican.

2.3.5 Partisi Chi-Square

Partisi chi-square membagi statistik uji total menjadi komponen ortogonal yang dapat diinterpretasikan secara terpisah.

2.3.5.1 Partisi 1: Democrat vs Republican

# Partisi 1: hanya Democrat dan Republican
tabel2_DR <- tabel2[, c("Democrat", "Republican")]
chi2_DR <- chisq.test(tabel2_DR, correct = FALSE)

cat("=== Partisi 1: Democrat vs Republican ===\n")
## === Partisi 1: Democrat vs Republican ===
print(tabel2_DR)
##        Democrat Republican
## Female      495        272
## Male        330        265
cat(sprintf("\nFrekuensi Harapan:\n"))
## 
## Frekuensi Harapan:
print(round(chi2_DR$expected, 2))
##        Democrat Republican
## Female   464.59     302.41
## Male     360.41     234.59
cat(sprintf("\nχ²_DR = %.4f, df = %d, p-value = %.6f\n",
            chi2_DR$statistic, chi2_DR$parameter, chi2_DR$p.value))
## 
## χ²_DR = 11.5545, df = 1, p-value = 0.000676
cat(sprintf("Keputusan: %s H0 pada α = 0.05\n",
            ifelse(chi2_DR$p.value < 0.05, "Tolak", "Gagal Tolak")))
## Keputusan: Tolak H0 pada α = 0.05

Interpretasi: Karena p-value < 0,05, maka keputusan yang diambil adalah menolak hipotesis nol (\(H_0\)). Artinya, terdapat perbedaan distribusi yang signifikan antara kelompok Democrat dan Republican berdasarkan jenis kelamin.

2.3.5.2 Partisi 2: (Democrat + Republican) vs Independent

# Partisi 2: (Democrat + Republican) vs Independent
# Gabungkan Democrat dan Republican
tabel2_DRvsI <- cbind(
  "Dem+Rep" = rowSums(tabel2[, c("Democrat", "Republican")]),
  "Independent" = tabel2[, "Independent"]
)

chi2_DRvsI <- chisq.test(tabel2_DRvsI, correct = FALSE)

cat("=== Partisi 2: (Democrat + Republican) vs Independent ===\n")
## === Partisi 2: (Democrat + Republican) vs Independent ===
print(tabel2_DRvsI)
##        Dem+Rep Independent
## Female     767         590
## Male       595         498
cat(sprintf("\nFrekuensi Harapan:\n"))
## 
## Frekuensi Harapan:
print(round(chi2_DRvsI$expected, 2))
##        Dem+Rep Independent
## Female  754.38      602.62
## Male    607.62      485.38
cat(sprintf("\nχ²_(DR)vsI = %.4f, df = %d, p-value = %.6f\n",
            chi2_DRvsI$statistic, chi2_DRvsI$parameter, chi2_DRvsI$p.value))
## 
## χ²_(DR)vsI = 1.0654, df = 1, p-value = 0.301979
cat(sprintf("Keputusan  : %s H0 pada α = 0.05\n",
            ifelse(chi2_DRvsI$p.value < 0.05, "Tolak", "Gagal Tolak")))
## Keputusan  : Gagal Tolak H0 pada α = 0.05

Interpretasi: Karena p-value >= 0,05, maka keputusan yang diambil adalah gagal menolak hipotesis nol (\(H_0\)). Artinya, tidak terdapat perbedaan distribusi yang signifikan antara kelompok (Democrat + Republican) dan Independent berdasarkan jenis kelamin.

2.3.5.3 Perbandingan Partisi dengan Chi-Square Keseluruhan

# Hitung penjumlahan
chi2_sum <- chi2_DR$statistic + chi2_DRvsI$statistic
df_sum   <- chi2_DR$parameter + chi2_DRvsI$parameter
p_sum    <- pchisq(chi2_sum, df_sum, lower.tail = FALSE)

html <- paste0(
'<p style="text-align:center; font-weight:bold; margin-bottom:5px;">
Tabel 2.10 Perbandingan Chi-Square Keseluruhan dan Hasil Partisi
</p>

<table class="my-table">
  <tr>
    <th>Metode</th>
    <th>χ²</th>
    <th>df</th>
    <th>p-value</th>
    <th>Keputusan</th>
  </tr>

  <tr style="background-color:#dfe6e9; font-weight:bold;">
    <td>Chi-Square Keseluruhan</td>
    <td>', round(chi2_full$statistic,4), '</td>
    <td>', chi2_full$parameter, '</td>
    <td>', formatC(chi2_full$p.value, format="f", digits=6), '</td>
    <td>Tolak H₀</td>
  </tr>

  <tr>
    <td>Partisi 1: Democrat vs Republican</td>
    <td>', round(chi2_DR$statistic,4), '</td>
    <td>', chi2_DR$parameter, '</td>
    <td>', formatC(chi2_DR$p.value, format="f", digits=6), '</td>
    <td>Tolak H₀</td>
  </tr>

  <tr>
    <td>Partisi 2: (Dem+Rep) vs Independent</td>
    <td>', round(chi2_DRvsI$statistic,4), '</td>
    <td>', chi2_DRvsI$parameter, '</td>
    <td>', formatC(chi2_DRvsI$p.value, format="f", digits=6), '</td>
    <td>Gagal Tolak H₀</td>
  </tr>

  <tr style="font-style:italic; background-color:#f8f9fa;">
    <td>Jumlah Partisi</td>
    <td>', round(chi2_sum,4), '</td>
    <td>', df_sum, '</td>
    <td>', formatC(p_sum, format="f", digits=6), '</td>
    <td>—</td>
  </tr>

</table>

<p><b>Catatan:</b> Perbedaan kecil dapat terjadi karena partisi tidak selalu ortogonal sempurna.</p>'
)

cat(html)

Tabel 2.10 Perbandingan Chi-Square Keseluruhan dan Hasil Partisi

Metode χ² df p-value Keputusan
Chi-Square Keseluruhan 12.5693 2 0.001865 Tolak H₀
Partisi 1: Democrat vs Republican 11.5545 1 0.000676 Tolak H₀
Partisi 2: (Dem+Rep) vs Independent 1.0654 1 0.301979 Gagal Tolak H₀
Jumlah Partisi 12.62 2 0.001818

Catatan: Perbedaan kecil dapat terjadi karena partisi tidak selalu ortogonal sempurna.

2.3.6 Visualisasi

par(mfrow = c(2, 2))

# 1. Mosaic Plot
mosaicplot(tabel2,
           main = "Mosaic Plot: Gender vs Partai Politik",
           color = c("#e74c3c", "#3498db", "#2ecc71"),
           xlab = "Gender",
           ylab = "Identifikasi Partai",
           cex.axis = 0.9,
           las = 1)

# 2. Bar plot proporsi per gender
prop_gender <- prop.table(tabel2, margin = 1)
barplot(t(prop_gender),
        beside = FALSE,
        col = c("#e74c3c", "#3498db", "#2ecc71"),
        main = "Proporsi Afiliasi Partai per Gender",
        xlab = "Gender",
        ylab = "Proporsi",
        legend.text = colnames(tabel2),
        args.legend = list(x = "topright", bty = "n", cex = 0.8),
        ylim = c(0, 1.2),
        border = NA)

# 3. Grouped bar chart
prop_partai <- prop.table(tabel2, margin = 2)
barplot(t(prop_gender) * 100,
        beside = TRUE,
        col = c("#e74c3c", "#3498db", "#2ecc71"),
        main = "Persentase Afiliasi Partai per Gender",
        xlab = "Gender",
        ylab = "Persentase (%)",
        legend.text = colnames(tabel2),
        args.legend = list(x = "topright", bty = "n", cex = 0.8),
        border = NA,
        ylim = c(0, 60))

# 4. Heatmap standardized residual
std_res_mat <- std_res
image(1:ncol(std_res_mat), 1:nrow(std_res_mat),
      t(abs(std_res_mat)),
      col = colorRampPalette(c("white", "#f39c12", "#e74c3c"))(20),
      xlab = "Partai", ylab = "Gender",
      main = "Heatmap |Standardized Residual|",
      axes = FALSE)
axis(1, at = 1:3, labels = colnames(std_res_mat))
axis(2, at = 1:2, labels = rownames(std_res_mat))
for (i in 1:ncol(std_res_mat)) {
  for (j in 1:nrow(std_res_mat)) {
    text(i, j, round(std_res_mat[j, i], 2), cex = 1.2, font = 2)
  }
}
box()

par(mfrow = c(1, 1))

2.3.7 Kategori yang Paling Berkontribusi

Kategori yang paling berkontribusi terhadap hubungan antar variabel ditentukan berdasarkan nilai pearson residual absolut terbesar, yaitu sel dengan deviasi paling besar antara frekuensi observasi dan ekspektasi.

# Data kontribusi
Sel <- c("Female-Democrat", "Female-Republican", "Female-Independent",
         "Male-Democrat", "Male-Republican", "Male-Independent")

Std <- round(as.vector(t(std_res)), 4)
Chi <- round(as.vector(t(kontribusi)), 4)

kontrib_df <- data.frame(Sel, Std, Chi)

# Urutkan berdasarkan |Std|
kontrib_df <- kontrib_df[order(-abs(kontrib_df$Std)), ]

# BUILD HTML
html <- '
<p style="text-align:center; font-weight:bold; margin-bottom:5px;">
Tabel 2.11 Ranking Kontribusi Setiap Sel terhadap Asosiasi
</p>

<table class="my-table">
  <tr>
    <th>Sel</th>
    <th>Standardized Residual</th>
    <th>Kontribusi χ²</th>
  </tr>
'

for (i in 1:nrow(kontrib_df)) {
  signif_flag <- abs(kontrib_df$Std[i]) > 1.96
  highlight <- ifelse(signif_flag,
                      "style=\"background-color:#ffeeba; font-weight:bold;\"",
                      "")
  
  html <- paste0(html,
    '<tr ', highlight, '>
      <td>', kontrib_df$Sel[i], '</td>
      <td>', kontrib_df$Std[i], '</td>
      <td>', kontrib_df$Chi[i], '</td>
    </tr>'
  )
}

html <- paste0(html, '
</table>
')

cat(html)

Tabel 2.11 Ranking Kontribusi Setiap Sel terhadap Asosiasi

Sel Standardized Residual Kontribusi χ²
Female-Democrat 3.2724 3.1686
Male-Democrat -3.2724 3.9339
Female-Republican -2.4986 2.1746
Male-Republican 2.4986 2.6999
Female-Independent -1.0322 0.2642
Male-Independent 1.0322 0.3281

Interpretasi:

Berdasarkan standardized residual dan hasil partisi:

  • Republican merupakan kategori yang paling membedakan gender: perempuan under-represented, laki-laki over-represented (atau sebaliknya tergantung tanda residual).
  • Partisi 1 (Democrat vs Republican) menghasilkan chi-square yang signifikan, menunjukkan bahwa perbedaan gender terutama terletak pada preferensi antara dua partai besar.
  • Partisi 2 (Dem+Rep vs Independent) tidak signifikan, artinya gender tidak membedakan secara bermakna apakah seseorang memilih jalur independen vs dua partai besar.

2.3.8 Kesimpulan

  1. Uji chi-square keseluruhan: \(\chi^2 = 12.5693\), \(df = 2\), \(p = 0.001865\)Tolak \(H_0\). Terdapat asosiasi signifikan antara gender dan identifikasi partai politik.

  2. Residual analisis: Sel-sel yang paling berkontribusi terhadap asosiasi adalah yang terkait dengan afiliasi Republican dan Democrat, bukan Independent.

  3. Partisi chi-square:

    • Democrat vs Republican: signifikan (\(p < 0.05\)) — gender membedakan pilihan antara dua partai besar.
    • (Dem+Rep) vs Independent: tidak signifikan (\(p > 0.05\)) — gender tidak membedakan kecenderungan untuk menjadi Independen.
  4. Kesimpulan substantif: Asosiasi antara gender dan afiliasi partai terutama didorong oleh perbedaan preferensi antara Partai Demokrat dan Republik. Perempuan cenderung lebih ke Demokrat, sementara laki-laki relatif lebih ke Republik — konsisten dengan “gender gap” yang sering ditemukan dalam literatur ilmu politik.


2.4 Kesimpulan Umum

Berikut merupakan kesimpulan umum yang diperoleh dari kasus 1 dan 2:

Tabel 2.12 Perbandingan Kasus 1 dan Kasus 2

Aspek Kasus 1 (2×2) Kasus 2 (2×3)
Metode utama Z-test, χ², G², Fisher χ², Residual, Partisi
Hasil uji Signifikan (semua metode) Signifikan (keseluruhan)
Ukuran asosiasi RD, RR, OR Standardized Residual
Kesimpulan Merokok → meningkatkan risiko kanker paru Gender → berhubungan dengan afiliasi partai

Kedua kasus menunjukkan bahwa tabel kontingensi dan inferensi statistik yang tepat mampu mengungkap pola asosiasi antar variabel kategorik secara sistematis dan reproducible.

3 Tugas 3: Generalized Linear Model

3.1 Regresi Logistik Biner

3.1.1 Deskripsi Kasus

Penyakit gagal jantung (heart failure) merupakan kondisi ketika jantung tidak mampu memompa darah secara efektif untuk memenuhi kebutuhan tubuh. Kondisi ini dapat meningkatkan risiko komplikasi serius hingga kematian. Oleh karena itu, diperlukan analisis untuk mengidentifikasi faktor-faktor yang berpengaruh terhadap risiko kematian pasien gagal jantung.

Pada penelitian ini digunakan data pasien gagal jantung yang memuat informasi karakteristik klinis dan hasil pemeriksaan laboratorium. Analisis dilakukan menggunakan regresi logistik biner untuk memodelkan hubungan antara karakteristik pasien dengan kemungkinan terjadinya kematian. Model ini diharapkan dapat membantu mengidentifikasi faktor-faktor paling berpengaruh dalam menentukan risiko kematian pasien gagal jantung.

Variabel yang digunakan ditampilkan pada tabel berikut:

Tabel 3.1 Variabel Regresi Logistik Biner

Variabel Tipe Keterangan
DEATH_EVENT Kategorik (Biner) Status kematian pasien (0 = hidup, 1 = meninggal)
age Numerik Usia pasien (tahun)
ejection_fraction Numerik Persentase darah yang dipompa keluar oleh jantung (%)
serum_creatinine Numerik Kadar kreatinin dalam darah (mg/dL)
serum_sodium Numerik Kadar natrium dalam darah (mEq/L)
time Numerik Lama waktu observasi pasien (hari)

3.1.2 Preparasi Data

Tahap awal dalam analisis ini adalah preparasi data yang bertujuan untuk memastikan dataset siap digunakan dalam analisis.

#-------------------------
# LOAD PACKAGE
#-------------------------
library(readr)
library(dplyr)
library(ggplot2)
library(tidyr)
library(caret)
## Warning: package 'caret' was built under R version 4.4.3
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
## The following object is masked from 'package:survival':
## 
##     cluster
library(car)
## Warning: package 'car' was built under R version 4.4.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.4.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
library(pscl)
## Warning: package 'pscl' was built under R version 4.4.3
## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.
library(ResourceSelection)
## Warning: package 'ResourceSelection' was built under R version 4.4.3
## ResourceSelection 0.3-6   2023-06-27
library(pROC)
## Warning: package 'pROC' was built under R version 4.4.3
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
library(tibble)
library(readxl)
## Warning: package 'readxl' was built under R version 4.4.3
library(knitr)
library(kableExtra)
library(rcompanion)
## Warning: package 'rcompanion' was built under R version 4.4.3
#-------------------------
# IMPORT DATA
#-------------------------
data1 <- read_excel("C:/Users/Asus/Downloads/DATA ADK.xlsx",
                    sheet = "Binary")

data1$DEATH_EVENT <- as.factor(data1$DEATH_EVENT)

3.1.3 Eksplorasi Data

Eksplorasi data dilakukan untuk memberikan gambaran awal mengenai karakteristik data serta hubungan antara variabel penjelas dengan variabel respon.

#-------------------------
# DESKRIPTIF 
#-------------------------

stat_desk <- data1 %>%
  group_by(DEATH_EVENT) %>%
  summarise(
    n = n(),
    Mean_Age = mean(age),
    Mean_EjectionFraction = mean(ejection_fraction),
    Mean_SerumCreatinine = mean(serum_creatinine),
    Mean_SerumSodium = mean(serum_sodium),
    Mean_Time = mean(time)
  )
#-------------------------
# DISTRIBUSI KELAS
#-------------------------

ggplot(
  data1,
  aes(
    x = DEATH_EVENT,
    fill = DEATH_EVENT
  )
) +
  geom_bar() +
  geom_text(
    stat = "count",
    aes(label = after_stat(count)),
    vjust = -0.5
  ) +
  theme_minimal() +
  labs(
    title = "Distribusi Death Event",
    x = "DEATH_EVENT",
    y = "Frekuensi"
  )

Berdasarkan grafik distribusi:

  • DEATH_EVENT = 0 (tidak meninggal): 203 pasien (≈ 67,9%)
  • DEATH_EVENT = 1 (meninggal): 96 pasien (≈ 32,1%)

Artinya, mayoritas pasien dalam dataset tidak mengalami kematian selama periode observasi. Namun proporsi pasien yang meninggal juga cukup besar (sekitar sepertiga data), sehingga masalah klasifikasi masih relatif seimbang dan tidak terlalu mengalami ketimpangan kelas yang ekstrem.

#-------------------------
# BOXPLOT VARIABEL PREDIKTOR
#-------------------------

data_long <- pivot_longer(
  data1,
  cols = c(
    age,
    ejection_fraction,
    serum_creatinine,
    serum_sodium,
    time
  )
)

ggplot(
  data_long,
  aes(
    x = DEATH_EVENT,
    y = value,
    fill = DEATH_EVENT
  )
) +
  geom_boxplot() +
  facet_wrap(~name, scales = "free") +
  theme_minimal() +
  labs(
    title = "Boxplot Variabel Prediktor terhadap DEATH_EVENT"
  ) 

Hasil eksplorasi data menggunakan boxplot menunjukkan adanya perbedaan distribusi beberapa variabel prediktor berdasarkan status DEATH_EVENT. Kelompok pasien yang mengalami kematian cenderung memiliki usia lebih tinggi, nilai ejection fraction lebih rendah, kadar serum creatinine lebih tinggi, kadar serum sodium lebih rendah, serta waktu observasi yang lebih pendek dibandingkan kelompok yang tidak mengalami kematian. Temuan ini mengindikasikan bahwa variabel-variabel tersebut berpotensi memiliki hubungan dengan kejadian kematian pada pasien gagal jantung dan layak dipertimbangkan dalam proses pemodelan klasifikasi.

3.1.4 Pembagian Data

Sebelum dilakukan pemodelan, dataset dibagi menjadi training data dan testing data dengan proporsi 80:20. Training data digunakan untuk membentuk model, sedangkan testing data digunakan untuk mengevaluasi performa model yang diperoleh.

#-------------------------
# SPLIT DATA
#-------------------------
set.seed(123)

train_index <- createDataPartition(
  data1$DEATH_EVENT,
  p = 0.8,
  list = FALSE
)

train <- data1[train_index, ]
test  <- data1[-train_index, ]

3.1.5 Pengujian Asumsi Multikolinearitas

Uji multikolinearitas untuk mengetahui apakah terdapat hubungan yang kuat antarvariabel prediktor. Multikolinearitas yang tinggi dapat menyebabkan ketidakstabilan dalam estimasi parameter model sehingga interpretasi koefisien regresi menjadi kurang reliabel. Pada penelitian ini, uji multikolinearitas dilakukan menggunakan nilai Variance Inflation Factor (VIF).

Kriteria yang digunakan adalah:

  • Jika VIF < 10, maka tidak terdapat masalah multikolinearitas.

  • Jika VIF > 10, maka terdapat masalah multikolinearitas.

#-------------------------
# MODEL LOGISTIK
#-------------------------
model <- glm(
  DEATH_EVENT ~ age + ejection_fraction + serum_creatinine +
    serum_sodium + time,
  family = binomial(link = "logit"),
  data = train
)

#-------------------------
# MULTIKOLINEARITAS
#-------------------------
vif_result <- vif(model)

df_vif <- data.frame(
  Variabel = names(vif_result),
  VIF      = round(as.numeric(vif_result), 4),
  Status   = ifelse(as.numeric(vif_result) < 10,
                    "✔ Tidak ada multikolinearitas",
                    "⚠ Multikolinearitas")
)

kable(df_vif,
      caption = "Tabel 3.2 Hasil Uji Multikolinearitas (VIF)",
      align = c("c","c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.2 Hasil Uji Multikolinearitas (VIF)
Variabel VIF Status
age 1.0970 ✔ Tidak ada multikolinearitas
ejection_fraction 1.2150 ✔ Tidak ada multikolinearitas
serum_creatinine 1.1208 ✔ Tidak ada multikolinearitas
serum_sodium 1.0267 ✔ Tidak ada multikolinearitas
time 1.0637 ✔ Tidak ada multikolinearitas

Berdasarkan hasil perhitungan Variance Inflation Factor (VIF), seluruh variabel prediktor memiliki nilai VIF kurang dari 10. Hal ini menunjukkan bahwa tidak terdapat permasalahan multikolinearitas antar variabel independen, sehingga model regresi yang dibangun stabil dan layak untuk diinterpretasikan.

3.1.6 Pemodelan Regresi Logistik Biner

Model regresi logistik biner digunakan untuk memodelkan hubungan antara variabel respon yang bersifat dikotomis dengan sejumlah variabel prediktor. Bentuk umum model regresi logistik biner adalah

\[ g(x)=\ln\left(\frac{\pi(x)}{1-\pi(x)}\right) =\beta_0+\beta_1X_1+\beta_2X_2+\cdots+\beta_pX_p \]

dengan:

  • \(\pi(x)\) = peluang suatu observasi masuk ke kategori yang diamati
  • \(X_i\) = variabel prediktor ke-i
  • \(\beta_0\) = merupakan intersep
  • \(\beta_i\) = merupakan koefisien regresi untuk variabel prediktor ke-i.
summary(model)
## 
## Call:
## glm(formula = DEATH_EVENT ~ age + ejection_fraction + serum_creatinine + 
##     serum_sodium + time, family = binomial(link = "logit"), data = train)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        9.155477   6.652360   1.376  0.16874    
## age                0.048619   0.016986   2.862  0.00421 ** 
## ejection_fraction -0.093347   0.019531  -4.779 1.76e-06 ***
## serum_creatinine   0.728157   0.185388   3.928 8.57e-05 ***
## serum_sodium      -0.062650   0.047300  -1.325  0.18533    
## time              -0.017803   0.003096  -5.751 8.87e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 301.20  on 239  degrees of freedom
## Residual deviance: 178.17  on 234  degrees of freedom
## AIC: 190.17
## 
## Number of Fisher Scoring iterations: 6
df_model <- data.frame(
  Variabel = rownames(summary(model)$coefficients),
  Estimate = round(summary(model)$coefficients[, "Estimate"], 4)
)

kable(df_model,
      caption = "Tabel 3.3 Hasil Estimasi Regresi Logistik",
      align = c("c", "c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.3 Hasil Estimasi Regresi Logistik
Variabel Estimate
(Intercept) 9.1555
age 0.0486
ejection_fraction -0.0933
serum_creatinine 0.7282
serum_sodium -0.0626
time -0.0178

3.1.7 Pengujian Hipotesis

3.1.7.1 Uji Simultan (Likelihood Ratio Test)

Uji signifikansi simultan dilakukan menggunakan Likelihood Ratio Test (LRT) untuk mengetahui apakah variabel prediktor secara bersama-sama berpengaruh terhadap variabel respon.

Hipotesis yang digunakan adalah:

  • \(H_0\) : \(\beta_1 = \beta_2 = \beta_3 = \cdots = \beta_k = 0\) (tidak terdapat pengaruh signifikan secara simultan antara variabel prediktor terhadap variabel respon)

  • \(H_1\) : paling sedikit terdapat satu \(\beta_j \neq 0\) (terdapat minimal satu variabel prediktor yang berpengaruh signifikan terhadap variabel respon)

#-------------------------
# UJI SIMULTAN
#-------------------------
model_null <- glm(DEATH_EVENT ~ 1,
                  family = binomial,
                  data = train)

anova(model_null, model, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: DEATH_EVENT ~ 1
## Model 2: DEATH_EVENT ~ age + ejection_fraction + serum_creatinine + serum_sodium + 
##     time
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       239     301.20                          
## 2       234     178.17  5   123.03 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df_simultan <- data.frame(
  `Statistik G` = 123.03,
  df = 5,
  `p-value` = "<0.001"
)

kable(df_simultan,
      caption = "Tabel 3.4 Hasil Uji Simultan Regresi Logistik Biner",
      align = c("c", "c", "c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.4 Hasil Uji Simultan Regresi Logistik Biner
Statistik.G df p.value
123.03 5 <0.001

Berdasarkan hasil uji simultan, diperoleh nilai statistik G sebesar 123,03 dengan derajat bebas (df) sebesar 5 dan p-value < 0,001. Karena p-value < 0,05, maka \(H_0\) ditolak. Dengan demikian, variabel age, ejection_fraction, serum_creatinine, serum_sodium, dan time secara simultan berpengaruh signifikan terhadap DEATH_EVENT.

3.1.7.2 Uji Parsial (Wald Test)

Uji signifikansi parsial dilakukan menggunakan uji Wald untuk mengetahui pengaruh masing-masing variabel prediktor terhadap variabel respon.

Hipotesis (untuk setiap variabel prediktor):

  • \(H_0 : \beta_i = 0\) (variabel ke-\(i\) tidak berpengaruh terhadap variabel respon)

  • \(H_1 : \beta_i \neq 0\) (variabel ke-\(i\) berpengaruh terhadap variabel respon)

# UJI PARSIAL
summary(model)
## 
## Call:
## glm(formula = DEATH_EVENT ~ age + ejection_fraction + serum_creatinine + 
##     serum_sodium + time, family = binomial(link = "logit"), data = train)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        9.155477   6.652360   1.376  0.16874    
## age                0.048619   0.016986   2.862  0.00421 ** 
## ejection_fraction -0.093347   0.019531  -4.779 1.76e-06 ***
## serum_creatinine   0.728157   0.185388   3.928 8.57e-05 ***
## serum_sodium      -0.062650   0.047300  -1.325  0.18533    
## time              -0.017803   0.003096  -5.751 8.87e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 301.20  on 239  degrees of freedom
## Residual deviance: 178.17  on 234  degrees of freedom
## AIC: 190.17
## 
## Number of Fisher Scoring iterations: 6
df_parsial <- data.frame(
  Variabel = rownames(summary(model)$coefficients),
  Estimate = round(summary(model)$coefficients[, "Estimate"], 4),
  `p-value` = round(summary(model)$coefficients[, "Pr(>|z|)"], 4)
)

kable(
  df_parsial,
  caption = "Tabel 3.5 Hasil Uji Parsial Regresi Logistik Biner",
  align = c("c", "c", "c"),
  row.names = FALSE
) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a") %>%
  row_spec(5, background = "#fdeaea")
Tabel 3.5 Hasil Uji Parsial Regresi Logistik Biner
Variabel Estimate p.value
(Intercept) 9.1555 0.1687
age 0.0486 0.0042
ejection_fraction -0.0933 0.0000
serum_creatinine 0.7282 0.0001
serum_sodium -0.0626 0.1853
time -0.0178 0.0000

Berdasarkan hasil uji parsial regresi logistik, diperoleh bahwa variabel age, ejection_fraction, serum_creatinine, dan time memiliki nilai p-value kurang dari 0,05 sehingga berpengaruh signifikan terhadap kejadian kematian (DEATH_EVENT). Sementara itu, variabel serum_sodium memiliki nilai p-value sebesar 0,1853 (> 0,05), sehingga tidak berpengaruh signifikan terhadap DEATH_EVENT. Dengan demikian, pada taraf signifikansi 5%, terdapat empat variabel yang secara parsial berpengaruh signifikan terhadap peluang terjadinya DEATH_EVENT, yaitu age, ejection_fraction, serum_creatinine, dan time.

3.1.8 Odds Ratio

Odds Ratio (OR) digunakan untuk mengukur besar pengaruh variabel prediktor terhadap variabel respon dalam model regresi logistik. Nilai Odds Ratio diperoleh dengan melakukan eksponensial terhadap koefisien regresi.

Secara matematis, Odds Ratio dinyatakan sebagai:

\[ OR = e^{\beta_j} \]

dengan \(\beta_j\) adalah koefisien regresi dari variabel prediktor ke-j.

Pada regresi logistik biner, Odds Ratio mengukur perubahan peluang terjadinya suatu kejadian (kategori 1 dibandingkan 0).

Interpretasi Odds Ratio:

  • Jika \(OR > 1\), variabel prediktor meningkatkan peluang kejadian.
  • Jika \(OR < 1\), variabel prediktor menurunkan peluang kejadian.
  • Jika \(OR = 1\), variabel tidak berpengaruh terhadap peluang kejadian.
#-------------------------
# ODDS RATIO
#-------------------------
exp(coef(model))
##       (Intercept)               age ejection_fraction  serum_creatinine 
##      9466.1478575         1.0498206         0.9108778         2.0712604 
##      serum_sodium              time 
##         0.9392723         0.9823541
exp(confint(model))
## Waiting for profiling to be done...
##                        2.5 %       97.5 %
## (Intercept)       0.02499116 6.166002e+09
## age               1.01638798 1.086742e+00
## ejection_fraction 0.87443604 9.443976e-01
## serum_creatinine  1.45015569 3.068003e+00
## serum_sodium      0.85436097 1.029564e+00
## time              0.97599478 9.879781e-01
df_or <- data.frame(
  Variabel = c("Age", "Ejection Fraction", "Serum Creatinine", "Serum Sodium", "Time"),
  OR = round(exp(coef(model))[-1], 4),
  CI_lower = round(exp(confint(model))[-1,1], 4),
  CI_upper = round(exp(confint(model))[-1,2], 4)
)
## Waiting for profiling to be done...
## Waiting for profiling to be done...
kable(df_or,
      caption = "Tabel 3.6 Odds Ratio dan 95% Confidence Interval",
      align = c("c","c","c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a") %>%
  row_spec(4, background = "#fdeaea")
Tabel 3.6 Odds Ratio dan 95% Confidence Interval
Variabel OR CI_lower CI_upper
Age 1.0498 1.0164 1.0867
Ejection Fraction 0.9109 0.8744 0.9444
Serum Creatinine 2.0713 1.4502 3.0680
Serum Sodium 0.9393 0.8544 1.0296
Time 0.9824 0.9760 0.9880

Contoh Interpretasi:

  • Variabel age memiliki nilai OR sebesar 1,050. Hal ini menunjukkan bahwa setiap peningkatan usia sebesar satu tahun akan meningkatkan odds terjadinya kematian sebesar 1,050 kali atau sekitar 4,98%, dengan asumsi variabel lain tetap. Selang kepercayaan 95% tidak memuat nilai 1 sehingga pengaruh usia signifikan terhadap kejadian kematian.

  • Variabel ejection fraction memiliki nilai OR sebesar 0,911. Artinya, setiap peningkatan ejection fraction sebesar 1% akan menurunkan odds terjadinya kematian sebesar 8,9%. Nilai OR yang kurang dari satu menunjukkan bahwa ejection fraction merupakan faktor protektif terhadap kejadian kematian. Selang kepercayaan 95% yang seluruhnya berada di bawah 1 mengindikasikan bahwa pengaruh tersebut signifikan.

Berdasarkan nilai odds ratio, serum_creatinine merupakan variabel yang memiliki pengaruh paling kuat terhadap kejadian kematian karena memiliki OR terbesar, yaitu 2,071. Sementara itu, ejection_fraction dan time berperan sebagai faktor protektif karena memiliki OR kurang dari satu. Hasil ini mengindikasikan bahwa peningkatan fungsi jantung (ejection fraction) serta lamanya waktu bertahan selama masa observasi berkaitan dengan penurunan peluang kematian pasien gagal jantung. Di sisi lain, peningkatan usia dan kadar serum creatinine cenderung meningkatkan peluang terjadinya kematian.

3.1.9 Evaluasi Model

3.1.9.1 Uji Hosmer-Lemeshow

Uji Hosmer-Lemeshow digunakan untuk mengevaluasi kesesuaian model regresi logistik dengan data observasi. Hipotesis yang digunakan adalah:

\(H_0\) : Model sesuai dengan data (good fit).

\(H_1\) : Model tidak sesuai dengan data.

#-------------------------
# HOSMER-LEMESHOW
#-------------------------
hoslem.test(
  as.numeric(train$DEATH_EVENT) - 1,
  fitted(model)
)
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  as.numeric(train$DEATH_EVENT) - 1, fitted(model)
## X-squared = 5.8366, df = 8, p-value = 0.6655
df_hl <- data.frame(
  `Chi-Square` = 5.8366,
  df = 8,
  `p-value` = 0.6655
)

kable(df_hl,
      caption = "Tabel 3.7 Uji Hosmer-Lemeshow Goodness of Fit",
      align = c("c","c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.7 Uji Hosmer-Lemeshow Goodness of Fit
Chi.Square df p.value
5.8366 8 0.6655

Berdasarkan hasil pengujian diperoleh nilai statistik Hosmer-Lemeshow sebesar 5.8366 dengan p-value sebesar 0.6655. Karena p-value > 0.05 maka gagal menolak \(H_0\).

Dengan demikian dapat disimpulkan bahwa model regresi logistik yang dibentuk memiliki kesesuaian yang baik dengan data observasi sehingga layak digunakan untuk analisis lebih lanjut.

3.1.9.2 Koefisisen Determinasi

Koefisien determinasi digunakan untuk mengevaluasi kemampuan model regresi logistik dalam menjelaskan variabilitas variabel respon.

#-------------------------
# PSEUDO R2
#-------------------------
nagelkerke(model)
## $Models
##                                                                                                                                
## Model: "glm, DEATH_EVENT ~ age + ejection_fraction + serum_creatinine + serum_sodium + time, binomial(link = \"logit\"), train"
## Null:  "glm, DEATH_EVENT ~ 1, binomial(link = \"logit\"), train"                                                               
## 
## $Pseudo.R.squared.for.model.vs.null
##                              Pseudo.R.squared
## McFadden                             0.408475
## Cox and Snell (ML)                   0.401083
## Nagelkerke (Cragg and Uhler)         0.561017
## 
## $Likelihood.ratio.test
##  Df.diff LogLik.diff  Chisq    p.value
##       -5     -61.516 123.03 7.1508e-25
## 
## $Number.of.observations
##           
## Model: 240
## Null:  240
## 
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
## 
## $Warnings
## [1] "None"
df_r2 <- data.frame(
  Indikator = c("McFadden",
                "Cox and Snell (ML)",
                "Nagelkerke (Cragg and Uhler)"),
  Nilai = c(0.408475,
            0.401083,
            0.561017)
)

kable(df_r2,
      caption = "Tabel 3.8 Nilai Pseudo R-Square Model Regresi Logistik Biner",
      align = c("c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.8 Nilai Pseudo R-Square Model Regresi Logistik Biner
Indikator Nilai
McFadden 0.408475
Cox and Snell (ML) 0.401083
Nagelkerke (Cragg and Uhler) 0.561017

Nilai pseudo R-square menunjukkan bahwa model memiliki kemampuan penjelasan yang cukup baik. McFadden R² sebesar 0,408 menunjukkan kecocokan model yang baik, Cox and Snell R² sebesar 0,401 menunjukkan kecocokan model yang moderat, dan Nagelkerke R² sebesar 0,561 menunjukkan bahwa model memiliki kemampuan penjelasan yang cukup tinggi terhadap variabel dependen.

3.1.9.3 ROC dan AUC

Selanjutnya, evaluasi dapat dilakukan dengan menggunakan Receiver Operating Characteristic (ROC) curve dan Area Under the Curve (AUC). ROC curve digunakan untuk menggambarkan kemampuan model dalam membedakan antara kelas positif dan negatif pada berbagai nilai ambang (threshold), sedangkan AUC digunakan sebagai ukuran ringkasan dari performa klasifikasi model.

#-------------------------
# PREDIKSI
#-------------------------
prob <- predict(model, newdata = test, type = "response")

prediksi <- ifelse(prob > 0.5, 1, 0)

#-------------------------
# ROC & AUC
#-------------------------
roc_obj <- roc(
  response = as.numeric(test$DEATH_EVENT) - 1,
  predictor = prob,
  quiet = TRUE
)

auc_value <- auc(roc_obj)

auc_value
## Area under the curve: 0.8724
roc_df <- tibble(
  Specificity = roc_obj$specificities,
  Sensitivity = roc_obj$sensitivities
) %>%
  mutate(FPR = 1 - Specificity)

ggplot(roc_df, aes(x = FPR, y = Sensitivity)) +
  geom_line(linewidth = 1.4, color = "#2f7a78") +
  geom_abline(linetype = "dashed", color = "grey60") +
  
  labs(
    title = "ROC Curve",
    subtitle = paste0("AUC = ", round(auc_value, 4)),
    x = "False Positive Rate (1 - Specificity)",
    y = "True Positive Rate (Sensitivity)"
  ) +
  
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

Kurva ROC digunakan untuk mengevaluasi kemampuan model dalam membedakan pasien yang mengalami kematian dan yang tidak mengalami kematian. Hasil analisis menunjukkan nilai Area Under Curve (AUC) sebesar 0.8724. Berdasarkan kriteria klasifikasi AUC, nilai tersebut berada pada kategori baik (good classification). Dengan demikian model regresi logistik yang dibangun memiliki kemampuan diskriminasi yang baik dalam memprediksi kejadian kematian pasien gagal jantung.

opt <- coords(
  roc_obj,
  x = "best",
  best.method = "youden",
  ret = c("threshold", "sensitivity", "specificity")
)

opt
##   threshold sensitivity specificity
## 1 0.3827618   0.7894737        0.85
cutoff <- as.numeric(opt$threshold)

df_roc <- data.frame(
  Threshold = 0.3827618,
  Sensitivity = 0.7894737,
  Specificity = 0.85
)

kable(df_roc,
      caption = "Tabel 3.9 Treshold Optimal",
      align = c("c","c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.9 Treshold Optimal
Threshold Sensitivity Specificity
0.3827618 0.7894737 0.85

Berdasarkan hasil analisis ROC, diperoleh nilai treshold optimal sebesar 0,3828 dengan sensitivity sebesar 0,7895 dan specificity sebesar 0,8500. Hal ini menunjukkan bahwa model memiliki kemampuan yang baik dalam membedakan kelas positif dan negatif.

3.1.9.4 Confusion Matrix dan Akurasi

Kinerja klasifikasi model regresi logistik dievaluasi menggunakan confusion matrix dan nilai akurasi. Confusion matrix digunakan untuk membandingkan hasil prediksi model dengan kondisi aktual sehingga dapat diketahui jumlah observasi yang diklasifikasikan dengan benar maupun salah.

#-------------------------
# CONFUSION MATRIX
#-------------------------
pred_opt <- ifelse(prob >= cutoff, 1, 0)

table(
  Actual = test$DEATH_EVENT,
  Predicted = pred_opt
)
##       Predicted
## Actual  0  1
##      0 34  6
##      1  4 15
df_cm <- data.frame(
  Actual_0 = c(34, 6),
  Actual_1 = c(4, 15)
)
rownames(df_cm) <- c("Predicted 0", "Predicted 1")

kable(df_cm,
      caption = "Tabel 3.10 Confusion Matrix Model Regresi Logistik Biner",
      align = c("c","c"),
      row.names = TRUE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.10 Confusion Matrix Model Regresi Logistik Biner
Actual_0 Actual_1
Predicted 0 34 4
Predicted 1 6 15

Interpretasi:

  • 34 pasien yang sebenarnya tidak meninggal (0) berhasil diprediksi tidak meninggal (True Negative).

  • 15 pasien yang sebenarnya meninggal (1) berhasil diprediksi meninggal (True Positive).

  • 6 pasien yang sebenarnya tidak meninggal diprediksi meninggal (False Positive).

  • 4 pasien yang sebenarnya meninggal diprediksi tidak meninggal (False Negative).

confusionMatrix(
  factor(pred_opt),
  factor(test$DEATH_EVENT),
  positive = "1"
)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 34  4
##          1  6 15
##                                           
##                Accuracy : 0.8305          
##                  95% CI : (0.7103, 0.9156)
##     No Information Rate : 0.678           
##     P-Value [Acc > NIR] : 0.006651        
##                                           
##                   Kappa : 0.6223          
##                                           
##  Mcnemar's Test P-Value : 0.751830        
##                                           
##             Sensitivity : 0.7895          
##             Specificity : 0.8500          
##          Pos Pred Value : 0.7143          
##          Neg Pred Value : 0.8947          
##              Prevalence : 0.3220          
##          Detection Rate : 0.2542          
##    Detection Prevalence : 0.3559          
##       Balanced Accuracy : 0.8197          
##                                           
##        'Positive' Class : 1               
## 
df_eval <- data.frame(
  Metrik = c("Accuracy",
             "Kappa",
             "Sensitivity",
             "Specificity",
             "Balanced Accuracy",
             "No Information Rate",
             "P-Value (Acc > NIR)"),
  Nilai = c(0.8305,
            0.6223,
            0.7895,
            0.8500,
            0.8197,
            0.6780,
            0.006651)
)

kable(df_eval,
      caption = "Tabel 3.11 Evaluasi Kinerja Model Regresi Logistik Biner",
      align = c("c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.11 Evaluasi Kinerja Model Regresi Logistik Biner
Metrik Nilai
Accuracy 0.830500
Kappa 0.622300
Sensitivity 0.789500
Specificity 0.850000
Balanced Accuracy 0.819700
No Information Rate 0.678000
P-Value (Acc > NIR) 0.006651

Berdasarkan hasil evaluasi kinerja model regresi logistik, diperoleh nilai accuracy sebesar 0,8305 dengan 95% CI (0,7103–0,9156), yang menunjukkan bahwa model memiliki tingkat ketepatan klasifikasi yang baik. Nilai Kappa sebesar 0,6223 mengindikasikan tingkat kesepakatan yang cukup kuat antara prediksi model dan data aktual. Selain itu, nilai sensitivity sebesar 0,7895 menunjukkan kemampuan model dalam mengidentifikasi kasus positif dengan baik, sedangkan specificity sebesar 0,8500 menunjukkan kemampuan yang baik dalam mengidentifikasi kasus negatif. Nilai balanced accuracy sebesar 0,8197 juga menegaskan bahwa performa model cukup seimbang dalam mengklasifikasikan kedua kelas.

Nilai No Information Rate (0,678) yang lebih rendah dibandingkan accuracy serta p-value 0,006651 (< 0,05) menunjukkan bahwa akurasi model secara signifikan lebih baik dibandingkan model tanpa informasi (model acak/majority class). Secara keseluruhan, model regresi logistik yang dibangun memiliki performa yang baik dalam memprediksi DEATH_EVENT.

3.1.10 Kesimpulan

Berdasarkan hasil analisis regresi logistik biner, variabel age, ejection fraction, serum creatinine, dan time berpengaruh signifikan terhadap kejadian kematian pasien gagal jantung, sedangkan serum sodium tidak berpengaruh signifikan. Model yang dibangun memiliki kesesuaian yang baik dengan data berdasarkan uji Hosmer–Lemeshow (p-value = 0,6655), yang menunjukkan tidak terdapat perbedaan signifikan antara nilai observasi dan prediksi model sehingga model layak digunakan.

Selain itu, model menunjukkan performa klasifikasi yang baik dengan akurasi sebesar 83,05%, Kappa sebesar 0,6223, serta balanced accuracy sebesar 81,97%, yang mengindikasikan bahwa model memiliki kemampuan klasifikasi yang cukup kuat dan seimbang antara kelas positif dan negatif. Nilai AUC sebesar 0,8724 juga menunjukkan bahwa model memiliki kemampuan diskriminasi yang sangat baik dalam membedakan pasien yang mengalami dan tidak mengalami kejadian kematian.

Berdasarkan odds ratio, variabel serum creatinine merupakan faktor yang paling dominan dalam meningkatkan risiko kematian, sedangkan ejection fraction dan time berperan sebagai faktor protektif yang menurunkan risiko kejadian kematian pada pasien gagal jantung.

3.2 Regresi Logistik Multinomial

3.2.1 Deskripsi Kasus

Status kelulusan mahasiswa merupakan salah satu indikator penting dalam evaluasi keberhasilan pendidikan tinggi. Mahasiswa dapat mengalami berbagai kemungkinan status akhir studi, seperti menyelesaikan studi (Graduate), masih melanjutkan perkuliahan (Enrolled), atau tidak menyelesaikan studi (Dropout). Kondisi ini dapat dipengaruhi oleh berbagai faktor akademik maupun administratif, seperti prestasi akademik, kemampuan adaptasi perkuliahan, hingga kepatuhan dalam memenuhi kewajiban administrasi.

Pada penelitian ini digunakan data Student Performance yang memuat informasi karakteristik mahasiswa selama masa studi. Data tersebut digunakan untuk menganalisis faktor-faktor yang berpotensi memengaruhi status akhir mahasiswa dengan tiga kategori keluaran.

Analisis dilakukan menggunakan regresi logistik multinomial untuk memodelkan hubungan antara variabel-variabel prediktor dengan kemungkinan status mahasiswa, yaitu Dropout, Enrolled, dan Graduate.

Variabel yang digunakan ditampilkan pada tabel berikut:

Tabel 3.12 Variabel Regresi Logistik Multinomial

Variabel Tipe Keterangan
Target Kategorik (Multinomial) Status mahasiswa: Dropout, Enrolled, Graduate
Admission grade Numerik Nilai kelulusan/masuk mahasiswa
Tuition fees up to date Kategorik (Biner) Status pembayaran UKT (ya/tidak)
Gender Kategorik (Biner) Jenis kelamin mahasiswa
Scholarship holder Kategorik (Biner) Status penerima beasiswa
Age at enrollment Numerik Usia saat masuk kuliah
Curricular units 1st sem (approved) Numerik Jumlah mata kuliah lulus semester 1
Curricular units 1st sem (grade) Numerik Rata-rata nilai semester 1

3.2.2 Preparasi Data

Tahap awal dalam analisis ini adalah preparasi data yang bertujuan untuk memastikan dataset siap digunakan dalam analisis.

# IMPORT DATA
library(readxl)
library(dplyr)
library(ggplot2)
library(tidyr)
library(caret)
library(nnet)
library(psych)
## Warning: package 'psych' was built under R version 4.4.3
## 
## Attaching package: 'psych'
## The following object is masked from 'package:rcompanion':
## 
##     phi
## The following object is masked from 'package:car':
## 
##     logit
## The following objects are masked from 'package:scales':
## 
##     alpha, rescale
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(knitr)
library(mlogit)
## Registered S3 method overwritten by 'micsr':
##   method          from
##   summary.effects lava
data2 <- read_excel(
  "C:/Users/Asus/Downloads/DATA ADK.xlsx",
  sheet = "Multinomial"
)

data2$Target <- trimws(data2$Target)

data2$Target <- factor(
  data2$Target,
  levels = c("Dropout", "Enrolled", "Graduate")
)

3.2.3 Eksplorasi Data

Eksplorasi data dilakukan untuk memberikan gambaran awal mengenai karakteristik data serta hubungan antara variabel penjelas dengan variabel respon.

# STATISTIK DESKRIPTIF

desc <- data2 %>%
  group_by(Target) %>%
  summarise(
    n = n(),
    Mean_Admission = mean(`Admission grade`),
    Mean_Approved = mean(`Curricular units 1st sem (approved)`),
    Mean_Grade = mean(`Curricular units 1st sem (grade)`),
    Mean_Age = mean(`Age at enrollment`)
  )
# DISTRIBUSI KELAS 
ggplot(data2, aes(x = Target, fill = Target)) +
  geom_bar() +
  geom_text(stat = "count", aes(label = after_stat(count)), vjust = -0.5) +
  theme_minimal() +
  labs(
    title = "Distribusi Status Mahasiswa",
    x = "Status",
    y = "Frekuensi"
  )

Distribusi status mahasiswa menunjukkan bahwa kategori Graduate memiliki jumlah observasi paling tinggi, diikuti oleh Dropout dan Enrolled. Hal ini menunjukkan adanya ketidakseimbangan ringan, namun masih dalam batas yang dapat diterima untuk pemodelan.

# BOXPLOT

data_long <- pivot_longer(
  data2,
  cols = c(
    `Admission grade`,
    `Tuition fees up to date`,
    `Gender`,
    `Scholarship holder`,
    `Age at enrollment`,
    `Curricular units 1st sem (approved)`,
    `Curricular units 1st sem (grade)`
  )
)

ggplot(data_long, aes(x = Target, y = value, fill = Target)) +
  geom_boxplot() +
  facet_wrap(~name, scales = "free") +
  theme_minimal() +
  labs(title = "Boxplot Variabel Penjelas")

Berdasarkan boxplot, terlihat bahwa mahasiswa Graduate memiliki nilai akademik yang lebih tinggi dibandingkan kelompok lainnya, terutama pada variabel Curricular units approved dan grade. Sebaliknya, kelompok Dropout memiliki nilai yang lebih rendah dan usia masuk yang relatif lebih tinggi.

3.2.4 Pembagian Data

Sebelum dilakukan pemodelan, dataset dibagi menjadi training data dan testing data dengan proporsi 80:20. Training data digunakan untuk membentuk model, sedangkan testing data digunakan untuk mengevaluasi performa model yang diperoleh

set.seed(123)

train_index <- createDataPartition(data2$Target, p = 0.8, list = FALSE)

train <- data2[train_index, ]
test  <- data2[-train_index, ]

3.2.5 Pemodelan Regresi Logistik Multinomial

Regresi logistik multinomial digunakan untuk memodelkan hubungan antara variabel respon yang memiliki lebih dari dua kategori dengan sejumlah variabel prediktor. Pada penelitian ini, variabel respon terdiri atas tiga kategori, yaitu Dropout, Enrolled, dan Graduate.

Secara umum, model regresi logistik multinomial dapat dituliskan sebagai berikut:

\[ \log\left(\frac{P(Y = k)}{P(Y = \text{baseline})}\right) = \beta_{0k} + \beta_{1k}X_1 + \beta_{2k}X_2 + \cdots + \beta_{pk}X_p \]

dengan:

  • \(Y\) = variabel respon (status mahasiswa)
  • \(k\) = kategori selain baseline (dalam kasus ini: Enrolled dan Graduate)
  • \(X_i\) = variabel prediktor ke-i
  • \(\beta_{ik}\) = koefisien untuk kategori ke-k

Kategori Dropout digunakan sebagai kategori referensi (baseline), sehingga model akan membandingkan peluang Enrolled dan Graduate terhadap Dropout.

model <- multinom(
  Target ~ `Admission grade` +
    `Tuition fees up to date` +
    `Gender` +
    `Scholarship holder` +
    `Age at enrollment` +
    `Curricular units 1st sem (approved)` +
    `Curricular units 1st sem (grade)`,
  data = train
)
## # weights:  27 (16 variable)
## initial  value 3890.186114 
## iter  10 value 3158.039192
## iter  20 value 2618.443945
## final  value 2597.125469 
## converged
summary(model)
## Call:
## multinom(formula = Target ~ `Admission grade` + `Tuition fees up to date` + 
##     Gender + `Scholarship holder` + `Age at enrollment` + `Curricular units 1st sem (approved)` + 
##     `Curricular units 1st sem (grade)`, data = train)
## 
## Coefficients:
##          (Intercept) `Admission grade` `Tuition fees up to date`     Gender
## Enrolled   -2.560408       0.002398544                  2.043255 -0.1580618
## Graduate   -6.815370       0.020577113                  3.383690 -0.5481313
##          `Scholarship holder` `Age at enrollment`
## Enrolled            0.1508801         -0.04843702
## Graduate            1.1840470         -0.06145175
##          `Curricular units 1st sem (approved)`
## Enrolled                             0.1119614
## Graduate                             0.4402572
##          `Curricular units 1st sem (grade)`
## Enrolled                         0.07910346
## Graduate                         0.08398049
## 
## Std. Errors:
##          (Intercept) `Admission grade` `Tuition fees up to date`    Gender
## Enrolled   0.5662427       0.003826835                 0.2002750 0.1131168
## Graduate   0.5826152       0.003616929                 0.2511561 0.1095417
##          `Scholarship holder` `Age at enrollment`
## Enrolled            0.1612823         0.007701978
## Graduate            0.1386512         0.007157951
##          `Curricular units 1st sem (approved)`
## Enrolled                            0.03485786
## Graduate                            0.03233681
##          `Curricular units 1st sem (grade)`
## Enrolled                         0.01729699
## Graduate                         0.01938267
## 
## Residual Deviance: 5194.251 
## AIC: 5226.251
df_mn <- data.frame(
  Variabel = c("Admission grade",
               "Tuition fees up to date",
               "Gender",
               "Scholarship holder",
               "Age at enrollment",
               "Curricular units 1st sem (approved)",
               "Curricular units 1st sem (grade)"),
  
  Enrolled = c(0.002398544,
               2.043255,
               -0.1580618,
               0.1508801,
               -0.04843702,
               0.1119614,
               0.07910346),
  
  Graduate = c(0.020577113,
               3.383690,
               -0.5481313,
               1.1840470,
               -0.06145175,
               0.4402572,
               0.08398049)
)

kable(df_mn,
      caption = "Tabel 3.13 Koefisien Model Regresi Logistik Multinomial",
      align = c("c","c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.13 Koefisien Model Regresi Logistik Multinomial
Variabel Enrolled Graduate
Admission grade 0.0023985 0.0205771
Tuition fees up to date 2.0432550 3.3836900
Gender -0.1580618 -0.5481313
Scholarship holder 0.1508801 1.1840470
Age at enrollment -0.0484370 -0.0614518
Curricular units 1st sem (approved) 0.1119614 0.4402572
Curricular units 1st sem (grade) 0.0791035 0.0839805

3.2.6 Pengujian Hipotesis

3.2.6.1 Uji Simultan (Likelihood Ratio Test)

Uji signifikansi simultan dilakukan menggunakan Likelihood Ratio Test (LRT) untuk mengetahui apakah variabel prediktor secara bersama-sama berpengaruh terhadap variabel respon.

Hipotesis yang digunakan adalah:

  • \(H_0\) : \(\beta_1 = \beta_2 = \beta_3 = \cdots = \beta_k = 0\) (tidak terdapat pengaruh signifikan secara simultan antara variabel prediktor terhadap variabel respon)

  • \(H_1\) : paling sedikit terdapat satu \(\beta_j \neq 0\) (terdapat minimal satu variabel prediktor yang berpengaruh signifikan terhadap variabel respon)

# UJI SIMULTAN
model_null <- multinom(
  Target ~ 1,
  data = train
)
## # weights:  6 (2 variable)
## initial  value 3890.186114 
## final  value 3611.622899 
## converged
anova(model_null, model, test = "Chisq")
## Likelihood ratio tests of Multinomial Models
## 
## Response: Target
##                                                                                                                                                                              Model
## 1                                                                                                                                                                                1
## 2 `Admission grade` + `Tuition fees up to date` + Gender + `Scholarship holder` + `Age at enrollment` + `Curricular units 1st sem (approved)` + `Curricular units 1st sem (grade)`
##   Resid. df Resid. Dev   Test    Df LR stat. Pr(Chi)
## 1      7080   7223.246                              
## 2      7066   5194.251 1 vs 2    14 2028.995       0
df_simul_multi <- data.frame(
  Statistik_G = 2028.995,
  Df = 14,   
  P_value = "< 0.001",
  Keputusan = "Tolak H0"
)

kable(
  df_simul_multi,
  caption = "Tabel 3.14 Hasil Uji Simultan Model Regresi Logistik Multinomial",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1a3a2a"
  )
Tabel 3.14 Hasil Uji Simultan Model Regresi Logistik Multinomial
Statistik_G Df P_value Keputusan
2028.995 14 < 0.001 Tolak H0

Berdasarkan hasil pengujian, diperoleh nilai p-value < 0.05, sehingga \(H_0\) ditolak.

Artinya, variabel Admission grade, Tuition fees up to date, Gender, Scholarship holder, Age at enrollment, Curricular units approved, dan Curricular units grade secara simultan berpengaruh signifikan terhadap status mahasiswa (Dropout, Enrolled, Graduate).

Hal ini menunjukkan bahwa model multinomial yang dibangun lebih baik dibanding model tanpa prediktor.

3.2.6.2 Uji Parsial (Wald Test)

Uji signifikansi parsial dilakukan menggunakan uji Wald untuk mengetahui pengaruh masing-masing variabel prediktor terhadap variabel respon.

Hipotesis (untuk setiap variabel prediktor):

  • \(H_0 : \beta_i = 0\) (variabel ke-\(i\) tidak berpengaruh terhadap variabel respon)

  • \(H_1 : \beta_i \neq 0\) (variabel ke-\(i\) berpengaruh terhadap variabel respon)

z <- summary(model)$coefficients /
  summary(model)$standard.errors

p_value <- 2 * (1 - pnorm(abs(z)))
p_value
##           (Intercept) `Admission grade` `Tuition fees up to date`       Gender
## Enrolled 6.133051e-06      5.308102e-01                         0 1.623137e-01
## Graduate 0.000000e+00      1.277023e-08                         0 5.619368e-07
##          `Scholarship holder` `Age at enrollment`
## Enrolled            0.3495291          3.1971e-10
## Graduate            0.0000000          0.0000e+00
##          `Curricular units 1st sem (approved)`
## Enrolled                           0.001318411
## Graduate                           0.000000000
##          `Curricular units 1st sem (grade)`
## Enrolled                       4.802157e-06
## Graduate                       1.472509e-05
df_pval <- data.frame(
  Variabel = c("Intercept",
               "Admission grade",
               "Tuition fees up to date",
               "Gender",
               "Scholarship holder",
               "Age at enrollment",
               "Curricular units 1st sem (approved)",
               "Curricular units 1st sem (grade)"),
  
  Enrolled = c("<0.001",
               "0.531",
               "<0.001",
               "0.162",
               "0.350",
               "<0.001",
               "0.001",
               "<0.001"),
  
  Graduate = c("<0.001",
               "<0.001",
               "<0.001",
               "<0.001",
               "<0.001",
               "<0.001",
               "<0.001",
               "<0.001")
)

kable(df_pval,
      caption = "Tabel 3.15 Uji Signifikansi Parameter Model Regresi Logistik Multinomial",
      align = c("c","c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.15 Uji Signifikansi Parameter Model Regresi Logistik Multinomial
Variabel Enrolled Graduate
Intercept <0.001 <0.001
Admission grade 0.531 <0.001
Tuition fees up to date <0.001 <0.001
Gender 0.162 <0.001
Scholarship holder 0.350 <0.001
Age at enrollment <0.001 <0.001
Curricular units 1st sem (approved) 0.001 <0.001
Curricular units 1st sem (grade) <0.001 <0.001

Berdasarkan hasil uji signifikansi parameter pada model regresi logistik multinomial, diketahui bahwa pada kategori Enrolled, variabel Admission grade, Tuition fees up to date, Age at enrollment, Curricular units 1st semester (approved), dan Curricular units 1st semester (grade) memiliki nilai p-value < 0,05 sehingga berpengaruh signifikan terhadap peluang mahasiswa berada pada kategori tersebut dibandingkan kategori referensi (Dropout). Sementara itu, variabel Gender dan Scholarship holder tidak menunjukkan pengaruh yang signifikan pada kategori Enrolled.

Pada kategori Graduate, seluruh variabel yang digunakan dalam model menunjukkan nilai p-value < 0,05 sehingga seluruhnya berpengaruh signifikan terhadap peluang mahasiswa berada pada kategori Graduate dibandingkan Dropout. Hal ini menunjukkan bahwa variabel-variabel yang digunakan dalam model memiliki kontribusi yang lebih kuat dalam menjelaskan peluang kelulusan mahasiswa.

3.2.7 Odds Ratio

Odds Ratio (OR) digunakan untuk mengukur besar pengaruh variabel prediktor terhadap variabel respon dalam model regresi logistik. Nilai Odds Ratio diperoleh dengan melakukan eksponensial terhadap koefisien regresi.

Secara matematis, Odds Ratio dinyatakan sebagai:

\[ OR = e^{\beta_j} \]

dengan \(\beta_j\) adalah koefisien regresi dari variabel prediktor ke-j.

Pada regresi logistik multinomial, Odds Ratio menunjukkan perbandingan peluang suatu kategori terhadap kategori referensi (baseline).

Interpretasi Odds Ratio:

  • Jika \(OR > 1\), variabel prediktor meningkatkan peluang kejadian.
  • Jika \(OR < 1\), variabel prediktor menurunkan peluang kejadian.
  • Jika \(OR = 1\), variabel tidak berpengaruh terhadap peluang kejadian.
exp(coef(model))
##          (Intercept) `Admission grade` `Tuition fees up to date`   Gender
## Enrolled 0.077273216          1.002401                   7.71568 0.853797
## Graduate 0.001096787          1.020790                  29.47936 0.578029
##          `Scholarship holder` `Age at enrollment`
## Enrolled             1.162857           0.9527173
## Graduate             3.267571           0.9403983
##          `Curricular units 1st sem (approved)`
## Enrolled                              1.118470
## Graduate                              1.553107
##          `Curricular units 1st sem (grade)`
## Enrolled                           1.082316
## Graduate                           1.087608
df_or_multi <- data.frame(
  Variabel = c("Intercept",
               "Admission grade",
               "Tuition fees up to date",
               "Gender",
               "Scholarship holder",
               "Age at enrollment",
               "Curricular units 1st sem (approved)",
               "Curricular units 1st sem (grade)"),
  
  Enrolled = c(0.077273216,
               1.002401,
               7.715680,
               0.853797,
               1.162857,
               0.9527173,
               1.118470,
               1.082316),
  
  Graduate = c(0.001096787,
               1.020790,
               29.47936,
               0.578029,
               3.267571,
               0.9403983,
               1.553107,
               1.087608)
)

kable(df_or_multi,
      caption = "Tabel 3.16 Odds Ratio (Exp(β)) Model Regresi Logistik Multinomial",
      align = c("c","c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.16 Odds Ratio (Exp(β)) Model Regresi Logistik Multinomial
Variabel Enrolled Graduate
Intercept 0.0772732 0.0010968
Admission grade 1.0024010 1.0207900
Tuition fees up to date 7.7156800 29.4793600
Gender 0.8537970 0.5780290
Scholarship holder 1.1628570 3.2675710
Age at enrollment 0.9527173 0.9403983
Curricular units 1st sem (approved) 1.1184700 1.5531070
Curricular units 1st sem (grade) 1.0823160 1.0876080

Contoh Interpretasi:

  • Tuition fees up to date (Enrolled) (OR = 7,716): mahasiswa yang memiliki status pembayaran UKT tepat waktu memiliki peluang sekitar 7,7 kali lebih besar untuk berada pada kategori Enrolled dibandingkan Dropout, dengan variabel lain dianggap konstan.
  • Scholarship holder (Graduate) (OR = 3,268): mahasiswa yang menerima beasiswa memiliki peluang sekitar 3,27 kali lebih besar untuk berada pada kategori Graduate dibandingkan Dropout, dengan variabel lain dianggap konstan.

Secara keseluruhan, variabel-variabel akademik seperti Curricular units 1st semester (approved) dan Curricular units 1st semester (grade) menunjukkan kecenderungan meningkatkan peluang mahasiswa untuk berada pada kategori Graduate. Selain itu, variabel Scholarship holder juga memiliki pengaruh positif terhadap peluang kelulusan. Sementara itu, Age at enrollment menunjukkan kecenderungan pengaruh negatif terhadap peluang mahasiswa menjadi Graduate.

Dengan demikian, faktor-faktor yang berkaitan dengan performa akademik dan dukungan pendidikan cenderung berasosiasi dengan peningkatan peluang mahasiswa untuk menyelesaikan studi dibandingkan dengan mahasiswa yang tidak lulus (Dropout).

3.2.8 Pengujian Asumsi

3.2.8.1 Respons Bersifat Nominal

Variabel respon berupa status mahasiswa (Dropout, Enrolled, Graduate) yang merupakan kategori nominal tanpa urutan alami sehingga memenuhi asumsi regresi logistik multinomial.

3.2.8.2 Observasi Independen

Observasi diasumsikan saling independen karena setiap mahasiswa hanya tercatat satu kali dalam dataset.

3.2.8.3 Tidak Ada Multikolinearitas Berat

Uji multikolinearitas untuk mengetahui apakah terdapat hubungan yang kuat antarvariabel prediktor. Multikolinearitas yang tinggi dapat menyebabkan ketidakstabilan dalam estimasi parameter model sehingga interpretasi koefisien regresi menjadi kurang reliabel. Pada penelitian ini, uji multikolinearitas dilakukan menggunakan nilai Variance Inflation Factor (VIF).

Kriteria yang digunakan adalah:

  • Jika VIF < 10, maka tidak terdapat masalah multikolinearitas.

  • Jika VIF > 10, maka terdapat masalah multikolinearitas.

vif(
  lm(
    `Admission grade` ~
      `Tuition fees up to date` +
      Gender +
      `Scholarship holder` +
      `Age at enrollment` +
      `Curricular units 1st sem (approved)` +
      `Curricular units 1st sem (grade)`,
    data = train
  )
)
##             `Tuition fees up to date`                                Gender 
##                              1.107542                              1.067235 
##                  `Scholarship holder`                   `Age at enrollment` 
##                              1.086922                              1.092907 
## `Curricular units 1st sem (approved)`    `Curricular units 1st sem (grade)` 
##                              1.977469                              2.010186
df_vif <- data.frame(
  Variabel = c("Tuition fees up to date",
               "Gender",
               "Scholarship holder",
               "Age at enrollment",
               "Curricular units 1st sem (approved)",
               "Curricular units 1st sem (grade)"),
  
  VIF = c(1.107542,
          1.067235,
          1.086922,
          1.092907,
          1.977469,
          2.010186)
)

kable(df_vif,
      caption = "Tabel 3.17 Uji Multikolinearitas Regresi Logistik Multinomial",
      align = c("c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.17 Uji Multikolinearitas Regresi Logistik Multinomial
Variabel VIF
Tuition fees up to date 1.107542
Gender 1.067235
Scholarship holder 1.086922
Age at enrollment 1.092907
Curricular units 1st sem (approved) 1.977469
Curricular units 1st sem (grade) 2.010186

Nilai VIF seluruh variabel berada pada rentang 1,06 hingga 2,01. Karena seluruh nilai VIF < 10, maka dapat disimpulkan bahwa tidak terdapat masalah multikolinearitas antar variabel independen dalam model regresi logistik multinomial yang dibangun.

3.2.8.4 Tidak Ada Perfect Separation

Tidak ditemukan indikasi perfect separation karena model berhasil diestimasi tanpa peringatan dan tidak terdapat koefisien yang bernilai sangat besar.

3.2.8.5 Linearitas Logit

# DROPOUT VS GRADUATE
data_bg <- subset(data2, Target != "Enrolled")
data_bg$Target <- factor(data_bg$Target)

glm_bg <- glm(
  Target ~ `Admission grade` +
    `Age at enrollment` +
    `Curricular units 1st sem (approved)` +
    `Curricular units 1st sem (grade)` +
    I(`Admission grade`*log(`Admission grade`)) +
    I(`Age at enrollment`*log(`Age at enrollment`)) +
    I(`Curricular units 1st sem (approved)`*
        log(`Curricular units 1st sem (approved)`)) +
    I(`Curricular units 1st sem (grade)`*
        log(`Curricular units 1st sem (grade)`)),
  family = binomial,
  data = data_bg
)

summary(glm_bg)
## 
## Call:
## glm(formula = Target ~ `Admission grade` + `Age at enrollment` + 
##     `Curricular units 1st sem (approved)` + `Curricular units 1st sem (grade)` + 
##     I(`Admission grade` * log(`Admission grade`)) + I(`Age at enrollment` * 
##     log(`Age at enrollment`)) + I(`Curricular units 1st sem (approved)` * 
##     log(`Curricular units 1st sem (approved)`)) + I(`Curricular units 1st sem (grade)` * 
##     log(`Curricular units 1st sem (grade)`)), family = binomial, 
##     data = data_bg)
## 
## Coefficients:
##                                                                                       Estimate
## (Intercept)                                                                           -5.37041
## `Admission grade`                                                                     -0.09381
## `Age at enrollment`                                                                   -1.11218
## `Curricular units 1st sem (approved)`                                                  2.78506
## `Curricular units 1st sem (grade)`                                                     1.66865
## I(`Admission grade` * log(`Admission grade`))                                          0.01575
## I(`Age at enrollment` * log(`Age at enrollment`))                                      0.24097
## I(`Curricular units 1st sem (approved)` * log(`Curricular units 1st sem (approved)`)) -0.80262
## I(`Curricular units 1st sem (grade)` * log(`Curricular units 1st sem (grade)`))       -0.37536
##                                                                                       Std. Error
## (Intercept)                                                                              9.92633
## `Admission grade`                                                                        0.26001
## `Age at enrollment`                                                                      0.19331
## `Curricular units 1st sem (approved)`                                                    0.17887
## `Curricular units 1st sem (grade)`                                                       2.29122
## I(`Admission grade` * log(`Admission grade`))                                            0.04435
## I(`Age at enrollment` * log(`Age at enrollment`))                                        0.04382
## I(`Curricular units 1st sem (approved)` * log(`Curricular units 1st sem (approved)`))    0.05690
## I(`Curricular units 1st sem (grade)` * log(`Curricular units 1st sem (grade)`))          0.64483
##                                                                                       z value
## (Intercept)                                                                            -0.541
## `Admission grade`                                                                      -0.361
## `Age at enrollment`                                                                    -5.753
## `Curricular units 1st sem (approved)`                                                  15.570
## `Curricular units 1st sem (grade)`                                                      0.728
## I(`Admission grade` * log(`Admission grade`))                                           0.355
## I(`Age at enrollment` * log(`Age at enrollment`))                                       5.499
## I(`Curricular units 1st sem (approved)` * log(`Curricular units 1st sem (approved)`)) -14.105
## I(`Curricular units 1st sem (grade)` * log(`Curricular units 1st sem (grade)`))        -0.582
##                                                                                       Pr(>|z|)
## (Intercept)                                                                              0.588
## `Admission grade`                                                                        0.718
## `Age at enrollment`                                                                   8.74e-09
## `Curricular units 1st sem (approved)`                                                  < 2e-16
## `Curricular units 1st sem (grade)`                                                       0.466
## I(`Admission grade` * log(`Admission grade`))                                            0.722
## I(`Age at enrollment` * log(`Age at enrollment`))                                     3.82e-08
## I(`Curricular units 1st sem (approved)` * log(`Curricular units 1st sem (approved)`))  < 2e-16
## I(`Curricular units 1st sem (grade)` * log(`Curricular units 1st sem (grade)`))          0.560
##                                                                                          
## (Intercept)                                                                              
## `Admission grade`                                                                        
## `Age at enrollment`                                                                   ***
## `Curricular units 1st sem (approved)`                                                 ***
## `Curricular units 1st sem (grade)`                                                       
## I(`Admission grade` * log(`Admission grade`))                                            
## I(`Age at enrollment` * log(`Age at enrollment`))                                     ***
## I(`Curricular units 1st sem (approved)` * log(`Curricular units 1st sem (approved)`)) ***
## I(`Curricular units 1st sem (grade)` * log(`Curricular units 1st sem (grade)`))          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3566.9  on 2982  degrees of freedom
## Residual deviance: 2451.8  on 2974  degrees of freedom
##   (647 observations deleted due to missingness)
## AIC: 2469.8
## 
## Number of Fisher Scoring iterations: 5
df_bt <- data.frame(
  Variabel = c("Admission grade",
               "Age at enrollment",
               "Curricular units 1st sem (approved)",
               "Curricular units 1st sem (grade)"),
  
  `p-value (interaksi)` = c(0.722,
                            "<0.001",
                            "<0.001",
                            0.560),
  
  Kesimpulan = c("Terpenuhi",
                 "Tidak terpenuhi",
                 "Tidak terpenuhi",
                 "Terpenuhi")
)

kable(df_bt,
      caption = "Tabel 3.18 Uji Linearitas Logit (Box–Tidwell) Dropout vs Graduate",
      align = c("c","c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a") %>%
  row_spec(2, background = "#fdeaea") %>%
  row_spec(3, background = "#fdeaea")
Tabel 3.18 Uji Linearitas Logit (Box–Tidwell) Dropout vs Graduate
Variabel p.value..interaksi. Kesimpulan
Admission grade 0.722 Terpenuhi
Age at enrollment <0.001 Tidak terpenuhi
Curricular units 1st sem (approved) <0.001 Tidak terpenuhi
Curricular units 1st sem (grade) 0.56 Terpenuhi

Hasil uji Box–Tidwell menunjukkan bahwa variabel Age at enrollment dan Curricular units 1st semester (approved) memiliki p-value < 0,05 sehingga tidak memenuhi asumsi linearitas terhadap logit. Sementara itu, variabel Admission grade dan Curricular units 1st semester (grade) memiliki p-value > 0,05 sehingga memenuhi asumsi linearitas logit.

Dengan demikian, tidak semua variabel memenuhi asumsi linearitas, sehingga interpretasi model dilakukan dengan mempertimbangkan adanya potensi penyimpangan pada sebagian variabel.

data_be <- subset(data2, Target != "Graduate")
data_be$Target <- factor(data_be$Target)

glm_be <- glm(
  Target ~
    `Admission grade` +
    `Age at enrollment` +
    `Curricular units 1st sem (approved)` +
    `Curricular units 1st sem (grade)` +
    I(`Admission grade` * log(`Admission grade`)) +
    I(`Age at enrollment` * log(`Age at enrollment`)) +
    I(`Curricular units 1st sem (approved)` *
        log(`Curricular units 1st sem (approved)`)) +
    I(`Curricular units 1st sem (grade)` *
        log(`Curricular units 1st sem (grade)`)),
  family = binomial,
  data = data_be
)

summary(glm_be)
## 
## Call:
## glm(formula = Target ~ `Admission grade` + `Age at enrollment` + 
##     `Curricular units 1st sem (approved)` + `Curricular units 1st sem (grade)` + 
##     I(`Admission grade` * log(`Admission grade`)) + I(`Age at enrollment` * 
##     log(`Age at enrollment`)) + I(`Curricular units 1st sem (approved)` * 
##     log(`Curricular units 1st sem (approved)`)) + I(`Curricular units 1st sem (grade)` * 
##     log(`Curricular units 1st sem (grade)`)), family = binomial, 
##     data = data_be)
## 
## Coefficients:
##                                                                                        Estimate
## (Intercept)                                                                            1.215738
## `Admission grade`                                                                      0.023696
## `Age at enrollment`                                                                   -0.477299
## `Curricular units 1st sem (approved)`                                                  1.064119
## `Curricular units 1st sem (grade)`                                                     0.079137
## I(`Admission grade` * log(`Admission grade`))                                         -0.004418
## I(`Age at enrollment` * log(`Age at enrollment`))                                      0.097934
## I(`Curricular units 1st sem (approved)` * log(`Curricular units 1st sem (approved)`)) -0.360045
## I(`Curricular units 1st sem (grade)` * log(`Curricular units 1st sem (grade)`))       -0.030415
##                                                                                       Std. Error
## (Intercept)                                                                             9.645028
## `Admission grade`                                                                       0.260735
## `Age at enrollment`                                                                     0.207295
## `Curricular units 1st sem (approved)`                                                   0.178386
## `Curricular units 1st sem (grade)`                                                      2.128206
## I(`Admission grade` * log(`Admission grade`))                                           0.044544
## I(`Age at enrollment` * log(`Age at enrollment`))                                       0.047240
## I(`Curricular units 1st sem (approved)` * log(`Curricular units 1st sem (approved)`))   0.065964
## I(`Curricular units 1st sem (grade)` * log(`Curricular units 1st sem (grade)`))         0.602278
##                                                                                       z value
## (Intercept)                                                                             0.126
## `Admission grade`                                                                       0.091
## `Age at enrollment`                                                                    -2.303
## `Curricular units 1st sem (approved)`                                                   5.965
## `Curricular units 1st sem (grade)`                                                      0.037
## I(`Admission grade` * log(`Admission grade`))                                          -0.099
## I(`Age at enrollment` * log(`Age at enrollment`))                                       2.073
## I(`Curricular units 1st sem (approved)` * log(`Curricular units 1st sem (approved)`))  -5.458
## I(`Curricular units 1st sem (grade)` * log(`Curricular units 1st sem (grade)`))        -0.050
##                                                                                       Pr(>|z|)
## (Intercept)                                                                             0.8997
## `Admission grade`                                                                       0.9276
## `Age at enrollment`                                                                     0.0213
## `Curricular units 1st sem (approved)`                                                 2.44e-09
## `Curricular units 1st sem (grade)`                                                      0.9703
## I(`Admission grade` * log(`Admission grade`))                                           0.9210
## I(`Age at enrollment` * log(`Age at enrollment`))                                       0.0382
## I(`Curricular units 1st sem (approved)` * log(`Curricular units 1st sem (approved)`)) 4.81e-08
## I(`Curricular units 1st sem (grade)` * log(`Curricular units 1st sem (grade)`))         0.9597
##                                                                                          
## (Intercept)                                                                              
## `Admission grade`                                                                        
## `Age at enrollment`                                                                   *  
## `Curricular units 1st sem (approved)`                                                 ***
## `Curricular units 1st sem (grade)`                                                       
## I(`Admission grade` * log(`Admission grade`))                                            
## I(`Age at enrollment` * log(`Age at enrollment`))                                     *  
## I(`Curricular units 1st sem (approved)` * log(`Curricular units 1st sem (approved)`)) ***
## I(`Curricular units 1st sem (grade)` * log(`Curricular units 1st sem (grade)`))          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2171.6  on 1573  degrees of freedom
## Residual deviance: 2047.8  on 1565  degrees of freedom
##   (641 observations deleted due to missingness)
## AIC: 2065.8
## 
## Number of Fisher Scoring iterations: 4
df_bt2 <- data.frame(
  Variabel = c("Admission grade",
               "Age at enrollment",
               "Curricular units 1st sem (approved)",
               "Curricular units 1st sem (grade)"),
  
  `p-value (interaksi)` = c(0.921,
                            0.038,
                            "<0.001",
                            0.960),
  
  Kesimpulan = c("Terpenuhi",
                 "Tidak terpenuhi",
                 "Tidak terpenuhi",
                 "Terpenuhi")
)

kable(df_bt2,
      caption = "Tabel 3.19 Uji Linearitas Logit (Box–Tidwell) Dropout vs Enrolled",
      align = c("c","c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a") %>%
  row_spec(2, background = "#fdeaea") %>%
  row_spec(3, background = "#fdeaea")
Tabel 3.19 Uji Linearitas Logit (Box–Tidwell) Dropout vs Enrolled
Variabel p.value..interaksi. Kesimpulan
Admission grade 0.921 Terpenuhi
Age at enrollment 0.038 Tidak terpenuhi
Curricular units 1st sem (approved) <0.001 Tidak terpenuhi
Curricular units 1st sem (grade) 0.96 Terpenuhi

Hasil uji Box–Tidwell pada logit Dropout vs Enrolled menunjukkan bahwa variabel Age at enrollment dan Curricular units 1st semester (approved) tidak memenuhi asumsi linearitas logit karena memiliki p-value < 0,05. Sementara itu, variabel Admission grade dan Curricular units 1st semester (grade) memenuhi asumsi linearitas logit dengan p-value > 0,05. Hal ini menunjukkan bahwa terdapat pelanggaran asumsi pada sebagian variabel dalam model, sehingga interpretasi model perlu dilakukan dengan mempertimbangkan kondisi tersebut.

3.2.8.6 Independence of Irrelevant Alternatives (IIA)

Uji Independence of Irrelevant Alternatives (IIA) dilakukan untuk menguji apakah penambahan atau pengurangan alternatif pilihan tidak memengaruhi rasio peluang antar alternatif lain dalam model multinomial logit.

Hipotesis:

  • \(H_0:\): Asumsi IIA terpenuhi (tidak ada pengaruh dari alternatif lain terhadap perbandingan dua alternatif).
  • \(H_1:\): Asumsi IIA tidak terpenuhi (terdapat pengaruh alternatif lain terhadap perbandingan dua alternatif).
data_mlogit <- mlogit.data(
  data2,
  choice = "Target",
  shape = "wide"
)

model_mlogit <- mlogit(
  Target ~ 1 |
    Admission.grade +
    Age.at.enrollment +
    Gender +
    Scholarship.holder +
    Curricular.units.1st.sem..approved. +
    Curricular.units.1st.sem..grade.,
  data = data_mlogit
)

model_restricted1 <- update(
  model_mlogit,
  alt.subset = c("Dropout", "Graduate")
)

hmftest(model_mlogit, model_restricted1)
## 
##  Hausman-McFadden test
## 
## data:  data_mlogit
## chisq = -156.55, df = 7, p-value = 1
## alternative hypothesis: IIA is rejected
model_restricted2 <- update(
  model_mlogit,
  alt.subset = c("Dropout", "Enrolled")
)

hmftest(model_mlogit, model_restricted2)
## 
##  Hausman-McFadden test
## 
## data:  data_mlogit
## chisq = 23.579, df = 7, p-value = 0.001351
## alternative hypothesis: IIA is rejected
df_iia <- data.frame(
  Perbandingan = c(
    "Full vs Restricted (Dropout vs Enrolled)",
    "Full vs Restricted (subset)"
  ),
  
  `Chi-square` = c(-156.55, 23.579),
  df = c(7, 7),
  `p-value` = c(1, 0.001351),
  
  Kesimpulan = c(
    "Tidak ada bukti pelanggaran IIA",
    "Asumsi IIA tidak terpenuhi"
  )
)

kable(df_iia,
      caption = "Tabel 3.20 Uji Asumsi IIA (Hausman–McFadden Test)",
      align = c("c","c","c","c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a") %>%
  row_spec(2, background = "#fdeaea")
Tabel 3.20 Uji Asumsi IIA (Hausman–McFadden Test)
Perbandingan Chi.square df p.value Kesimpulan
Full vs Restricted (Dropout vs Enrolled) -156.550 7 1.000000 Tidak ada bukti pelanggaran IIA
Full vs Restricted (subset) 23.579 7 0.001351 Asumsi IIA tidak terpenuhi

Berdasarkan Hausman-McFadden Test, diperoleh hasil yang berbeda pada beberapa subset alternatif. Pada subset Dropout–Graduate, tidak terdapat bukti pelanggaran IIA (p-value > 0,05), sedangkan pada subset Dropout–Enrolled terdapat indikasi pelanggaran IIA (p-value < 0,05). Oleh karena itu, asumsi IIA tidak sepenuhnya terpenuhi pada model multinomial logit yang digunakan.

3.2.9 Evaluasi Model

3.2.9.1 Koefisien Determinasi

Koefisien determinasi digunakan untuk mengevaluasi kemampuan model regresi logistik dalam menjelaskan variabilitas variabel respon.

# PSEUDO R2
pR2(model)
## fitting null model for pseudo-r2
## # weights:  6 (2 variable)
## initial  value 3890.186114 
## final  value 3611.622899 
## converged
##           llh       llhNull            G2      McFadden          r2ML 
## -2597.1254685 -3611.6228993  2028.9948615     0.2808979     0.4361689 
##          r2CU 
##     0.5013684
df_r2_simple <- data.frame(
  Indikator = c("McFadden R²",
                "Cox & Snell (ML)",
                "Nagelkerke (Cragg & Uhler)"),
  
  Nilai = c(0.2808979,
            0.4361689,
            0.5013684)
)

kable(df_r2_simple,
      caption = "Tabel 3.21 Nilai Pseudo R-Square Model Regresi Logistik Multinomial",
      align = c("c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.21 Nilai Pseudo R-Square Model Regresi Logistik Multinomial
Indikator Nilai
McFadden R² 0.2808979
Cox & Snell (ML) 0.4361689
Nagelkerke (Cragg & Uhler) 0.5013684

Nilai pseudo R-square menunjukkan bahwa model memiliki kemampuan penjelasan yang cukup baik, dengan McFadden R² sebesar 0,2809, Cox & Snell R² sebesar 0,4362, dan Nagelkerke R² sebesar 0,5014 yang mengindikasikan bahwa model mampu menjelaskan sekitar 50% variasi pada data secara relatif baik.

3.2.9.2 Confusion Matrix dan Akurasi

Kinerja klasifikasi model regresi logistik dievaluasi menggunakan confusion matrix dan nilai akurasi. Confusion matrix digunakan untuk membandingkan hasil prediksi model dengan kondisi aktual sehingga dapat diketahui jumlah observasi yang diklasifikasikan dengan benar maupun salah.

prediksi <- predict(model, newdata = test)

table(Aktual = test$Target, Prediksi = prediksi)
##           Prediksi
## Aktual     Dropout Enrolled Graduate
##   Dropout      223        2       59
##   Enrolled      54        6       98
##   Graduate      41        3      397
df_cm_multi <- data.frame(
  Aktual = c("Dropout", "Enrolled", "Graduate"),
  Dropout = c(223, 54, 41),
  Enrolled = c(2, 6, 3),
  Graduate = c(59, 98, 397)
)

kable(df_cm_multi,
      caption = "Tabel 3.22 Confusion Matrix Model Regresi Logistik Multinomial",
      align = c("c","c","c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.22 Confusion Matrix Model Regresi Logistik Multinomial
Aktual Dropout Enrolled Graduate
Dropout 223 2 59
Enrolled 54 6 98
Graduate 41 3 397

Berdasarkan hasil confusion matrix, model mampu mengklasifikasikan sebagian besar data dengan benar, terutama pada kategori Dropout dan Graduate. Namun, terdapat kesalahan klasifikasi yang cukup terlihat pada kategori Enrolled, yang cenderung diprediksi ke kategori lain seperti Graduate atau Dropout. Hal ini menunjukkan bahwa performa model dalam membedakan kategori Enrolled masih relatif lebih rendah dibandingkan dua kategori lainnya.

confusionMatrix(prediksi, test$Target)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Dropout Enrolled Graduate
##   Dropout      223       54       41
##   Enrolled       2        6        3
##   Graduate      59       98      397
## 
## Overall Statistics
##                                           
##                Accuracy : 0.7089          
##                  95% CI : (0.6778, 0.7387)
##     No Information Rate : 0.4994          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.4881          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
## 
## Statistics by Class:
## 
##                      Class: Dropout Class: Enrolled Class: Graduate
## Sensitivity                  0.7852        0.037975          0.9002
## Specificity                  0.8414        0.993103          0.6448
## Pos Pred Value               0.7013        0.545455          0.7166
## Neg Pred Value               0.8920        0.825688          0.8663
## Prevalence                   0.3216        0.178935          0.4994
## Detection Rate               0.2525        0.006795          0.4496
## Detection Prevalence         0.3601        0.012458          0.6274
## Balanced Accuracy            0.8133        0.515539          0.7725
df_eval_multi <- data.frame(
  Metrik = c("Accuracy",
             "95% CI Lower",
             "95% CI Upper",
             "No Information Rate",
             "Kappa",
             "P-Value (Acc > NIR)"),
  
  Nilai = c(0.7089,
            0.6778,
            0.7387,
            0.4994,
            0.4881,
            "<0.001")
)

kable(df_eval_multi,
      caption = "Tabel 3.23 Evaluasi Kinerja Model Regresi Logistik Multinomial",
      align = c("c","c"),
      row.names = FALSE) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.23 Evaluasi Kinerja Model Regresi Logistik Multinomial
Metrik Nilai
Accuracy 0.7089
95% CI Lower 0.6778
95% CI Upper 0.7387
No Information Rate 0.4994
Kappa 0.4881
P-Value (Acc > NIR) <0.001

Model regresi logistik multinomial menghasilkan nilai accuracy sebesar 0,7089 dengan Kappa sebesar 0,4881, yang menunjukkan tingkat kesesuaian prediksi yang moderat. Nilai accuracy yang lebih tinggi dibandingkan No Information Rate (0,4994) serta p-value < 0,001 menunjukkan bahwa model secara signifikan lebih baik dibandingkan model tanpa informasi. Secara keseluruhan, model memiliki performa klasifikasi yang cukup baik, meskipun masih terdapat ketidakseimbangan performa pada masing-masing kelas, terutama pada kategori Enrolled yang memiliki sensitivitas rendah.

3.2.10 Kesimpulan

Berdasarkan hasil analisis regresi logistik multinomial terhadap status mahasiswa (Dropout, Enrolled, dan Graduate), diperoleh bahwa model yang dibangun memiliki performa klasifikasi yang cukup baik dengan nilai akurasi sebesar 0,7089 dan nilai Kappa sebesar 0,4881, yang menunjukkan tingkat kesesuaian prediksi pada kategori sedang hingga cukup baik.

Hasil uji signifikansi parameter menunjukkan bahwa beberapa variabel seperti Age at enrollment, Curricular units 1st semester (approved), serta Curricular units 1st semester (grade) berpengaruh signifikan terhadap status mahasiswa, sementara variabel lain seperti Gender dan sebagian variabel akademik tertentu tidak selalu signifikan pada semua logit.

Dari hasil Odds Ratio, variabel yang berkaitan dengan performa akademik seperti jumlah unit semester pertama yang disetujui dan nilai semester pertama menunjukkan kecenderungan peningkatan peluang mahasiswa untuk berada pada kategori Graduate, dibandingkan dengan Dropout.

Namun, hasil uji asumsi menunjukkan bahwa tidak semua asumsi terpenuhi. Uji linearitas logit (Box–Tidwell) menunjukkan adanya pelanggaran pada beberapa variabel, sementara uji independensi alternatif (IIA) menunjukkan bahwa asumsi IIA tidak terpenuhi. Selain itu, uji multikolinearitas menunjukkan tidak terdapat masalah korelasi tinggi antar variabel independen.

Secara keseluruhan, model regresi logistik multinomial yang dibangun mampu memberikan gambaran yang cukup baik dalam mengklasifikasikan status mahasiswa, meskipun terdapat beberapa keterbatasan pada pemenuhan asumsi model yang perlu diperhatikan dalam interpretasi hasil.

3.3 Regresi Logistik Ordinal

3.3.1 Deskripsi Kasus

Kualitas white wine merupakan salah satu indikator penting dalam industri minuman fermentasi karena berpengaruh terhadap preferensi konsumen dan nilai jual produk. Penilaian kualitas wine biasanya didasarkan pada karakteristik sensorik yang dapat dipengaruhi oleh berbagai komponen kimia di dalamnya, seperti kadar alkohol, keasaman, kandungan gula, dan senyawa lainnya.

Analisis ini bertujuan untuk mengidentifikasi faktor-faktor kimia yang memengaruhi tingkat kualitas wine serta menjelaskan bagaimana variabel-variabel tersebut berkontribusi terhadap peningkatan atau penurunan kualitas. Oleh karena itu, digunakan metode regresi logistik ordinal untuk mengevaluasi pengaruh variabel-variabel kimia terhadap peluang suatu wine berada pada kategori kualitas rendah, sedang, maupun tinggi, sekaligus menilai kemampuan model dalam mengklasifikasikan kualitas wine berdasarkan karakteristik tersebut.

Dataset yang digunakan dalam analisis ini adalah White Wine Quality Dataset yang berasal dari UCI Machine Learning Repository. Dataset ini berisi hasil pengukuran kimia dari berbagai sampel white wine yang diproduksi di Portugal, dengan total 4.898 observasi.

Setiap sampel wine memiliki karakteristik kimia yang digunakan sebagai variabel penjelas (independent variables). Sementara itu, variabel respon adalah quality, yaitu skor kualitas wine yang diberikan oleh panel sensorik dengan rentang nilai 0 hingga 10.

Dalam penelitian ini, variabel kualitas tersebut dikategorikan menjadi tiga tingkat ordinal, yaitu:

  • Low (rendah): kualitas rendah (≤ 5)
  • Medium (sedang): kualitas sedang (= 6)
  • High (tinggi): kualitas baik (≥ 7)

Pengelompokan ini dilakukan untuk memungkinkan pemodelan menggunakan regresi logistik ordinal, karena variabel respon memiliki struktur berurutan (ordered categories).

Variabel yang digunakan ditampilkan pada tabel berikut:

Tabel 3.24 Variabel Regresi Logistik Ordinal

Variabel Tipe Keterangan
quality Numerik (Diskrit) Skor kualitas wine (0–10) berdasarkan penilaian sensorik
alcohol Numerik Kadar alkohol dalam wine
density Numerik Massa jenis wine
volatile acidity Numerik Tingkat keasaman volatil
residual sugar Numerik Sisa gula dalam wine
sulphates Numerik Kandungan sulfat dalam wine

3.3.2 Preparasi Data

Tahap awal dalam analisis ini adalah preparasi data yang bertujuan untuk memastikan dataset siap digunakan dalam analisis.

# LOAD PACKAGES
library(readr)
library(dplyr)
library(MASS)
library(pscl)
library(readxl)
library(dplyr)
library(ggplot2)
library(tidyr)
library(psych)
library(caret)
library(brant)
## Warning: package 'brant' was built under R version 4.4.3
# IMPORT DATA
wine <- read_excel(
  "C:/Users/Asus/Downloads/DATA ADK.xlsx",
  sheet = "Ordinal"
)

wine <- wine %>%
  dplyr::select(
    quality,
    alcohol,
    density,
    `volatile acidity`,
    `residual sugar`,
    sulphates
  )

# KLASIFIKASI Y
wine <- wine %>%
  mutate(
    quality_ord = case_when(
      quality <= 5 ~ "low",
      quality == 6 ~ "medium",
      quality >= 7 ~ "high"
    )
  )

wine$quality_ord <- factor(
  wine$quality_ord,
  levels = c("low", "medium", "high"),
  ordered = TRUE
)

3.3.3 Eksplorasi Data

Eksplorasi data dilakukan untuk memberikan gambaran awal mengenai karakteristik data serta hubungan antara variabel penjelas dengan variabel respon.

# STATISTIK DESKRIPTIF 
stat_desc3 <- wine %>%
  summarise(
    n = n(),
    Mean_Alcohol = mean(alcohol, na.rm = TRUE),
    Mean_Density = mean(density, na.rm = TRUE),
    Mean_VA = mean(`volatile acidity`, na.rm = TRUE),
    Mean_RS = mean(`residual sugar`, na.rm = TRUE),
    Mean_Sulphates = mean(sulphates, na.rm = TRUE),
    Mean_Quality = mean(quality, na.rm = TRUE)
  )
# DISTRIBUSI KELAS
ggplot(wine, aes(x = quality_ord, fill = quality_ord)) +
  geom_bar() +
  geom_text(
    stat = "count",
    aes(label = after_stat(count)),
    vjust = -0.5
  ) +
  theme_minimal() +
  labs(
    title = "Distribusi Kualitas White Wine (Ordinal)",
    x = "Kualitas",
    y = "Frekuensi"
  )

Grafik distribusi menunjukkan bahwa kategori kualitas medium memiliki jumlah observasi terbanyak, yaitu 2198 sampel (44.9%), diikuti kategori low sebanyak 1640 sampel (33.5%), dan kategori high sebanyak 1060 sampel (21.6%).

Hal ini menunjukkan bahwa sebagian besar white wine dalam dataset memiliki kualitas sedang. Meskipun terdapat perbedaan jumlah antar kategori, distribusi data masih cukup representatif untuk dilakukan analisis regresi ordinal karena seluruh kategori memiliki jumlah observasi yang memadai.

# BOXPLOT
wine_long <- wine %>%
  pivot_longer(
    cols = c(alcohol, density, `volatile acidity`,
             `residual sugar`, sulphates),
    names_to = "variabel",
    values_to = "nilai"
  )

ggplot(
  wine_long,
  aes(x = quality_ord, y = nilai, fill = quality_ord)
) +
  geom_boxplot() +
  facet_wrap(~variabel, scales = "free") +
  theme_minimal() +
  labs(title = "Boxplot Variabel Prediktor")

Berdasarkan boxplot, variabel yang menunjukkan perbedaan paling jelas antar kategori kualitas wine adalah alcohol dan volatile acidity. Kualitas wine cenderung meningkat pada kadar alcohol yang lebih tinggi dan volatile acidity yang lebih rendah. Sementara itu, density dan sulphates juga menunjukkan kecenderungan hubungan dengan kualitas, sedangkan residual sugar memperlihatkan perbedaan yang relatif lemah.

3.3.4 Pemodelan Regresi Logistik Ordinal

Regresi logistik ordinal digunakan untuk memodelkan hubungan antara variabel respon yang bersifat kategorik berurutan (ordered), yaitu terdiri atas lebih dari dua kategori dengan adanya tingkatan atau urutan. Pada penelitian ini, variabel respon terdiri atas tiga kategori, yaitu Dropout, Enrolled, dan Graduate, dengan asumsi terdapat urutan dari kondisi terendah hingga tertinggi.

Secara umum, model regresi logistik ordinal (cumulative logit model) dapat dituliskan sebagai berikut:

\[ \log\left(\frac{P(Y \le j)}{P(Y > j)}\right) = \alpha_j - \left(\beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p\right) \]

dengan:

  • \(Y\) = variabel respon (status mahasiswa)
  • \(j\) = titik batas kategori (cutpoint)
  • \(X_i\) = variabel prediktor ke-i
  • \(\beta_i\) = koefisien regresi untuk masing-masing prediktor
  • \(\alpha_j\) = parameter ambang (threshold/intercept) untuk kategori ke-j

Dalam model ini digunakan pendekatan cumulative logit, yaitu membandingkan peluang kumulatif \((P(Y \le j)\) terhadap \(P(Y > j)\). Dengan demikian, kategori Dropout, Enrolled, dan Graduate diperlakukan sebagai variabel ordinal yang memiliki tingkatan berurutan.

# PEMODELAN REGRESI LOGISTIK ORDINAL
model <- polr(
  quality_ord ~ alcohol + density + `volatile acidity` +
    `residual sugar` + sulphates,
  data = wine,
  Hess = TRUE
)

summary(model)
## Call:
## polr(formula = quality_ord ~ alcohol + density + `volatile acidity` + 
##     `residual sugar` + sulphates, data = wine, Hess = TRUE)
## 
## Coefficients:
##                        Value Std. Error   t value
## alcohol               0.7479   0.030016    24.918
## density            -211.5570   0.177692 -1190.586
## `volatile acidity`   -5.4606   0.323756   -16.866
## `residual sugar`      0.1405   0.006567    21.401
## sulphates             1.8171   0.251157     7.235
## 
## Intercepts:
##             Value      Std. Error t value   
## low|medium   -203.0532     0.1773 -1145.4639
## medium|high  -200.4898     0.1989 -1007.8503
## 
## Residual Deviance: 8739.417 
## AIC: 8753.417
coef_df <- data.frame(
  Variabel = c("Alcohol", "Density", "Volatile acidity", "Residual sugar", "Sulphates"),
  Estimate = c(0.7479, -211.5570, -5.4606, 0.1405, 1.8171)
)

kable(coef_df,
      caption = "Tabel 3.25 Hasil Estimasi Model Regresi Logistik Ordinal",
      align = "c") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.25 Hasil Estimasi Model Regresi Logistik Ordinal
Variabel Estimate
Alcohol 0.7479
Density -211.5570
Volatile acidity -5.4606
Residual sugar 0.1405
Sulphates 1.8171

3.3.5 Pengujian Hipotesis

3.3.5.1 Uji Simultan (Likelihood Ratio Test)

Uji signifikansi simultan dilakukan menggunakan Likelihood Ratio Test (LRT) untuk mengetahui apakah variabel prediktor secara bersama-sama berpengaruh terhadap variabel respon.

Hipotesis yang digunakan adalah:

  • \(H_0\) : \(\beta_1 = \beta_2 = \beta_3 = \cdots = \beta_k = 0\) (tidak terdapat pengaruh signifikan secara simultan antara variabel prediktor terhadap variabel respon)

  • \(H_1\) : paling sedikit terdapat satu \(\beta_j \neq 0\) (terdapat minimal satu variabel prediktor yang berpengaruh signifikan terhadap variabel respon)

# UJI SIMULTAN
model_null <- polr(
  quality_ord ~ 1,
  data = wine,
  Hess = TRUE
)

anova(model_null, model, test = "Chisq")
## Likelihood ratio tests of ordinal regression models
## 
## Response: quality_ord
##                                                                   Model
## 1                                                                     1
## 2 alcohol + density + `volatile acidity` + `residual sugar` + sulphates
##   Resid. df Resid. Dev   Test    Df LR stat. Pr(Chi)
## 1      4896  10355.955                              
## 2      4891   8739.417 1 vs 2     5 1616.538       0
lrt_df <- data.frame(
  Statistik_G = 1616.538,
  Df = 5,
  P_value = "< 0.001",
  Keputusan = "Tolak H0"
)

kable(
  lrt_df,
  caption = "Tabel 3.26 Hasil Uji Simultan Model Regresi Logistik Ordinal",
  align = "c"
) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1a3a2a"
  )
Tabel 3.26 Hasil Uji Simultan Model Regresi Logistik Ordinal
Statistik_G Df P_value Keputusan
1616.538 5 < 0.001 Tolak H0

Hasil uji likelihood ratio menunjukkan bahwa nilai statistik uji sebesar 1616.538 dengan derajat bebas 5 menghasilkan p-value yang sangat kecil (p < 0.001). Hal ini mengindikasikan bahwa model dengan variabel prediktor secara signifikan lebih baik dibandingkan model tanpa prediktor.

Dengan demikian, dapat disimpulkan bahwa secara simultan variabel alcohol, density, volatile acidity, residual sugar, dan sulphates berpengaruh signifikan terhadap kategori kualitas wine.

3.3.5.2 Uji Parsial (Wald Test)

Uji signifikansi parsial dilakukan menggunakan uji Wald untuk mengetahui pengaruh masing-masing variabel prediktor terhadap variabel respon.

Hipotesis (untuk setiap variabel prediktor):

  • \(H_0 : \beta_i = 0\) (variabel ke-\(i\) tidak berpengaruh terhadap variabel respon)

  • \(H_1 : \beta_i \neq 0\) (variabel ke-\(i\) berpengaruh terhadap variabel respon)

# UJI SIGNIFIKANSI
ctable <- coef(summary(model))
p_value <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable <- cbind(ctable, "p value" = p_value)
ctable
##                           Value  Std. Error     t value       p value
## alcohol               0.7479495 0.030016275    24.91813 4.732806e-137
## density            -211.5570480 0.177691549 -1190.58587  0.000000e+00
## `volatile acidity`   -5.4605934 0.323756278   -16.86637  7.953556e-64
## `residual sugar`      0.1405422 0.006566996    21.40129 1.299497e-101
## sulphates             1.8171017 0.251157466     7.23491  4.658373e-13
## low|medium         -203.0532477 0.177267254 -1145.46395  0.000000e+00
## medium|high        -200.4898156 0.198928175 -1007.85027  0.000000e+00
ctable_df <- data.frame(
  Variabel = c("Alcohol", "Density", "Volatile acidity", "Residual sugar", "Sulphates",
               "low | medium", "medium | high"),
  Estimate = c(0.7479495, -211.5570480, -5.4605934, 0.1405422, 1.8171017,
               -203.0532477, -200.4898156),
  Std_Error = c(0.030016275, 0.177691549, 0.323756278, 0.006566996, 0.251157466,
                0.177267254, 0.198928175),
  t_value = c(24.91813, -1190.58587, -16.86637, 21.40129, 7.23491,
              -1145.46395, -1007.85027),
  p_value = c(4.732806e-137, 0, 7.953556e-64, 1.299497e-101, 4.658373e-13,
              0, 0)
)

ctable_df$Signif <- ifelse(ctable_df$p_value < 0.05, "Signifikan", "Tidak signifikan")

kable(ctable_df,
      caption = "Tabel 3.27 Hasil Uji Parsial Model Regresi Logistik Ordinal",
      align = "c") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.27 Hasil Uji Parsial Model Regresi Logistik Ordinal
Variabel Estimate Std_Error t_value p_value Signif
Alcohol 0.7479495 0.0300163 24.91813 0 Signifikan
Density -211.5570480 0.1776915 -1190.58587 0 Signifikan
Volatile acidity -5.4605934 0.3237563 -16.86637 0 Signifikan
Residual sugar 0.1405422 0.0065670 21.40129 0 Signifikan
Sulphates 1.8171017 0.2511575 7.23491 0 Signifikan
low &#124; medium -203.0532477 0.1772673 -1145.46395 0 Signifikan
medium &#124; high -200.4898156 0.1989282 -1007.85027 0 Signifikan

Secara parsial, seluruh variabel prediktor yaitu alcohol, density, volatile acidity, residual sugar, dan sulphates berpengaruh signifikan terhadap kualitas wine, yang ditunjukkan dengan nilai p-value < 0.05 berdasarkan uji Wald untuk seluruh variabel.

3.3.6 Odds Ratio

Pada regresi logistik ordinal, Odds Ratio (OR) digunakan untuk mengukur perubahan odds suatu variabel prediktor terhadap peluang berpindah ke kategori yang lebih tinggi dari variabel respon.

Odds Ratio diperoleh dari eksponensial koefisien regresi:

\[ OR = e^{\beta_j} \]

dengan \(\beta_j\) adalah koefisien regresi dari variabel prediktor ke-j.

Interpretasi Odds Ratio:

  • Jika \(OR > 1\), maka peningkatan nilai variabel prediktor akan meningkatkan peluang berada pada kategori yang lebih tinggi dari variabel respon.
  • Jika \(OR < 1\), maka peningkatan nilai variabel prediktor akan menurunkan peluang berada pada kategori yang lebih tinggi.
  • Jika \(OR = 1\), maka variabel tidak memiliki pengaruh terhadap peluang kategori respon.
# ODDS RATIO
exp(coef(model))
##            alcohol            density `volatile acidity`   `residual sugar` 
##       2.112664e+00       1.324163e-92       4.251033e-03       1.150898e+00 
##          sulphates 
##       6.153997e+00
or_df <- data.frame(
  Variabel = c("Alcohol", "Density", "Volatile acidity", "Residual sugar", "Sulphates"),
  Odds_Ratio = c(2.112664, 1.324163e-92, 4.251033e-03, 1.150898, 6.153997)
)

kable(or_df,
      caption = "Tabel 3.28 Odds Ratio Model Regresi Logistik Ordinal",
      align = "c") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.28 Odds Ratio Model Regresi Logistik Ordinal
Variabel Odds_Ratio
Alcohol 2.112664
Density 0.000000
Volatile acidity 0.004251
Residual sugar 1.150898
Sulphates 6.153997

Interpretasi:

  • Alcohol (OR = 2.11) → setiap kenaikan 1 unit alcohol meningkatkan peluang masuk ke kategori kualitas yang lebih tinggi sebesar 2.11 kali, dengan asumsi variabel lain konstan.
  • Density (OR ≈ 1.32e-92) → nilai OR sangat kecil (~0), menunjukkan bahwa peningkatan density secara kuat menurunkan peluang kualitas lebih tinggi.
  • Volatile acidity (OR = 0.00425) → meningkatkan volatile acidity akan menurunkan peluang masuk kategori kualitas yang lebih tinggi secara signifikan.
  • Residual sugar (OR = 1.15) → setiap kenaikan residual sugar meningkatkan peluang kualitas lebih tinggi sebesar 1.15 kali, meskipun efeknya kecil.
  • Sulphates (OR = 6.15) → variabel ini memiliki pengaruh paling kuat positif, yaitu meningkatkan peluang kualitas lebih tinggi hingga 6.15 kali.

3.3.7 Pengujian Asumsi Proportional Odds

Uji Brant digunakan untuk menguji asumsi proportional odds pada model regresi logistik ordinal.

Hipotesis:

  • \(H_0:\) Model memenuhi asumsi proportional odds (tidak ada perbedaan efek variabel antar kategori).
  • \(H_1:\) Model tidak memenuhi asumsi proportional odds (terdapat perbedaan efek variabel antar kategori).
wine <- wine %>%
  rename(
    volatile_acidity = `volatile acidity`,
    residual_sugar = `residual sugar`
  )

model_brant <- polr(
  quality_ord ~ alcohol + density + volatile_acidity +
    residual_sugar + sulphates,
  data = wine,
  Hess = TRUE
)

brant(model_brant)
## ---------------------------------------------------- 
## Test for     X2  df  probability 
## ---------------------------------------------------- 
## Omnibus          27.12   5   0
## alcohol          2.22    1   0.14
## density          0.01    1   0.94
## volatile_acidity 25.78   1   0
## residual_sugar       0.75    1   0.39
## sulphates        0.13    1   0.72
## ---------------------------------------------------- 
## 
## H0: Parallel Regression Assumption holds
brant_df <- data.frame(
  Variabel = c("Omnibus", "Alcohol", "Density", "Volatile acidity", "Residual sugar", "Sulphates"),
  Chi_Square = c(27.12, 2.22, 0.01, 25.78, 0.75, 0.13),
  df = c(5, 1, 1, 1, 1, 1),
  p_value = c(0, 0.14, 0.94, 0, 0.39, 0.72)
)

kable(brant_df,
      caption = "Tabel 3.29 Hasil Uji Brant",
      align = "c") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a") %>%
  row_spec(which(brant_df$p_value < 0.05), background = "#ffcccc")
Tabel 3.29 Hasil Uji Brant
Variabel Chi_Square df p_value
Omnibus 27.12 5 0.00
Alcohol 2.22 1 0.14
Density 0.01 1 0.94
Volatile acidity 25.78 1 0.00
Residual sugar 0.75 1 0.39
Sulphates 0.13 1 0.72

Hasil uji Brant menunjukkan bahwa secara omnibus diperoleh p-value < 0,05, sehingga secara keseluruhan terdapat indikasi pelanggaran asumsi parallel regression. Namun, secara parsial hanya variabel volatile acidity yang menunjukkan pelanggaran asumsi (p < 0,05), sedangkan variabel lainnya seperti alcohol, density, residual sugar, dan sulphates memenuhi asumsi.

Dengan demikian, dapat disimpulkan bahwa sebagian besar variabel masih memenuhi asumsi proportional odds, meskipun terdapat satu variabel yang menunjukkan pelanggaran.

3.3.8 Evaluasi Model

3.3.8.1 Koefisien Determinasi

Koefisien determinasi digunakan untuk mengevaluasi kemampuan model regresi logistik dalam menjelaskan variabilitas variabel respon.

# PSEUDO R2
pR2(model)
## fitting null model for pseudo-r2
##           llh       llhNull            G2      McFadden          r2ML 
## -4369.7082983 -5177.9773763  1616.5381559     0.1560975     0.2811054 
##          r2CU 
##     0.3196976
pR2_df <- data.frame(
  Indeks = c("Log-Likelihood (Full Model)", "Log-Likelihood (Null Model)",
             "Likelihood Ratio (G2)", "McFadden R²", "Cox & Snell R²", "Nagelkerke R²"),
  Nilai = c(-4369.7083, -5177.9774, 1616.5382, 0.1561, 0.2811, 0.3197)
)

kable(pR2_df,
      caption = "Tabel 3.30 Hasil Goodness of Fit Model Regresi Logistik Ordinal",
      align = "c") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.30 Hasil Goodness of Fit Model Regresi Logistik Ordinal
Indeks Nilai
Log-Likelihood (Full Model) -4369.7083
Log-Likelihood (Null Model) -5177.9774
Likelihood Ratio (G2) 1616.5382
McFadden R² 0.1561
Cox & Snell R² 0.2811
Nagelkerke R² 0.3197

Nilai pseudo R-square menunjukkan kemampuan model dalam menjelaskan variasi data. Nilai McFadden R² sebesar 0.1561 menunjukkan bahwa model memiliki kemampuan penjelasan yang cukup baik dalam konteks regresi logistik. Selain itu, nilai Nagelkerke R² sebesar 0.3197 menunjukkan bahwa sekitar 31.97% variasi pada variabel respon dapat dijelaskan oleh variabel prediktor dalam model.

3.3.8.2 Confusion Matrix dan Akurasi

Kinerja klasifikasi model regresi logistik dievaluasi menggunakan confusion matrix dan nilai akurasi. Confusion matrix digunakan untuk membandingkan hasil prediksi model dengan kondisi aktual sehingga dapat diketahui jumlah observasi yang diklasifikasikan dengan benar maupun salah.

# CONFUSION MATRIX
pred <- predict(model, type = "class")

cm <- table(Predicted = pred, Actual = wine$quality_ord)
cm
##          Actual
## Predicted  low medium high
##    low     936    466   58
##    medium  682   1456  683
##    high     22    276  319
cm_df <- data.frame(
  Predicted = c("Low", "Medium", "High"),
  Low = c(936, 682, 22),
  Medium = c(466, 1456, 276),
  High = c(58, 683, 319)
)

kable(cm_df,
      caption = "Tabel 3.31 Confusion Matrix Model Regresi Logistik Ordinal",
      align = "c") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.31 Confusion Matrix Model Regresi Logistik Ordinal
Predicted Low Medium High
Low 936 466 58
Medium 682 1456 683
High 22 276 319

Berdasarkan confusion matrix, model menunjukkan kemampuan prediksi yang cukup baik terutama pada kelas medium yang memiliki jumlah prediksi benar tertinggi (1456 observasi). Namun, masih terdapat kesalahan klasifikasi pada seluruh kelas, terutama pada kelas medium yang cukup sering diprediksi sebagai low maupun high. Hal ini menunjukkan bahwa meskipun model telah mampu menangkap pola umum data, masih terdapat tumpang tindih antar kategori.

confusionMatrix(cm)
## Confusion Matrix and Statistics
## 
##          Actual
## Predicted  low medium high
##    low     936    466   58
##    medium  682   1456  683
##    high     22    276  319
## 
## Overall Statistics
##                                           
##                Accuracy : 0.5535          
##                  95% CI : (0.5394, 0.5675)
##     No Information Rate : 0.4488          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.2733          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
## 
## Statistics by Class:
## 
##                      Class: low Class: medium Class: high
## Sensitivity              0.5707        0.6624     0.30094
## Specificity              0.8392        0.4944     0.92236
## Pos Pred Value           0.6411        0.5161     0.51702
## Neg Pred Value           0.7952        0.6428     0.82691
## Prevalence               0.3348        0.4488     0.21641
## Detection Rate           0.1911        0.2973     0.06513
## Detection Prevalence     0.2981        0.5759     0.12597
## Balanced Accuracy        0.7049        0.5784     0.61165
stats_df <- data.frame(
  Metric = c("Accuracy", "Kappa", "No Information Rate", "P-Value (Acc > NIR)",
             "Sensitivity (Low)", "Sensitivity (Medium)", "Sensitivity (High)",
             "Specificity (Low)", "Specificity (Medium)", "Specificity (High)",
             "Balanced Accuracy (Low)", "Balanced Accuracy (Medium)", "Balanced Accuracy (High)"),
  Value = c(0.5535, 0.2733, 0.4488, 2.2e-16,
            0.5707, 0.6624, 0.30094,
            0.8392, 0.4944, 0.92236,
            0.7049, 0.5784, 0.61165)
)

kable(stats_df,
      caption = "Tabel 3.32 Evaluasi Kinerja Model Regresi Logistik Ordinal",
      align = "c") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1a3a2a")
Tabel 3.32 Evaluasi Kinerja Model Regresi Logistik Ordinal
Metric Value
Accuracy 0.55350
Kappa 0.27330
No Information Rate 0.44880
P-Value (Acc > NIR) 0.00000
Sensitivity (Low) 0.57070
Sensitivity (Medium) 0.66240
Sensitivity (High) 0.30094
Specificity (Low) 0.83920
Specificity (Medium) 0.49440
Specificity (High) 0.92236
Balanced Accuracy (Low) 0.70490
Balanced Accuracy (Medium) 0.57840
Balanced Accuracy (High) 0.61165

Berdasarkan hasil evaluasi model, diperoleh nilai accuracy sebesar 0.5535, yang menunjukkan bahwa model mampu mengklasifikasikan sekitar 55.35% data dengan benar. Nilai ini lebih tinggi dibandingkan No Information Rate (0.4488), sehingga model memiliki kinerja yang lebih baik dibandingkan model tebakan acak.

Nilai Kappa sebesar 0.2733 menunjukkan tingkat kesepakatan yang tergolong fair antara hasil prediksi dan data aktual. Selain itu, performa model pada setiap kelas menunjukkan variasi, di mana kelas medium memiliki sensitivitas tertinggi (0.6624), sedangkan kelas high memiliki sensitivitas terendah (0.30094), yang menunjukkan bahwa model masih kurang optimal dalam mengenali kelas high.

Secara keseluruhan, model memiliki performa moderat dengan kemampuan prediksi yang lebih baik pada kelas low dan medium dibandingkan kelas high.

3.3.9 Kesimpulan

Berdasarkan hasil analisis regresi logistik ordinal, model dengan variabel prediktor terbukti signifikan dalam menjelaskan variabel quality, serta seluruh variabel independen juga berpengaruh secara parsial.

Berdasarkan odds ratio, alcohol dan sulphates meningkatkan peluang kualitas yang lebih tinggi, sedangkan density dan volatile acidity cenderung menurunkan peluang tersebut.

Hasil uji Brant menunjukkan bahwa asumsi proportional odds belum sepenuhnya terpenuhi, sehingga terdapat indikasi bahwa pengaruh salah satu variabel tidak konsisten pada seluruh kategori respon. Namun, pelanggaran ini bersifat parsial sehingga model masih dapat digunakan dengan interpretasi yang lebih hati-hati.

Secara keseluruhan, model memiliki performa moderat dengan akurasi 55.35% dan Kappa 0.2733.

3.4 Regresi Poisson

3.4.1 Deskripsi Kasus

Penelitian ini menggunakan data sewa sepeda di Seoul Bike Sharing System yang merekam jumlah penyewaan sepeda per jam selama periode tertentu. Variabel respons dalam analisis ini adalah jumlah sepeda yang disewa (Rented Bike Count) yang merupakan data berbentuk hitungan (count data) dengan nilai nonnegatif.

Tujuan analisis ini adalah untuk mengetahui faktor-faktor yang memengaruhi jumlah penyewaan sepeda berdasarkan kondisi lingkungan dan waktu. Variabel yang digunakan ditampilkan pada tabel berikut:

Tabel 3.33 Variabel Regresi Poisson

Variabel Tipe Keterangan
Rented Bike Count Numerik (Count) Jumlah sepeda yang disewa (data hitungan nonnegatif)
Hour Kategorik / Numerik Diskrit Waktu dalam satu hari (0–23 jam)
Temperature (°C) Numerik Suhu udara
Humidity (%) Numerik Kelembapan udara
Wind speed (m/s) Numerik Kecepatan angin
Solar Radiation (MJ/m²) Numerik Intensitas radiasi matahari
Rainfall (mm) Numerik Curah hujan
Snowfall (cm) Numerik Ketebalan salju
Seasons Kategorik Musim (Spring, Summer, Autumn, Winter)

Data ini bersifat time-based dan dipengaruhi oleh kondisi cuaca serta pola aktivitas masyarakat yang berubah sepanjang hari dan musim. Oleh karena itu, pendekatan regresi Poisson digunakan sebagai metode awal untuk memodelkan hubungan antara variabel respons dan variabel prediktor, karena sesuai untuk data berbentuk hitungan.

Analisis ini diharapkan dapat memberikan gambaran mengenai faktor dominan yang memengaruhi fluktuasi jumlah penyewaan sepeda serta membantu memahami pola penggunaan layanan transportasi berbasis sepeda.

3.4.2 Preparasi Data

Tahap awal dalam analisis ini adalah preparasi data yang bertujuan untuk memastikan dataset siap digunakan dalam pemodelan regresi Poisson.

# LOAD PACKAGE
library(dplyr)
library(ggplot2)
library(tibble)
library(pscl)

# PREPARASI DATA 
data <- read_excel(
  "C:/Users/Asus/Downloads/DATA ADK.xlsx",
  sheet = "Poisson"
  )
## New names:
## • `FALSE` -> `FALSE...4`
## • `FALSE` -> `FALSE...8`
bike <- data %>%
  dplyr::rename(
    temperature = `FALSE...4`,
    humidity = `Humidity(%)`,
    wind = `Wind speed (m/s)`,
    solar = `Solar Radiation (MJ/m2)`,
    rainfall = `Rainfall(mm)`,
    snowfall = `Snowfall (cm)`,
    rented_bike = `Rented Bike Count`
  ) %>%
  dplyr::select(
    rented_bike,
    Hour,
    temperature,
    humidity,
    wind,
    solar,
    rainfall,
    snowfall,
    Seasons
  ) %>%
  mutate(
    Seasons = factor(Seasons)
  )

3.4.3 Eksplorasi Data

Eksplorasi data dilakukan untuk memberikan gambaran awal mengenai karakteristik data serta hubungan antara variabel penjelas dengan variabel respon.

# EKSPLORASI DATA
ggplot(bike, aes(x = rented_bike)) +
  geom_histogram(bins = 30, fill = "#2f7f73") +
  labs(
    title = "Distribusi Jumlah Penyewaan Sepeda",
    x = "Rented Bike Count",
    y = "Frekuensi"
  )

Distribusi Rented Bike Count cenderung miring ke kanan (right-skewed). Sebagian besar pengamatan memiliki jumlah penyewaan yang rendah hingga sedang, sementara hanya sedikit pengamatan dengan jumlah penyewaan yang sangat tinggi.

3.4.4 Pemodelan Regresi Poisson

Regresi Poisson digunakan untuk memodelkan hubungan antara variabel respon yang berbentuk data cacahan (count data) dengan satu atau lebih variabel prediktor. Model ini mengasumsikan bahwa variabel respon mengikuti distribusi Poisson, di mana nilai yang diamati merupakan jumlah kejadian yang terjadi dalam suatu interval waktu, ruang, atau kondisi tertentu. Hubungan antara nilai harapan variabel respon dan variabel prediktor dinyatakan melalui fungsi log link.

Secara umum, model regresi Poisson dapat dituliskan sebagai berikut:

\[ \log(\mu_i)=\beta_0+\beta_1X_{1i}+\beta_2X_{2i}+\cdots+\beta_pX_{pi} \]

dengan:

  • \(\mu_i = E(Y_i)\) = nilai harapan dari variabel respon ke-\(i\)
  • \(Y_i\) = variabel respon berupa data cacahan
  • \(X_{ji}\) = variabel prediktor ke-\(j\) pada pengamatan ke-\(i\)
  • \(\beta_0\) = intersep
  • \(\beta_j\) = koefisien regresi untuk variabel prediktor ke-\(j\)
fit_pois <- glm(
  rented_bike ~ Hour + temperature + humidity +
    wind + solar + rainfall + snowfall + Seasons,
  data = bike,
  family = poisson(link = "log")
)

summary(fit_pois)
## 
## Call:
## glm(formula = rented_bike ~ Hour + temperature + humidity + wind + 
##     solar + rainfall + snowfall + Seasons, family = poisson(link = "log"), 
##     data = bike)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    6.264e+00  2.471e-03 2535.43   <2e-16 ***
## Hour           4.535e-02  6.905e-05  656.82   <2e-16 ***
## temperature    2.879e-02  7.629e-05  377.41   <2e-16 ***
## humidity      -9.396e-03  2.818e-05 -333.42   <2e-16 ***
## wind           1.240e-02  4.515e-04   27.47   <2e-16 ***
## solar         -6.517e-02  5.661e-04 -115.12   <2e-16 ***
## rainfall      -5.344e-01  2.179e-03 -245.29   <2e-16 ***
## snowfall      -9.466e-02  1.971e-03  -48.02   <2e-16 ***
## SeasonsSpring -7.839e-02  1.110e-03  -70.64   <2e-16 ***
## SeasonsSummer -3.651e-02  1.341e-03  -27.24   <2e-16 ***
## SeasonsWinter -8.960e-01  2.096e-03 -427.50   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4979261  on 8759  degrees of freedom
## Residual deviance: 2162206  on 8749  degrees of freedom
## AIC: 2229321
## 
## Number of Fisher Scoring iterations: 7
coef_df <- data.frame(
  Variabel = c(
    "Intercept",
    "Hour",
    "Temperature",
    "Humidity",
    "Wind",
    "Solar",
    "Rainfall",
    "Snowfall",
    "Spring",
    "Summer",
    "Winter"
  ),
  Estimate = c(
    6.264,
    0.04535,
    0.02879,
    -0.009396,
    0.01240,
    -0.06517,
    -0.5344,
    -0.09466,
    -0.07839,
    -0.03651,
    -0.8960
  )
)

kable(
  coef_df,
  caption = "Tabel 3.34 Hasil Estimasi Parameter Regresi Poisson",
  align = "c"
) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1a3a2a"
  )
Tabel 3.34 Hasil Estimasi Parameter Regresi Poisson
Variabel Estimate
Intercept 6.264000
Hour 0.045350
Temperature 0.028790
Humidity -0.009396
Wind 0.012400
Solar -0.065170
Rainfall -0.534400
Snowfall -0.094660
Spring -0.078390
Summer -0.036510
Winter -0.896000

3.4.5 Pengujian Hipotesis

3.4.5.1 Uji Simultan (Likelihood Ratio Test)

Uji signifikansi simultan dilakukan menggunakan Likelihood Ratio Test (LRT) untuk mengetahui apakah variabel prediktor secara bersama-sama berpengaruh terhadap variabel respon.

Hipotesis yang digunakan adalah:

  • \(H_0\) : \(\beta_1 = \beta_2 = \beta_3 = \cdots = \beta_k = 0\) (tidak terdapat pengaruh signifikan secara simultan antara variabel prediktor terhadap variabel respon)

  • \(H_1\) : paling sedikit terdapat satu \(\beta_j \neq 0\) (terdapat minimal satu variabel prediktor yang berpengaruh signifikan terhadap variabel respon)

# UJI SIMULTAN 
model_null <- glm(
  rented_bike ~ 1,
  data = bike,
  family = poisson
)

anova(model_null, fit_pois, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: rented_bike ~ 1
## Model 2: rented_bike ~ Hour + temperature + humidity + wind + solar + 
##     rainfall + snowfall + Seasons
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      8759    4979261                          
## 2      8749    2162206 10  2817056 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrt_df <- data.frame(
  Statistik_G = 2817056,
  Df = 10,
  P_value = "< 0.001",
  Keputusan = "Tolak H0"
)

kable(
  lrt_df,
  caption = "Tabel 3.35 Hasil Uji Simultan Regresi Poisson",
  align = "c"
) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1a3a2a"
  )
Tabel 3.35 Hasil Uji Simultan Regresi Poisson
Statistik_G Df P_value Keputusan
2817056 10 < 0.001 Tolak H0

Berdasarkan hasil uji simultan diperoleh nilai statistik G sebesar 2.817.056 dengan derajat bebas 10 dan p-value < 0.001. Karena p-value < 0.05, maka \(H_0\) ditolak. Dengan demikian, variabel Hour, temperature, humidity, wind, solar, rainfall,snowfall, dan Seasons secara simultan berpengaruh signifikan terhadap jumlah rented bike.

3.4.5.2 Uji Parsial (Wald Test)

Uji signifikansi parsial dilakukan menggunakan uji Wald untuk mengetahui pengaruh masing-masing variabel prediktor terhadap variabel respon.

Hipotesis (untuk setiap variabel prediktor):

  • \(H_0 : \beta_i = 0\) (variabel ke-\(i\) tidak berpengaruh terhadap variabel respon)

  • \(H_1 : \beta_i \neq 0\) (variabel ke-\(i\) berpengaruh terhadap variabel respon)

# UJI PARSIAL
coef(summary(fit_pois))
##                   Estimate   Std. Error    z value      Pr(>|z|)
## (Intercept)    6.263887336 2.470546e-03 2535.42588  0.000000e+00
## Hour           0.045350641 6.904589e-05  656.81879  0.000000e+00
## temperature    0.028791791 7.628797e-05  377.40930  0.000000e+00
## humidity      -0.009396041 2.818126e-05 -333.41451  0.000000e+00
## wind           0.012402860 4.515529e-04   27.46712 4.338796e-166
## solar         -0.065169559 5.660893e-04 -115.12240  0.000000e+00
## rainfall      -0.534398096 2.178664e-03 -245.28707  0.000000e+00
## snowfall      -0.094664774 1.971300e-03  -48.02151  0.000000e+00
## SeasonsSpring -0.078393324 1.109687e-03  -70.64451  0.000000e+00
## SeasonsSummer -0.036513653 1.340598e-03  -27.23684 2.379751e-163
## SeasonsWinter -0.895960774 2.095794e-03 -427.50414  0.000000e+00
parsial_df <- data.frame(
  Variabel = c(
    "Intercept",
    "Hour",
    "Temperature",
    "Humidity",
    "Wind",
    "Solar",
    "Rainfall",
    "Snowfall",
    "Spring",
    "Summer",
    "Winter"
  ),
  Estimate = c(
    6.263887336,
    0.045350641,
    0.028791791,
    -0.009396041,
    0.012402860,
    -0.065169559,
    -0.534398096,
    -0.094664774,
    -0.078393324,
    -0.036513653,
    -0.895960774
  ),
  `p-value` = c(
    0,
    0,
    0,
    0,
    4.338796e-166,
    0,
    0,
    0,
    0,
    2.379751e-163,
    0
  )
)

kable(
  parsial_df,
  caption = "Tabel 3.36 Hasil Uji Parsial Regresi Poisson",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1a3a2a"
  ) %>%
  row_spec(
    which(parsial_df$`p-value` >= 0.05),
    background = "#ffcccc"
  )
Tabel 3.36 Hasil Uji Parsial Regresi Poisson
Variabel Estimate p.value
Intercept 6.2638873 0
Hour 0.0453506 0
Temperature 0.0287918 0
Humidity -0.0093960 0
Wind 0.0124029 0
Solar -0.0651696 0
Rainfall -0.5343981 0
Snowfall -0.0946648 0
Spring -0.0783933 0
Summer -0.0365137 0
Winter -0.8959608 0

Berdasarkan hasil uji parsial, seluruh variabel memiliki nilai p-value kurang dari 0.05 sehingga \(H_0\) ditolak. Hal ini menunjukkan bahwa variabel Hour, temperature, humidity, wind, `solar, rainfall, snowfall, dan Seasons berpengaruh signifikan terhadap jumlah rented bike. Dengan demikian, seluruh variabel prediktor yang digunakan dalam model memberikan kontribusi yang signifikan dalam menjelaskan variasi jumlah sepeda yang disewa.

3.4.6 Pengujian Asumsi

3.4.6.1 Respons berupa nilai non-negatif

Analisis regresi Poisson digunakan untuk memodelkan data dengan variabel respon berupa data hitungan yang bernilai non-negatif. Oleh karena itu, sebelum dilakukan pemodelan, perlu dipastikan bahwa karakteristik data sesuai dengan asumsi yang mendasari regresi Poisson.

# PENGUJIAN ASUMSI

## Respons berupa hitungan nonnegatif
any(bike$rented_bike < 0)
## [1] FALSE

Kondisi awal data menunjukkan bahwa tidak terdapat nilai negatif pada variabel respon rented_bike (any = FALSE), sehingga asumsi dasar regresi Poisson terkait data hitungan bernilai non-negatif telah terpenuhi.

3.4.6.2 Observasi Independen

Asumsi independensi pada regresi Poisson mengacu pada desain data yang digunakan, yaitu setiap observasi dianggap saling bebas sehingga tidak terdapat ketergantungan antar pengamatan. Pada data ini, jumlah penyewaan sepeda pada satu periode waktu tidak boleh dipengaruhi oleh periode waktu lainnya, sesuai dengan rancangan pengumpulan data yang menganggap setiap pengamatan sebagai unit yang independen.

3.4.6.4 Equidispersion

Asumsi utama pada regresi Poisson adalah equidispersion, yaitu kondisi di mana nilai rata-rata (mean) dari variabel respon sama dengan variansinya.

## Equidispersion 
mean(bike$rented_bike)
## [1] 704.6021
var(bike$rented_bike)
## [1] 416021.7
pearson_chi2 <- sum(residuals(fit_pois, type = "pearson")^2)
df_resid <- fit_pois$df.residual

dispersion <- pearson_chi2 / df_resid
dispersion
## [1] 528009.8
if (dispersion < 1.5) {
  cat("Tidak ada indikasi overdispersion berat")
} else if (dispersion < 2.5) {
  cat("Ada indikasi overdispersion sedang")
} else {
  cat("Ada indikasi overdispersion kuat")
}
## Ada indikasi overdispersion kuat
library(kableExtra)

disp_df <- data.frame(
  `Chi-Square` = c(528009.8),
  df = c(8759),
  Kesimpulan = ("Ada indikasi overdispersion kuat"),
  Nilai = c(528009.8, 8759, 60.28)
)

kable(
  disp_df,
  caption = "Tabel 3.37 Hasil Uji Overdispersion pada Model Regresi Poisson",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1a3a2a"
  )
Tabel 3.37 Hasil Uji Overdispersion pada Model Regresi Poisson
Chi.Square df Kesimpulan Nilai
528009.8 8759 Ada indikasi overdispersion kuat 528009.80
528009.8 8759 Ada indikasi overdispersion kuat 8759.00
528009.8 8759 Ada indikasi overdispersion kuat 60.28

Berdasarkan hasil uji overdispersion diperoleh nilai Pearson Chi-Square sebesar 528009.8 dengan derajat bebas 8759, sehingga menghasilkan rasio dispersi sebesar 60.28. Karena nilai tersebut jauh lebih besar dari 1, maka dapat disimpulkan bahwa terdapat overdispersion yang kuat pada model regresi Poisson, yang menunjukkan bahwa asumsi kesetaraan antara mean dan varians tidak terpenuhi.

Meskipun demikian, model regresi Poisson tetap digunakan sebagai model dasar (baseline model) untuk menggambarkan hubungan awal antara variabel prediktor dan jumlah rented bike, serta memberikan gambaran arah pengaruh variabel melalui koefisiennya. Selain itu, hasil model ini juga dapat digunakan sebagai pembanding untuk pengembangan model yang lebih sesuai, seperti regresi Negative Binomial, terutama dalam kondisi terdapat overdispersion.

3.4.7 IRR

Dalam regresi Poisson, interpretasi koefisien regresi lebih mudah dilakukan menggunakan Incidence Rate Ratio (IRR), yaitu hasil transformasi eksponensial dari koefisien model. IRR menunjukkan perubahan relatif pada rata-rata kejadian akibat perubahan satu satuan variabel prediktor.

# IRR
pois_coef <- as.data.frame(coef(summary(fit_pois))) %>%
  tibble::rownames_to_column("parameter") %>%
  mutate(
    IRR = exp(Estimate),
    CI_low = exp(Estimate - 1.96 * `Std. Error`),
    CI_high = exp(Estimate + 1.96 * `Std. Error`),
    persen_perubahan = 100 * (IRR - 1)
  )

pois_coef
##        parameter     Estimate   Std. Error    z value      Pr(>|z|)         IRR
## 1    (Intercept)  6.263887336 2.470546e-03 2535.42588  0.000000e+00 525.2568264
## 2           Hour  0.045350641 6.904589e-05  656.81879  0.000000e+00   1.0463947
## 3    temperature  0.028791791 7.628797e-05  377.40930  0.000000e+00   1.0292103
## 4       humidity -0.009396041 2.818126e-05 -333.41451  0.000000e+00   0.9906480
## 5           wind  0.012402860 4.515529e-04   27.46712 4.338796e-166   1.0124801
## 6          solar -0.065169559 5.660893e-04 -115.12240  0.000000e+00   0.9369086
## 7       rainfall -0.534398096 2.178664e-03 -245.28707  0.000000e+00   0.5860219
## 8       snowfall -0.094664774 1.971300e-03  -48.02151  0.000000e+00   0.9096778
## 9  SeasonsSpring -0.078393324 1.109687e-03  -70.64451  0.000000e+00   0.9246007
## 10 SeasonsSummer -0.036513653 1.340598e-03  -27.23684 2.379751e-163   0.9641449
## 11 SeasonsWinter -0.895960774 2.095794e-03 -427.50414  0.000000e+00   0.4082152
##         CI_low     CI_high persen_perubahan
## 1  522.7195386 527.8064302    52425.6826402
## 2    1.0462531   1.0465363        4.6394705
## 3    1.0290564   1.0293642        2.9210281
## 4    0.9905932   0.9907027       -0.9352036
## 5    1.0115844   1.0133766        1.2480095
## 6    0.9358696   0.9379487       -6.3091412
## 7    0.5835248   0.5885297      -41.3978087
## 8    0.9061698   0.9131994       -9.0322169
## 9    0.9225919   0.9266139       -7.5399313
## 10   0.9616149   0.9666816       -3.5855069
## 11   0.4065418   0.4098955      -59.1784792
irr_df <- data.frame(
  Variabel = c(
    "(Intercept)", "Hour", "Temperature", "Humidity", "Wind",
    "Solar", "Rainfall", "Snowfall",
    "Spring", "Summer", "Winter"
  ),
  IRR = c(
    525.2568, 1.0463947, 1.0292103, 0.9906480, 1.0124801,
    0.9369086, 0.5860219, 0.9096778,
    0.9246007, 0.9641449, 0.4082152
  ),
  CI_Low = c(
    522.7195, 1.0462531, 1.0290564, 0.9905932, 1.0115844,
    0.9358696, 0.5835248, 0.9061698,
    0.9225919, 0.9616149, 0.4065418
  ),
  CI_High = c(
    527.8064, 1.0465363, 1.0293642, 0.9907027, 1.0133766,
    0.9379487, 0.5885297, 0.9131994,
    0.9266139, 0.9666816, 0.4098955
  ),
  Persen_Perubahan = c(
    52425.6826, 4.6394705, 2.9210281, -0.9352036, 1.2480095,
    -6.3091412, -41.3978087, -9.0322169,
    -7.5399313, -3.5855069, -59.17848
  )
)

kable(
  irr_df,
  caption = "Tabel 3.38 Nilai Incidence Rate Ratio (IRR) Model Regresi Poisson",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1a3a2a"
  )
Tabel 3.38 Nilai Incidence Rate Ratio (IRR) Model Regresi Poisson
Variabel IRR CI_Low CI_High Persen_Perubahan
(Intercept) 525.2568000 522.7195000 527.8064000 52425.6826000
Hour 1.0463947 1.0462531 1.0465363 4.6394705
Temperature 1.0292103 1.0290564 1.0293642 2.9210281
Humidity 0.9906480 0.9905932 0.9907027 -0.9352036
Wind 1.0124801 1.0115844 1.0133766 1.2480095
Solar 0.9369086 0.9358696 0.9379487 -6.3091412
Rainfall 0.5860219 0.5835248 0.5885297 -41.3978087
Snowfall 0.9096778 0.9061698 0.9131994 -9.0322169
Spring 0.9246007 0.9225919 0.9266139 -7.5399313
Summer 0.9641449 0.9616149 0.9666816 -3.5855069
Winter 0.4082152 0.4065418 0.4098955 -59.1784800

Berdasarkan nilai IRR, variabel Hour, temperature, dan wind memiliki nilai IRR > 1 yang menunjukkan bahwa peningkatan variabel tersebut akan meningkatkan laju jumlah rented bike. Sebaliknya, variabel humidity, solar, rainfall, snowfall, serta kategori musim Spring, Summer, dan Winter memiliki nilai IRR < 1 yang menunjukkan penurunan laju jumlah penyewaan sepeda. Variabel rainfall dan Winter memiliki pengaruh penurunan paling besar terhadap jumlah rented bike.

3.4.8 Evaluasi Model

# Deviance
null_deviance <- fit_pois$null.deviance
residual_deviance <- fit_pois$deviance
df_resid <- fit_pois$df.residual

deviance_ratio <- residual_deviance / df_resid

null_deviance
## [1] 4979261
residual_deviance
## [1] 2162206
deviance_ratio
## [1] 247.1375
# AIC
AIC(fit_pois)
## [1] 2229321
# Pseudo R-square
pseudo_r2 <- pR2(fit_pois)
## fitting null model for pseudo-r2
pseudo_r2
##           llh       llhNull            G2      McFadden          r2ML 
## -1.114650e+06 -2.523178e+06  2.817056e+06  5.582357e-01  1.000000e+00 
##          r2CU 
##  1.000000e+00
gof_df <- data.frame(
  Indikator = c(
    "Null Deviance",
    "Residual Deviance",
    "Deviance Ratio",
    "AIC",
    "McFadden R²",
    "Cragg & Uhler R² (Nagelkerke)"
  ),
  Nilai = c(
    4979261,
    2162206,
    247.1375,
    2229321,
    0.5582357,
    1
  )
)

kable(
  gof_df,
  caption = "Tabel 3.39 Hasil Goodness of Fit Model Regresi Poisson",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = TRUE,
    position = "center"
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1a3a2a"
  )
Tabel 3.39 Hasil Goodness of Fit Model Regresi Poisson
Indikator Nilai
Null Deviance 4.979261e+06
Residual Deviance 2.162206e+06
Deviance Ratio 2.471375e+02
AIC 2.229321e+06
McFadden R² 5.582357e-01
Cragg & Uhler R² (Nagelkerke) 1.000000e+00

Nilai McFadden R² sebesar 0.5582 menunjukkan bahwa model memiliki kemampuan penjelasan yang cukup baik terhadap variasi data. Sementara itu, nilai Cragg & Uhler R² (Nagelkerke) mendekati 1 menunjukkan kecocokan model yang sangat tinggi, meskipun nilai ini perlu diinterpretasikan dengan hati-hati karena adanya indikasi overdispersion pada model.

3.4.9 Kesimpulan

Berdasarkan hasil regresi Poisson, variabel Hour, temperature, humidity, wind, solar radiation, rainfall, snowfall, dan Seasons berpengaruh signifikan terhadap jumlah penyewaan sepeda (p-value < 0.05). Secara simultan, seluruh variabel prediktor juga berpengaruh signifikan terhadap rented bike count dengan nilai likelihood ratio sebesar 2817056 dan p-value < 0.001.

Namun, hasil pengujian asumsi menunjukkan adanya overdispersion yang sangat kuat, sehingga asumsi equidispersion pada regresi Poisson tidak terpenuhi. Oleh karena itu, hasil inferensi dari model Poisson perlu diinterpretasikan dengan hati-hati dan model alternatif seperti Negative Binomial Regression lebih direkomendasikan untuk memperoleh estimasi yang lebih reliabel.

4 Referensi

  1. Agresti, A. (2013). Categorical Data Analysis. Wiley.

  2. Agresti, A., & Finlay, B. (2009). Statistical Methods for the Social Sciences. Pearson.

  3. Chicco, D., & Jurman, G. (2020). Heart Failure Clinical Records [Data set]. UCI Machine Learning Repository. https://doi.org/10.24432/C5Z89R

  4. Cho, S. (2020). Seoul Bike Sharing Demand [Data set]. UCI Machine Learning Repository. https://doi.org/10.24432/C5F62R

  5. Cortez, P., Cerdeira, A., Almeida, F., Matos, T., & Reis, J. (2009). Wine Quality [Data set]. UCI Machine Learning Repository. https://doi.org/10.24432/C56S3T

  6. McCullagh, P., & Nelder, J. A. (1989). Generalized linear models (2nd ed.). Chapman & Hall.

  7. Realinho, V., Vieira Martins, M., Machado, J., & Baptista, L. (2021). Predict Students’ Dropout and Academic Success [Data set]. UCI Machine Learning Repository. https://doi.org/10.24432/C5MC89


:::