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.
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.
Analisis data kategori banyak digunakan dalam berbagai bidang penelitian. Beberapa contoh penerapannya adalah sebagai berikut:
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.
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}\)
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.
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}. \]
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}} \]
Ukuran asosiasi digunakan untuk mengukur kekuatan hubungan antara dua variabel dalam tabel kontingensi.
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}. \]
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:
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:
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.}} \]
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:
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:
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:
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:
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.
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:
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)} \]
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))} \]
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))} \]
Uji independensi digunakan untuk mengetahui apakah terdapat hubungan signifikan antara dua variabel kategorik pada tabel kontingensi. Hipotesis yang diuji adalah:
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:
\[ E_{ij} = \frac{n_{i.} n_{.j}}{n} \]
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:
\[\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}\]
\[R_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}}\]
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)\)
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:
p-value uji Fisher didapatkan dari total probabilitas semua tabel yang sama atau lebih ekstrem daripada tabel yang diamati
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
4. Deteksi Outlier
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
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 |
\[ 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 \]
\[ 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 \]
\[ \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 \]
\[ \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 \]
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 \]
# ================================
# 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
## [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
## 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"
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.
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.
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:
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.
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:
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:
Hipotesis:
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
##
## 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.
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:
Perhitungan Manual:
Langkah 1: Lakukan Chi-Square Secara Keseluruhan
\[ \chi^2 = 12.57 \]
\[ 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:
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.
Misalkan diberikan contoh sebagai berikut:
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
## [1] 0.01380413
## [1] 7.26533e-05
## [1] 0.001598373
## [1] 0.01380413
## [1] 0.06211857
## [1] 0.162464
## [1] 0.2599423
## [1] 0.2599423
## [1] 0.162464
## [1] 0.06211857
## [1] 0.01380413
## [1] 0.001598373
## [1] 7.26533e-05
## [1] 0.01380413
## [1] 0.03094
##
## 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.
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:
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.
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 |
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
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
## 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:
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 |
# 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:
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:
##
## Attaching package: 'epitools'
## The following object is masked from 'package:survival':
##
## ratetable
## 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
## Warning: package 'ggplot2' was built under R version 4.4.3
## Warning: package 'kableExtra' was built under R version 4.4.3
## Warning: package 'MASS' was built under R version 4.4.3
## 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
## 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
## Warning: package 'corrplot' was built under R version 4.4.3
## corrplot 0.95 loaded
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:
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} \]
\[\hat{p}_1 = \frac{a}{n_1} = \frac{688}{1338}\]
\[\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%)
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).
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
[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 | p̂ | 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).
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) ===
## RD = 0.5142 - 0.2625 = 0.2517
## SE(RD) = 0.0511
## 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).
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) ===
## RR = 0.5142 / 0.2625 = 1.9589
## SE(ln RR) = 0.1893
## 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.
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) ===
## RR = 0.5142 / 0.2625 = 2.9738
## SE(ln OR) = 0.2599
## 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) |
Hipotesis:
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) ===
## Proporsi gabungan p̄ = (688 + 21) / 1418 = 0.5000
## SE gabungan = 0.0575
## Statistik Z = 4.3737
## p-value (2-sisi) = 0.000012
## 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() ---
##
## 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.
Hipotesis:
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}\]
## === Uji Chi-Square Independensi ===
##
## Frekuensi Observasi:
## [,1] [,2]
## [1,] 688 650
## [2,] 21 59
##
## Frekuensi Harapan (Expected):
## [,1] [,2]
## [1,] 669 669
## [2,] 40 40
##
## Statistik χ² = 19.1292
## df = 1
## 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.
Uji Likelihood Ratio merupakan alternatif dari uji Chi-Square untuk menguji asosiasi antar variabel kategorik.
Hipotesis:
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²) ===
##
## 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
##
## G² = 2 × Σ(O × ln(O/E)) = 2 × 9.9390 = 19.8780
## df = 1
## 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.
Fisher exact test digunakan ketika asumsi chi-square tidak terpenuhi (frekuensi harapan < 5), dengan menghitung probabilitas eksak berdasarkan distribusi hipergeometrik.
Hipotesis:
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 Exact Test ===
## Odds Ratio (MLE) = 2.9716
## 95% CI OR = (1.7556, 5.2107)
## 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.
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:
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")Berdasarkan seluruh analisis pada Kasus 1, diperoleh kesimpulan sebagai berikut:
Estimasi proporsi: Proporsi kanker paru pada kelompok Smoker (\(\hat{p}_1 = 0.5142\)) jauh lebih tinggi dibandingkan Non-Smoker (\(\hat{p}_2 = 0.2625\)).
Ukuran asosiasi:
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).
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.
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 |
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.
Hipotesis:
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.
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.
Partisi chi-square membagi statistik uji total menjadi komponen ortogonal yang dapat diinterpretasikan secara terpisah.
# 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 ===
## Democrat Republican
## Female 495 272
## Male 330 265
##
## Frekuensi Harapan:
## 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.
# 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 ===
## Dem+Rep Independent
## Female 767 590
## Male 595 498
##
## Frekuensi Harapan:
## 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.
# 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.
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()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:
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.
Residual analisis: Sel-sel yang paling berkontribusi terhadap asosiasi adalah yang terkait dengan afiliasi Republican dan Democrat, bukan Independent.
Partisi chi-square:
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.
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.
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) |
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
## 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
## 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.
## Warning: package 'ResourceSelection' was built under R version 4.4.3
## ResourceSelection 0.3-6 2023-06-27
## 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
## Warning: package 'readxl' was built under R version 4.4.3
## Warning: package 'rcompanion' was built under R version 4.4.3
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.
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.
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")| 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.
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:
##
## 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")| Variabel | Estimate |
|---|---|
| (Intercept) | 9.1555 |
| age | 0.0486 |
| ejection_fraction | -0.0933 |
| serum_creatinine | 0.7282 |
| serum_sodium | -0.0626 |
| time | -0.0178 |
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")| 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.
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)
##
## 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")| 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.
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:
## (Intercept) age ejection_fraction serum_creatinine
## 9466.1478575 1.0498206 0.9108778 2.0712604
## serum_sodium time
## 0.9392723 0.9823541
## 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")| 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.
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")| 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.
Koefisien determinasi digunakan untuk mengevaluasi kemampuan model regresi logistik dalam menjelaskan variabilitas variabel respon.
## $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")| 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.
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")| 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.
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")| 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).
## 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")| 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.
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.
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 |
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
## Registered S3 method overwritten by 'micsr':
## method from
## summary.effects lava
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.
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
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:
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
## 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")| 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 |
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)
## # weights: 6 (2 variable)
## initial value 3890.186114
## final value 3611.622899
## converged
## 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"
)| 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.
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")| 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.
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:
## (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")| 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).
Variabel respon berupa status mahasiswa (Dropout, Enrolled, Graduate) yang merupakan kategori nominal tanpa urutan alami sehingga memenuhi asumsi regresi logistik multinomial.
Observasi diasumsikan saling independen karena setiap mahasiswa hanya tercatat satu kali dalam dataset.
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")| 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.
Tidak ditemukan indikasi perfect separation karena model berhasil diestimasi tanpa peringatan dan tidak terdapat koefisien yang bernilai sangat besar.
# 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")| 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")| 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.
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:
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")| 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.
Koefisien determinasi digunakan untuk mengevaluasi kemampuan model regresi logistik dalam menjelaskan variabilitas variabel respon.
## 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")| 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.
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
## 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")| 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.
## 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")| 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.
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.
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:
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 |
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
)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.
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:
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")| Variabel | Estimate |
|---|---|
| Alcohol | 0.7479 |
| Density | -211.5570 |
| Volatile acidity | -5.4606 |
| Residual sugar | 0.1405 |
| Sulphates | 1.8171 |
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"
)| 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.
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")| 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 | medium | -203.0532477 | 0.1772673 | -1145.46395 | 0 | Signifikan |
| medium | 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.
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:
## 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")| 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.Uji Brant digunakan untuk menguji asumsi proportional odds pada model regresi logistik ordinal.
Hipotesis:
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")| 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.
Koefisien determinasi digunakan untuk mengevaluasi kemampuan model regresi logistik dalam menjelaskan variabilitas variabel respon.
## 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")| 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.
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")| 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.
## 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")| 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.
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.
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.
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)
)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.
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:
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"
)| 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 |
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"
)| 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.
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)
## 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"
)| 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.
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.
## [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.
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.
Asumsi bahwa mean variabel respon dimodelkan dengan fungsi link log telah terpenuhi melalui spesifikasi model regresi Poisson, yaitu dengan menghubungkan kovariat secara linear pada skala logaritmik.
Asumsi utama pada regresi Poisson adalah equidispersion, yaitu kondisi di mana nilai rata-rata (mean) dari variabel respon sama dengan variansinya.
## [1] 704.6021
## [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"
)| 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.
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"
)| 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.
# 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
## [1] 2162206
## [1] 247.1375
## [1] 2229321
## fitting null model for 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"
)| 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.
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.
Agresti, A. (2013). Categorical Data Analysis. Wiley.
Agresti, A., & Finlay, B. (2009). Statistical Methods for the Social Sciences. Pearson.
Chicco, D., & Jurman, G. (2020). Heart Failure Clinical Records [Data set]. UCI Machine Learning Repository. https://doi.org/10.24432/C5Z89R
Cho, S. (2020). Seoul Bike Sharing Demand [Data set]. UCI Machine Learning Repository. https://doi.org/10.24432/C5F62R
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
McCullagh, P., & Nelder, J. A. (1989). Generalized linear models (2nd ed.). Chapman & Hall.
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
:::