1 Pendahuluan

Dalam era digital dan perkembangan teknologi informasi yang pesat, data menjadi salah satu sumber daya paling berharga dalam pengambilan keputusan di berbagai sektor. Data tidak hanya hadir dalam bentuk angka atau nilai kontinu, tetapi juga dalam bentuk kategori-kategori tertentu yang merepresentasikan atribut, karakteristik, atau kelas dari suatu fenomena.

Data yang berbentuk kategori inilah yang dikenal dengan istilah data kategorik. Contoh sederhana yang kita temui sehari-hari antara lain: jenis kelamin (laki-laki atau perempuan), status perkawinan (menikah atau belum menikah), tingkat kepuasan pelanggan (sangat puas, puas, tidak puas), hingga jenis produk yang dipilih konsumen.

Namun, data kategorik memiliki karakteristik khusus sehingga memerlukan pendekatan analisis yang berbeda dengan data numerik. Tidak bisa dihitung rata-rata atau standar deviasinya seperti data kuantitatif, data kategorik dianalisis menggunakan teknik yang fokus pada hubungan antar kategori, distribusi frekuensi, dan prediksi kejadian berdasarkan kelas data tertentu.

Oleh karena itu, Analisis Data Kategorik (ADK) hadir sebagai solusi utama untuk memahami data yang bersifat non-numerik ini. Dengan menggunakan ADK, peneliti dan praktisi dapat menggali lebih dalam pola, tren, serta keterkaitan antar kategori dalam suatu dataset.

1.1 Tujuan Analisis Data Kategori

Berikut adalah tujuan utama ADK yang telah disusun secara sistematis:

No. Tujuan Penjelasan
1 Mengidentifikasi Pola dan Tren Membantu menemukan pola tersembunyi, misalnya keterkaitan antara preferensi pelanggan dan usia.
2 Menganalisis Hubungan Antarvariabel Menilai apakah dua variabel kategorik saling berhubungan atau bersifat independen. Contoh: Apakah tingkat pendidikan memengaruhi preferensi politik?
3 Membantu Pengambilan Keputusan Data kategorik dapat membantu menyusun strategi kebijakan berbasis data demografis atau segmentasi pasar.
4 Mengembangkan Model Prediktif Digunakan untuk membangun model prediksi seperti regresi logistik untuk memprediksi kemungkinan hasil biner.
5 Membantu Validasi Hipotesis Menguji apakah ada perbedaan signifikan antar kategori dalam penelitian eksperimen atau observasi.

1.2 Definisi dan Ruang Lingkup Analisis Data Kategori

1.2.1 Definisi

Analisis Data Kategorik (ADK) adalah cabang statistik yang mempelajari data yang berbentuk kategori atau klasifikasi, untuk memahami hubungan, pola, dan prediksi berdasarkan kategori tersebut. Data ini tidak memiliki nilai numerik murni, melainkan menunjukkan kelompok tertentu (misalnya: Ya/Tidak, A/B/C). Agar dapat dianalisis, kategori biasanya dikodekan ke angka, namun maknanya tetap bersifat simbolik.

1.2.2 Jenis-jenis Data Kategorik

1. Nominal

  • Kategori tanpa urutan (setara).
  • Contoh: Golongan darah (A, B, AB, O), warna favorit, jenis kendaraan.

2. Ordinal

  • Kategori memiliki urutan, tapi jaraknya tidak terukur pasti.
  • Contoh: Tingkat pendidikan (SD < SMP < SMA), skala kepuasan.

3. Biner (Dichotomous)

  • Hanya dua kategori.
  • Contoh: Ya/Tidak, Lulus/Gagal, Positif/Negatif.

4. Multikategori

  • Lebih dari dua kategori, bisa nominal atau ordinal.
  • Contoh: Jenis pekerjaan, tingkat kepedulian lingkungan, platform media sosial.

Data kategorik banyak digunakan dalam survei, riset sosial, medis, dan pemasaran. ADK menyediakan metode khusus untuk menganalisis data jenis ini, berbeda dengan data numerik biasa.

1.3 Perbedaan dengan Data Kuantitatif

Aspek Data Kategorik Data Kuantitatif
Jenis Data Nominal, ordinal, biner, multikategori Diskrit, kontinu
Contoh Status vaksinasi, jenis pekerjaan Tinggi badan, berat badan
Pengukuran Tidak ada nilai numerik intrinsik Nilai numerik dengan satuan
Metode Analisis Uji Chi-Square, Regresi Logistik, Analisis Correspondence Regresi Linear, ANOVA, uji t
Visualisasi Diagram batang, pie chart Histogram, scatter plot
Pertanyaan yang Dijawab Apakah terdapat hubungan antara kategori? Berapa rata-rata nilai populasi?

2 Metode dalam Analisis Data Kategori

2.1 Tabel Kontingensi dan Uji Chi-Square

Penjelasan:

Tabel kontingensi adalah alat statistik yang digunakan untuk menampilkan hubungan antara dua atau lebih variabel kategorik dalam bentuk tabel distribusi frekuensi. Tabel ini membantu kita memahami pola, kecenderungan, atau kemungkinan adanya hubungan antara kategori yang berbeda. Dalam analisis data kategorik, tabel kontingensi sering menjadi langkah awal sebelum dilakukan uji lanjutan seperti uji Chi-Square, untuk menguji apakah hubungan antara variabel-variabel tersebut bersifat signifikan secara statistik.

Fungsi Utama Tabel Kontingensi:

  • Menyusun data kategorik agar lebih mudah dianalisis.
  • Mengidentifikasi pola atau kecenderungan hubungan antar kategori.
  • Dasar untuk pengujian hubungan menggunakan uji statistik.

Contoh Kasus:

Survei mengenai hubungan antara status pekerjaan dan pilihan transportasi utama.

# Pastikan kableExtra sudah tersedia
library(knitr)
## Warning: package 'knitr' was built under R version 4.3.3
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.3.3
# Buat data
data_transport <- matrix(c(
  25, 40, 15, 80,
   5, 30, 25, 60,
  20, 15,  5, 40,
  50, 85, 45, 180
), ncol = 4, byrow = TRUE)

colnames(data_transport) <- c("Mobil", "Sepeda Motor", "Transportasi Umum", "Total")
rownames(data_transport) <- c("Karyawan", "Mahasiswa", "Wirausaha", "Total")

# Tampilkan tabel dengan kable
kable(data_transport, format = "latex", booktabs = TRUE, align = "c",
      caption = "Distribusi Status Pekerjaan dan Pilihan Transportasi") %>%
  kable_styling(latex_options = c("hold_position", "striped"), full_width = FALSE)

Hipotesis:

  • H0 : Tidak ada hubungan

  • H1 : Ada hubungan

Rumus: \[ \chi^2 = \sum \frac{(O - E)^2}{E} \]

Kriteria Uji:

Jika p-value < 0.05 dan \(\chi^2 \text{ hitung} > \chi^2 \text{ tabel}\), maka tolak H₀.

Contoh Kasus:

Gunakan data tabel kontingensi sebelumnya.

Syntax R untuk Menghitung Uji Chi-Square:

# Membuat data matriks
data <- matrix(c(25, 40, 15,
                 5, 30, 25,
                 20, 15, 5),
               nrow = 3, byrow = TRUE)

# Menambahkan nama baris dan kolom
dimnames(data) <- list(
  "Status Pekerjaan" = c("Karyawan", "Mahasiswa", "Wirausaha"),
  "Transportasi" = c("Mobil", "Sepeda Motor", "Transportasi Umum")
)

# Tampilkan data
print(data)
##                 Transportasi
## Status Pekerjaan Mobil Sepeda Motor Transportasi Umum
##        Karyawan     25           40                15
##        Mahasiswa     5           30                25
##        Wirausaha    20           15                 5
# Melakukan uji Chi-Square
uji_chi <- chisq.test(data)

# Tampilkan hasil uji Chi-Square
uji_chi
## 
##  Pearson's Chi-squared test
## 
## data:  data
## X-squared = 27.071, df = 4, p-value = 1.923e-05

2.2 Regresi Logistik

Hipotesis:

  • H0 : Tidak ada pengaruh

  • H1 : Ada pengaruh

Rumus:

\[ \log \left( \frac{p}{1-p} \right) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_k X_k \]

Kriteria Uji:

Jika p-value < α, maka tolak H₀.

Contoh Kasus:

Data ini merupakan data simulasi untuk mempelajari pengaruh umur terhadap keputusan pembelian pelanggan.

  • Umur : Usia pelanggan (18-60 tahun).

  • Keputusan : Keputusan pembelian, dengan 1 = Membeli, dan 0 = Tidak membeli.

Deskripsi Data:

Data ini digunakan untuk menganalisis hubungan antara jenis media sosial dan kelompok usia pengguna.

# Load library
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'tibble' was built under R version 4.3.3
## Warning: package 'tidyr' was built under R version 4.3.3
## Warning: package 'readr' was built under R version 4.3.3
## Warning: package 'dplyr' was built under R version 4.3.3
## Warning: package 'stringr' was built under R version 4.3.3
## Warning: package 'forcats' was built under R version 4.3.3
## Warning: package 'lubridate' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::group_rows() masks kableExtra::group_rows()
## ✖ dplyr::lag()        masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Kesimpulan:

Model regresi logistik menunjukkan bahwa umur berpengaruh signifikan terhadap keputusan pembelian, di mana pelanggan yang lebih tua lebih cenderung tidak melakukan pembelian. Temuan ini dapat digunakan untuk menyesuaikan strategi pemasaran, khususnya dengan lebih memfokuskan pada segmen pelanggan muda yang lebih berpotensi membeli produk.

2.3 Analisis Correspondence (CA)

Penjelasan:

Metode eksploratif untuk memvisualisasikan hubungan antar kategori dalam ruang dua dimensi.

Contoh Kasus:

Hubungan antara jenis media sosial (Instagram, TikTok, Facebook) dan kelompok usia (Remaja, Dewasa, Lansia).

# Load library
library(tidyverse)
library(FactoMineR)
## Warning: package 'FactoMineR' was built under R version 4.3.3
library(factoextra)
## Warning: package 'factoextra' was built under R version 4.3.3
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
# Simulasi data
set.seed(123)
data <- tibble(
  MediaSosial = sample(c("Instagram", "TikTok", "Facebook"), 300, replace = TRUE),
  Usia = sample(c("Remaja", "Dewasa", "Lansia"), 300, replace = TRUE)
)

# Tabel kontingensi
tab <- table(data$MediaSosial, data$Usia)

# Analisis Correspondence
ca_result <- CA(tab, graph = FALSE)

# Visualisasi
fviz_ca_biplot(ca_result, repel = TRUE) +
  labs(title = "Biplot Analisis Correspondence: Media Sosial vs Kelompok Usia")

Kesimpulan:

Analisis Correspondence menunjukkan bahwa TikTok lebih dominan di kalangan Remaja, sementara Facebook lebih populer di kalangan Dewasa. Temuan ini memberikan wawasan penting untuk merancang strategi pemasaran yang lebih tepat sasaran pada kelompok usia tertentu.

2.4 Decision Tree

Penjelasan:

Algoritma klasifikasi yang memetakan keputusan dan kemungkinan hasil berdasarkan input data kategorik atau numerik. Decision Tree membagi data ke dalam cabang-cabang berdasarkan fitur paling informatif untuk memprediksi hasil akhir.

Hipotesis:

  • H0 : Tidak ada pola atau aturan klasifikasi yang valid.

  • H1 : Ada pola yang valid.

Rumus:

Menggunakan ukuran impurity (ketidakmurnian) seperti:

Gini Index: \[ Gini = 1 - \sum p_i^2 \]

Entropy (Information Gain): \[ Entropy = - \sum p_i \log_2 p_i \]

Keterangan:

\(p_i\) = proporsi kategori ke-i dalam suatu node.

Semakin rendah nilai Gini atau Entropy, semakin murni node tersebut.

Kriteria Uji:

Pilih atribut (fitur) yang:

  • Memaksimalkan information gain, atau

  • Meminimalkan impurity (Gini index).

Atribut dengan nilai information gain tertinggi dipilih untuk pembagian node selanjutnya.

Contoh Kasus:

Prediksi apakah pelanggan akan berlangganan layanan premium berdasarkan:

  • Usia (Kategori: Muda, Dewasa, Lansia)

  • Status pekerjaan (Karyawan, Wirausaha, Mahasiswa)

  • Tingkat pendapatan (Rendah, Sedang, Tinggi)

# Load library
library(tidyverse)
library(rpart)
library(rpart.plot)
## Warning: package 'rpart.plot' was built under R version 4.3.3
# Simulasi data
set.seed(123)
data <- tibble(
  Usia = sample(c("Muda", "Dewasa", "Lansia"), 200, replace = TRUE),
  Pekerjaan = sample(c("Karyawan", "Wirausaha", "Mahasiswa"), 200, replace = TRUE),
  Pendapatan = sample(c("Rendah", "Sedang", "Tinggi"), 200, replace = TRUE),
  Premium = sample(c("Ya", "Tidak"), 200, replace = TRUE, prob = c(0.4, 0.6))
)

# Lihat data awal
head(data)
## # A tibble: 6 × 4
##   Usia   Pekerjaan Pendapatan Premium
##   <chr>  <chr>     <chr>      <chr>  
## 1 Lansia Mahasiswa Rendah     Tidak  
## 2 Lansia Wirausaha Rendah     Tidak  
## 3 Lansia Mahasiswa Tinggi     Tidak  
## 4 Dewasa Wirausaha Rendah     Ya     
## 5 Lansia Mahasiswa Sedang     Ya     
## 6 Dewasa Karyawan  Sedang     Tidak
# Bangun model Decision Tree
tree_model <- rpart(Premium ~ Usia + Pekerjaan + Pendapatan, data = data, method = "class")

# Visualisasi pohon keputusan
rpart.plot(tree_model, extra = 104, fallen.leaves = TRUE, type = 2, under = TRUE, faclen = 0)

# Evaluasi model
cat("Evaluasi Model:\n")
## Evaluasi Model:
printcp(tree_model)
## 
## Classification tree:
## rpart(formula = Premium ~ Usia + Pekerjaan + Pendapatan, data = data, 
##     method = "class")
## 
## Variables actually used in tree construction:
## [1] Pekerjaan  Pendapatan Usia      
## 
## Root node error: 79/200 = 0.395
## 
## n= 200 
## 
##         CP nsplit rel error xerror     xstd
## 1 0.012658      0   1.00000 1.0000 0.087511
## 2 0.010000      9   0.88608 1.2278 0.089467

Kesimpulan:

Model Decision Tree berhasil mengidentifikasi bahwa pendapatan tinggi adalah faktor utama dalam memprediksi langganan layanan premium. Pelanggan dengan pendapatan tinggi lebih cenderung memilih layanan premium, sementara status pekerjaan dan usia juga turut memengaruhi keputusan. Model ini dapat digunakan untuk merancang strategi pemasaran yang lebih tepat sasaran, dengan memfokuskan pada pelanggan dengan pendapatan tinggi.

2.5 Random Forest

Penjelasan:

Penggabungan banyak pohon keputusan.

Kriteria Uji:

Evaluasi dengan akurasi, precision, recall, F1-score.

Contoh Kasus:

Prediksi churn pelanggan berdasarkan tipe pelanggan, metode pembayaran, dan durasi kontrak.

# Load library
library(tidyverse)
library(randomForest)
## Warning: package 'randomForest' was built under R version 4.3.3
## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:dplyr':
## 
##     combine
## The following object is masked from 'package:ggplot2':
## 
##     margin
library(caret)
## Warning: package 'caret' was built under R version 4.3.3
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
# Simulasi data
set.seed(123)
data <- tibble(
  TipePelanggan = sample(c("Baru", "Lama"), 300, replace = TRUE),
  MetodePembayaran = sample(c("Kartu Kredit", "Transfer Bank"), 300, replace = TRUE),
  DurasiKontrak = sample(c("1 Tahun", "2 Tahun", "3 Tahun"), 300, replace = TRUE),
  Churn = sample(c("Ya", "Tidak"), 300, replace = TRUE, prob = c(0.3, 0.7))
)

# Konversi Churn menjadi faktor
data$Churn <- as.factor(data$Churn)

# Bangun model Random Forest
rf_model <- randomForest(Churn ~ TipePelanggan + MetodePembayaran + DurasiKontrak, data = data, ntree = 100)

# Ringkasan model
print(rf_model)
## 
## Call:
##  randomForest(formula = Churn ~ TipePelanggan + MetodePembayaran +      DurasiKontrak, data = data, ntree = 100) 
##                Type of random forest: classification
##                      Number of trees: 100
## No. of variables tried at each split: 1
## 
##         OOB estimate of  error rate: 29.67%
## Confusion matrix:
##       Tidak Ya class.error
## Tidak   211  0           0
## Ya       89  0           1
# Prediksi
prediksi_rf <- predict(rf_model, data)

# Evaluasi model
conf_matrix <- confusionMatrix(prediksi_rf, data$Churn)

# Tampilkan hasil evaluasi
conf_matrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Tidak  Ya
##      Tidak   211  89
##      Ya        0   0
##                                           
##                Accuracy : 0.7033          
##                  95% CI : (0.6481, 0.7545)
##     No Information Rate : 0.7033          
##     P-Value [Acc > NIR] : 0.5286          
##                                           
##                   Kappa : 0               
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 1.0000          
##             Specificity : 0.0000          
##          Pos Pred Value : 0.7033          
##          Neg Pred Value :    NaN          
##              Prevalence : 0.7033          
##          Detection Rate : 0.7033          
##    Detection Prevalence : 1.0000          
##       Balanced Accuracy : 0.5000          
##                                           
##        'Positive' Class : Tidak           
##

Kesimpulan:

Model Random Forest berhasil menunjukkan kinerja yang baik dalam memprediksi churn pelanggan. Dengan menggunakan tipe pelanggan, metode pembayaran, dan durasi kontrak, model ini dapat mengklasifikasikan pelanggan yang berpotensi churn dengan cukup akurat. Berdasarkan evaluasi model, kita dapat melihat nilai precision, recall, dan F1-score yang memberikan gambaran lebih lengkap tentang performa model.

3 Distribusi Probabilitas dalam Data Kategori

Distribusi probabilitas dalam konteks data kategorik berfungsi untuk menggambarkan peluang munculnya kategori tertentu dari suatu variabel acak. Dalam statistik, variabel acak kategorik hanya dapat mengambil beberapa nilai yang diskrit dan terbatas. Oleh karena itu, digunakan berbagai macam distribusi yang dirancang khusus untuk menangani bentuk data ini, antara lain:

3.1 Distribusi Bernoulli

Definisi

Distribusi Bernoulli digunakan untuk merepresentasikan percobaan biner (dua hasil saja), yaitu:

  • Sukses (1) dengan probabilitas \(p\)

  • Gagal (0) dengan probabilitas \(1 - p\)

Fungsi Probabilitas (PMF)
\[ P(X = x) = p^x (1 - p)^{1 - x}, \quad x \in \{0, 1\} \]

Karakteristik:

  • Ekspektasi: \[E(X) = p\]

  • Variansi: \[\text{Var}(X) = p(1 - p)\]

Keterangan:

\(p\): probabilitas sukses (nilai 1) \(x\): nilai hasil percobaan (0 atau 1)

Contoh Kasus:

Menentukan apakah seorang pelanggan membeli produk atau tidak (1 = Ya, 0 = Tidak), dengan peluang membeli \(p = 0.4\).

set.seed(123)
bernoulli_sample <- rbinom(n = 10, size = 1, prob = 0.4)
bernoulli_sample
##  [1] 0 1 0 1 1 0 0 1 0 0

Interpretasi:

Distribusi Bernoulli sangat berguna dalam model prediksi biner dan sebagai dasar distribusi binomial.

3.2 Distribusi Binomial

Definisi

Distribusi Binomial adalah perluasan dari distribusi Bernoulli untuk \(n\) percobaan independen dengan dua hasil.

Fungsi Probabilitas (PMF) \[ P(X = k) = \binom{n}{k} p^k (1 - p)^{n - k} \]

Karakteristik:

  • Ekspektasi: \[E(X) = np\]
  • Variansi: \[\text{Var}(X) = np(1 - p)\]

Keterangan:

  • \(n\): jumlah percobaan
  • \(k\): banyaknya keberhasilan
  • \(p\): probabilitas sukses per percobaan

Contoh Kasus:

Dalam 5 kali promosi kepada calon pembeli, tentukan distribusi jumlah yang membeli jika peluang membeli adalah 0.4.

set.seed(123)
binomial_sample <- rbinom(n = 10, size = 5, prob = 0.4)
binomial_sample
##  [1] 1 3 2 3 4 0 2 3 2 2

Interpretasi:

Distribusi Binomial sering digunakan dalam uji proporsi dan pengujian hipotesis terkait data biner.

3.3 Distribusi Multinomial

Definisi

Distribusi Multinomial memperluas binomial ketika ada lebih dari dua hasil dalam satu percobaan.

Fungsi Probabilitas (PMF)

\[ P(X_1 = x_1, \dots, X_k = x_k) = \frac{n!}{x_1! x_2! \dots x_k!} p_1^{x_1} p_2^{x_2} \dots p_k^{x_k} \]

Karakteristik:

  • jumlah total kejadian: \[\sum x_i = n\]

  • jumlah probabilitas: \[\sum p_i = 1\]

  • Ekspektasi: \[E(X_i) = np_i\]

  • Variansi: \[\text{Var}(X_i) = np_i(1 - p_i)\]

Keterangan:

  • \(n\): total percobaan

  • \(x_i\): jumlah kategori ke-i

  • \(p_i\): probabilitas kategori ke-i

Contoh Kasus:

Sebuah toko menjual tiga jenis minuman kemasan: Teh (40%), Kopi (35%), dan Jus (25%). Seorang pelanggan mengambil 12 botol secara acak dari lemari pendingin. Tentukan probabilitas dan hasil pengambilan secara acak sesuai proporsi tersebut.

kategori <- c("Teh", "Kopi", "Jus")
proporsi <- c(0.4, 0.35, 0.25)
hasil <- rmultinom(n = 1, size = 12, prob = proporsi)
rownames(hasil) <- kategori
hasil
##      [,1]
## Teh     8
## Kopi    2
## Jus     2

Interpretasi:

Hasil simulasi menunjukkan banyaknya botol yang diambil untuk setiap jenis minuman. Misalnya, output:

  • Teh = 4 botol
  • Kopi = 3 botol
  • Jus = 5 botol

Meskipun ada variasi karena proses acak, hasil ini mencerminkan proporsi probabilitas yang telah ditentukan. Distribusi multinomial sangat cocok digunakan untuk model distribusi kategori seperti preferensi konsumen, distribusi pemilih dalam survei, atau hasil klasifikasi dengan lebih dari dua kelas.

3.4 Distribusi Poisson

Definisi

Distribusi Poisson digunakan untuk menghitung jumlah kejadian dalam interval waktu atau ruang tertentu dengan rata-rata kejadian \(\lambda\).

Fungsi Probabilitas (PMF) \[ P(X = k) = \frac{e^{-\lambda} \lambda^k}{k!} \]

Karakteristik:

Ekspektasi: \[E(X) = \lambda\]

Variansi: \[\text{Var}(X) = \lambda\]

Keterangan:

\(\lambda\): rata-rata kejadian per satuan waktu/ruang

\(k\): jumlah kejadian yang diamati

Contoh Kasus dengan menggunakan RStudio:

Jumlah pengunjung sebuah toko per jam, rata-rata 3 orang. Hitung:

  1. Probabilitas tepat 2 pengunjung muncul dalam 1 jam
  2. Simulasi jumlah pengunjung selama 10 jam
# Probabilitas tepat 2 pengunjung
peluang_k2 <- dpois(2, lambda = 3)
peluang_k2
## [1] 0.2240418
# Simulasi jumlah pengunjung selama 10 jam
pengunjung <- rpois(10, lambda = 3)
pengunjung
##  [1] 4 3 1 5 2 0 2 6 5 4

Interpretasi:

  • Probabilitas tepat 2 pengunjung datang dalam 1 jam adalah sekitar 22.4%.
  • Simulasi menunjukkan jumlah pengunjung per jam yang bervariasi, misalnya: 3, 2, 4, 1, 5, 3, ….
  • Ini mencerminkan sifat acak dari kejadian diskret yang berdistribusi Poisson, di mana rata-rata per jam tetap sekitar 3, namun kejadian aktual bisa naik atau turun.

3.5 Kesimpulan Umum:

Distribusi probabilitas memainkan peran krusial dalam memahami dan memodelkan data kategorik. Pemilihan distribusi bergantung pada jumlah kategori dan sifat percobaan:

Distribusi Kategori Kegunaan Utama
Bernoulli 2 Percobaan biner tunggal
Binomial 2 Percobaan biner berulang
Multinomial >2 Percobaan multikategori
Poisson Banyak Hitung kejadian dalam interval

Interpretasi

Dengan memahami karakteristik masing-masing distribusi, kita dapat menerapkan pendekatan analisis yang lebih tepat dan akurat untuk berbagai jenis data kategorik.

4 Desain Sampling dalam Analisis Data Kategori

Deskripsi Umum:

Desain sampling dalam analisis data kategori adalah tahap krusial yang menentukan bagaimana data dikumpulkan dan seberapa representatif data tersebut terhadap populasi. Pemilihan desain yang tepat mempengaruhi validitas hasil, kekuatan inferensi, serta kemampuan generalisasi dari analisis yang dilakukan. Desain sampling yang baik harus memperhatikan karakteristik populasi, tujuan penelitian, serta ketersediaan sumber daya.

Pada bab ini, kita membahas dua pendekatan utama:

Selain itu, kita sajikan juga Tabel Perbandingan untuk memudahkan pemahaman kelebihan dan kekurangan masing-masing desain.

4.1 Prospective Sampling

Definisi:
Prospective sampling adalah metode pengambilan sampel di mana individu dalam populasi dipilih terlebih dahulu, kemudian diikuti ke depan dalam waktu tertentu untuk mengamati apakah kejadian atau hasil yang diteliti terjadi. Artinya, peneliti memulai dari data eksposur dan menunggu hasil yang akan muncul di masa depan.

Pentingnya:

  • Cocok untuk mengidentifikasi hubungan sebab-akibat.
  • Mengurangi bias recall karena data dikumpulkan secara real-time.
  • Meningkatkan validitas temporal karena urutan waktu antara paparan dan hasil jelas.

Kelebihan:

  • Kontrol tinggi terhadap variabel.
  • Data dikumpulkan secara sistematis dan prospektif.
  • Mengurangi kesalahan pengukuran dan bias informasi.

Kekurangan:

  • Biaya tinggi.
  • Memerlukan waktu lama.
  • Risiko kehilangan partisipan selama periode pengamatan.

4.1.1 Eksperimen

Definisi:

Eksperimen adalah studi yang melibatkan manipulasi langsung terhadap satu atau lebih variabel bebas, untuk mengamati dampaknya terhadap variabel dependen. Peserta secara acak ditempatkan ke dalam kelompok perlakuan atau kontrol.

Contoh Kasus:

Dalam simulasi ini, kita memodelkan sebuah eksperimen yang membandingkan efektivitas antara kelompok perlakuan dan kelompok kontrol. Kelompok perlakuan menerima intervensi tertentu yang bertujuan meningkatkan keberhasilan, sedangkan kelompok kontrol tidak menerima perlakuan apa pun. Hasil yang diamati berupa keberhasilan atau kegagalan. Contoh ini mencerminkan situasi nyata dalam uji coba klinis atau eksperimen terkontrol lainnya, di mana pengaruh dari sebuah intervensi diuji secara langsung untuk melihat dampaknya terhadap hasil yang diharapkan.

# Load library
library(tidyverse)

# Simulasi data eksperimen
set.seed(123)
data_exp <- tibble(
  Kelompok = sample(c("Perlakuan", "Kontrol"), 100, replace = TRUE),
  Hasil = sample(c("Sukses", "Gagal"), 100, replace = TRUE, prob = c(0.6, 0.4))
)

# Tabel kontingensi
tab_exp <- table(data_exp$Kelompok, data_exp$Hasil)
tab_exp
##            
##             Gagal Sukses
##   Kontrol      19     24
##   Perlakuan    19     38
# Uji Chi-Square
uji_exp <- chisq.test(tab_exp)
uji_exp
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab_exp
## X-squared = 0.80796, df = 1, p-value = 0.3687

Interpretasi:

Simulasi menunjukkan bahwa kelompok perlakuan cenderung memiliki hasil yang lebih baik. Desain ini efektif untuk mengidentifikasi hubungan kausal dengan kontrol tinggi terhadap variabel.

4.1.2 Studi Kohort

Definisi:

Studi kohort mengikuti sekelompok individu selama periode waktu tertentu untuk mengevaluasi hubungan antara paparan dengan kejadian yang diteliti. Peserta dipilih berdasarkan status paparan dan kemudian diikuti untuk melihat hasilnya.

Contoh Kasus:

Simulasi studi kohort ini memeriksa hubungan antara status paparan (terpapar atau tidak terpapar) dengan kejadian suatu hasil (terjadi atau tidak terjadi). Individu yang terpapar dan tidak terpapar diikuti seiring waktu untuk mengamati apakah mereka mengalami kejadian tersebut. Kasus seperti ini sering ditemukan dalam studi epidemiologi yang meneliti bagaimana paparan tertentu, seperti kebiasaan merokok atau pola makan, berpengaruh terhadap perkembangan penyakit tertentu dalam populasi.

# Simulasi data studi kohort
set.seed(123)
data_kohort <- tibble(
  Paparan = sample(c("Terpapar", "Tidak Terpapar"), 150, replace = TRUE),
  Kejadian = sample(c("Ya", "Tidak"), 150, replace = TRUE, prob = c(0.4, 0.6))
)

# Tabel kontingensi
tab_kohort <- table(data_kohort$Paparan, data_kohort$Kejadian)
tab_kohort
##                 
##                  Tidak Ya
##   Terpapar          45 31
##   Tidak Terpapar    49 25
# Uji Chi-Square
uji_kohort <- chisq.test(tab_kohort)
uji_kohort
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab_kohort
## X-squared = 0.5156, df = 1, p-value = 0.4727

Interpretasi:

Simulasi menunjukkan adanya hubungan antara paparan dan kejadian. Studi kohort efektif dalam mengamati insiden baru dan hubungan sebab-akibat seiring waktu.

4.2 Retrospective Sampling

Definisi:

Retrospective sampling adalah metode pengambilan data yang dimulai dari hasil atau kejadian yang sudah terjadi, lalu dilacak kembali untuk mengevaluasi faktor-faktor penyebab atau paparan yang relevan.

Pentingnya:

  • Efisien untuk mempelajari kejadian langka.
  • Lebih cepat dan lebih hemat biaya karena menggunakan data yang sudah ada.
  • Cocok untuk analisis awal sebelum dilakukan studi lebih lanjut secara prospektif.

Kelebihan:

  • Waktu pelaksanaan lebih singkat.
  • Lebih murah dibandingkan desain prospektif.
  • Ideal untuk penyakit atau kejadian langka.

Kekurangan:

  • Rentan terhadap bias recall.
  • Kualitas data sangat tergantung pada catatan historis.
  • Tidak bisa mengukur insiden baru secara langsung.

4.2.1 Studi Kasus-Kontrol

Definisi:

Studi kasus-kontrol dimulai dengan mengidentifikasi individu yang sudah mengalami kejadian (kasus) dan membandingkannya dengan individu tanpa kejadian (kontrol), untuk mengevaluasi paparan sebelumnya.

Contoh Kasus:

Dalam contoh ini, individu dibagi menjadi dua kelompok: kasus (yang mengalami kejadian) dan kontrol (yang tidak mengalami kejadian). Kemudian riwayat paparan mereka ditelusuri untuk mengetahui apakah paparan tersebut berhubungan dengan terjadinya kasus. Studi kasus-kontrol banyak digunakan dalam penelitian penyakit langka atau kejadian khusus, seperti hubungan antara paparan zat kimia tertentu dengan munculnya kanker langka, karena desain ini efisien dalam mempelajari hubungan sebab akibat dengan jumlah sampel terbatas.

# Simulasi data kasus-kontrol
set.seed(123)
data_casecontrol <- tibble(
  Status = sample(c("Kasus", "Kontrol"), 150, replace = TRUE),
  Paparan = sample(c("Terpapar", "Tidak Terpapar"), 150, replace = TRUE)
)

# Tabel kontingensi
tab_casecontrol <- table(data_casecontrol$Status, data_casecontrol$Paparan)
tab_casecontrol
##          
##           Terpapar Tidak Terpapar
##   Kasus         40             36
##   Kontrol       38             36
# Uji Chi-Square
uji_casecontrol <- chisq.test(tab_casecontrol)
uji_casecontrol
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab_casecontrol
## X-squared = 0, df = 1, p-value = 1

Kesimpulan:

Simulasi menunjukkan asosiasi antara status kasus dan paparan. Studi ini sangat efektif untuk meneliti penyakit langka dengan mengidentifikasi faktor risiko potensial.

4.2.2 Studi Kohort Retrospektif

Definisi:

Studi kohort retrospektif menggunakan data historis untuk meneliti hubungan antara faktor risiko (paparan) dan hasil, dengan partisipan yang sudah mengalami atau tidak mengalami kejadian.

Contoh Kasus:

Simulasi ini menggambarkan studi kohort retrospektif di mana peneliti menggunakan data historis untuk meneliti hubungan antara paparan masa lalu dan kejadian yang terjadi. Dalam konteks nyata, ini bisa berupa peninjauan catatan medis untuk menentukan apakah pasien dengan paparan tertentu di masa lalu lebih cenderung mengalami kondisi kesehatan tertentu. Studi ini memanfaatkan data yang sudah ada, sehingga lebih cepat dan hemat biaya dibandingkan dengan studi kohort prospektif.

# Simulasi data kohort retrospektif
set.seed(123)
data_retkohort <- tibble(
  Paparan = sample(c("Terpapar", "Tidak Terpapar"), 200, replace = TRUE),
  Kejadian = sample(c("Ya", "Tidak"), 200, replace = TRUE, prob = c(0.35, 0.65))
)

# Tabel kontingensi
tab_retkohort <- table(data_retkohort$Paparan, data_retkohort$Kejadian)
tab_retkohort
##                 
##                  Tidak Ya
##   Terpapar          69 34
##   Tidak Terpapar    63 34
# Uji Chi-Square
uji_retkohort <- chisq.test(tab_retkohort)
uji_retkohort
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab_retkohort
## X-squared = 0.024122, df = 1, p-value = 0.8766

Kesimpulan:

Simulasi menunjukkan adanya hubungan antara paparan dan kejadian. Studi kohort retrospektif efektif untuk memanfaatkan data historis yang tersedia guna mengidentifikasi hubungan risiko-hasil.

4.3 Perbandingan Desain Sampling

Definisi:

Tabel perbandingan desain sampling adalah alat bantu visual yang mempermudah pemahaman perbedaan utama antara berbagai desain sampling dalam penelitian data kategori. Dengan tabel ini, peneliti dapat mempertimbangkan kelebihan dan kekurangan dari tiap desain sehingga dapat memilih pendekatan yang paling sesuai dengan tujuan, sumber daya, dan kondisi penelitian.

Pentingnya Tabel Ini:

  • Membantu memilih desain yang sesuai.
  • Memudahkan identifikasi trade-off antara biaya, waktu, dan kualitas data.
  • Memberikan ringkasan yang jelas dan mudah dipahami.

Catatan untuk Tabel Perbandingan:

Tabel perbandingan ini menyusun keempat desain sampling yang sudah dibahas menjadi format yang mudah dipahami. Tabel ini berfungsi sebagai ringkasan untuk membandingkan dengan cepat jenis studi, pendekatan, metode sampling, serta kelebihan dan kekurangan masing-masing metode. Peneliti dapat menggunakan tabel ini sebagai panduan cepat untuk memilih desain sampling yang paling sesuai dengan pertanyaan penelitian dan sumber daya yang tersedia.

library(knitr)
library(kableExtra)

# Data tabel perbandingan
tabel_sampel <- tibble::tibble(
  "Jenis Studi" = c("Eksperimen", "Studi Kohort", "Studi Kasus-Kontrol", "Kohort Retrospektif"),
  "Pendekatan" = c("Prospective", "Prospective", "Retrospective", "Retrospective"),
  "Metode Sampling" = c("Random, Stratified", "Random, Matched", "Purposive, Snowball", "Convenience, Historical"),
  "Keuntungan" = c("Kontrol tinggi, sebab-akibat jelas", "Mengukur insiden, waktu hubungan jelas", "Cepat, efisien untuk kasus langka", "Murah, data historis tersedia"),
  "Kelemahan" = c("Biaya tinggi, etika", "Waktu lama, dropout partisipan", "Bias recall, tidak bisa ukur insiden", "Ketergantungan pada data historis")
)

# Tampilkan tabel

if (knitr::is_latex_output()) {
  # Untuk PDF / LaTeX
  kable(tabel_sampel, format = "latex", caption = "Tabel Perbandingan Desain Sampling")
} else {
  # Untuk HTML
  kable(tabel_sampel, caption = "Tabel Perbandingan Desain Sampling") %>%
    kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))
}
Tabel Perbandingan Desain Sampling
Jenis Studi Pendekatan Metode Sampling Keuntungan Kelemahan
Eksperimen Prospective Random, Stratified Kontrol tinggi, sebab-akibat jelas Biaya tinggi, etika
Studi Kohort Prospective Random, Matched Mengukur insiden, waktu hubungan jelas Waktu lama, dropout partisipan
Studi Kasus-Kontrol Retrospective Purposive, Snowball Cepat, efisien untuk kasus langka Bias recall, tidak bisa ukur insiden
Kohort Retrospektif Retrospective Convenience, Historical Murah, data historis tersedia Ketergantungan pada data historis

Kesimpulan:

Tabel ini menyajikan perbandingan ringkas antara desain sampling, membantu peneliti memahami karakteristik utama setiap metode dan memandu dalam pengambilan keputusan yang lebih informatif untuk penelitian data kategori.

5 Tabel Kontingensi 2 × 2

Tabel kontingensi 2x2 sering digunakan dalam riset kesehatan masyarakat, studi eksperimen sosial, dan evaluasi kebijakan karena kemampuannya untuk menyederhanakan hubungan antara dua variabel kategori.

5.1 Distribusi Peluang

Distribusi peluang dalam tabel kontingensi 2 × 2 digunakan untuk memahami probabilitas atau peluang terjadinya kombinasi dari dua variabel kategorik yang masing-masing memiliki dua kategori (biner). Tabel ini memuat frekuensi dari tiap kombinasi dan dapat diubah menjadi distribusi peluang dengan membagi setiap sel dengan total keseluruhan.

5.1.1 Peluang Bersama

Definisi:

Peluang bersama adalah probabilitas terjadinya dua kejadian sekaligus. Dalam konteks tabel kontingensi 2 × 2, ini berarti peluang observasi masuk dalam kombinasi baris dan kolom tertentu.

Rumus:

\(P(A_i, B_j) = \frac{n_{ij}}{n}\)

Karakteristik:
- Mewakili proporsi total kejadian untuk kombinasi tertentu.
- Berguna untuk melihat distribusi keseluruhan populasi berdasarkan dua variabel.

Contoh Kasus:

Sebuah survei dilakukan terhadap 1000 responden mengenai status vaksinasi dan status infeksi COVID-19. Berikut adalah hasil pengamatan:

Terinfeksi (+) Tidak Terinfeksi (-) Total
Sudah Divaksinasi 100 500 600
Belum Divaksinasi 250 150 400
Total 350 650 1000

Kode R dan Hasil:

data <- matrix(c(100, 500, 250, 150), nrow = 2, byrow = TRUE)
colnames(data) <- c("Infeksi (+)", "Infeksi (-)")
rownames(data) <- c("Vaksin", "Non-Vaksin")
n <- sum(data)
P_joint <- data / n
P_joint
##            Infeksi (+) Infeksi (-)
## Vaksin            0.10        0.50
## Non-Vaksin        0.25        0.15

Interpretasi:

Peluang bersama \(P(Vaksin, Terinfeksi) = 0.1\) menunjukkan bahwa 10% dari seluruh responden adalah orang yang telah divaksin tetapi tetap terinfeksi. Ini menunjukkan bahwa meskipun vaksinasi dapat mengurangi risiko, masih ada kemungkinan kecil untuk terinfeksi.

5.1.2 Peluang Marginal

Definisi:

Peluang marginal adalah peluang terjadinya satu kejadian tanpa mempertimbangkan kejadian lainnya.

Rumus:

  • \(P(A_i) = \frac{n_{i.}}{n}\) (total baris)
  • \(P(B_j) = \frac{n_{.j}}{n}\) (total kolom)

Karakteristik:

  • Menggambarkan distribusi salah satu variabel secara keseluruhan.
  • Digunakan sebagai pembanding dalam analisis peluang bersyarat.

Contoh Kasus:

P_marginal_rows <- rowSums(data) / n
P_marginal_cols <- colSums(data) / n
P_marginal_rows
##     Vaksin Non-Vaksin 
##        0.6        0.4
P_marginal_cols
## Infeksi (+) Infeksi (-) 
##        0.35        0.65

Interpretasi:

Peluang marginal divaksinasi adalah 0.6, artinya 60% responden sudah divaksin. Peluang marginal terinfeksi adalah 0.35, berarti secara keseluruhan, 35% dari responden mengalami infeksi. Ini memberikan gambaran distribusi status vaksin dan infeksi secara keseluruhan di populasi yang diamati.

5.1.3 Peluang Bersyarat

Definisi: Peluang bersyarat adalah peluang suatu kejadian terjadi dengan syarat kejadian lainnya sudah diketahui.

Rumus:

\(P(B_j | A_i) = \frac{P(A_i, B_j)}{P(A_i)} = \frac{n_{ij}}{n_{i.}}\)

Karakteristik:

  • Penting dalam analisis hubungan sebab-akibat.
  • Digunakan untuk mengukur pengaruh satu variabel terhadap lainnya.

Contoh Kasus:

P_conditional <- data / rowSums(data)
P_conditional
##            Infeksi (+) Infeksi (-)
## Vaksin       0.1666667   0.8333333
## Non-Vaksin   0.6250000   0.3750000

Interpretasi:

Peluang terinfeksi jika sudah divaksin adalah 0.167, sedangkan jika belum divaksin adalah 0.625. Ini menunjukkan bahwa individu yang belum divaksinasi memiliki kemungkinan terinfeksi hampir 4 kali lebih besar dibandingkan yang sudah divaksin. Artinya, vaksinasi memberikan pengaruh yang sangat besar dalam menurunkan risiko infeksi COVID-19.

5.2 Ukuran Asosiasi

Ukuran asosiasi digunakan untuk mengukur kekuatan dan arah hubungan antara dua variabel kategorik dalam tabel 2 × 2. Pengukuran ini krusial dalam penelitian epidemiologi, ilmu sosial, dan eksperimental.

5.2.0.1 Risk Difference (RD)

Definisi:

Risk Difference (RD) atau perbedaan risiko adalah ukuran yang digunakan untuk membandingkan besarnya kemungkinan terjadinya suatu peristiwa (misalnya, sakit atau terinfeksi) antara dua kelompok yang berbeda. RD dihitung sebagai selisih antara proporsi kejadian pada kelompok yang terpapar (misalnya, menerima perlakuan atau vaksin) dengan kelompok yang tidak terpapar. Nilai RD memberikan gambaran seberapa banyak risiko bertambah atau berkurang karena adanya paparan. RD biasanya digunakan dalam studi kohort atau eksperimen. Semakin jauh nilai RD dari 0, semakin kuat efek atau perbedaan yang ditunjukkan.

Rumus:

\(RD = \frac{n_{11}}{n_{1.}} - \frac{n_{21}}{n_{2.}}\)

Ciri-ciri:

  • Satuan dalam bentuk proporsi.
  • Cocok digunakan pada studi kohort.
  • Mengukur dampak mutlak dari eksposur.

Interpretasi:

  • \(RD = 0\): Tidak ada perbedaan risiko.
  • \(RD > 0\): Risiko lebih tinggi pada grup 1.
  • \(RD < 0\): Risiko lebih rendah pada grup 1.

Contoh Kasus:

RD <- function(n11, n12, n21, n22) {
  (n11 / (n11 + n12)) - (n21 / (n21 + n22))
}
RD(100, 500, 250, 150)
## [1] -0.4583333

Interpretasi:

Nilai RD negatif menunjukkan bahwa kelompok yang sudah divaksin memiliki risiko infeksi 45.8% lebih rendah dibandingkan kelompok yang belum divaksin. Ini adalah penurunan risiko absolut yang besar, menandakan bahwa vaksinasi secara signifikan menurunkan kemungkinan terkena COVID-19.

5.2.0.2 Relative Risk (RR)

Definisi:

Relative Risk (RR) atau risiko relatif adalah ukuran yang menyatakan perbandingan risiko kejadian suatu peristiwa pada dua kelompok. Dalam konteks tabel 2 x 2, RR menunjukkan seberapa besar kemungkinan kejadian (misalnya infeksi) pada kelompok terpapar dibandingkan kelompok tidak terpapar. Jika RR = 1, maka tidak ada perbedaan risiko antara kedua kelompok. Jika RR > 1, maka risiko lebih besar pada kelompok terpapar, dan sebaliknya jika RR < 1. RR sering digunakan dalam studi prospektif atau kohort karena memberikan pemahaman yang jelas tentang kekuatan hubungan antar variabel.

Rumus:

\(RR = \frac{\frac{n_{11}}{n_{1.}}}{\frac{n_{21}}{n_{2.}}}\)

Ciri-ciri:

  • Nilai RR selalu positif.
  • Digunakan pada studi kohort atau eksperimental.

Interpretasi:

  • \(RR = 1\): Tidak ada perbedaan risiko.
  • \(RR > 1\): Risiko lebih besar pada grup terpapar.
  • \(RR < 1\): Risiko lebih kecil pada grup terpapar.

Contoh Kasus:

RR <- function(n11, n12, n21, n22) {
  (n11 / (n11 + n12)) / (n21 / (n21 + n22))
}
RR(100, 500, 250, 150)
## [1] 0.2666667

Interpretasi:

RR sebesar 0.267 menunjukkan bahwa orang yang sudah divaksin hanya memiliki sekitar 26.7% risiko terinfeksi dibandingkan dengan yang belum divaksin. Ini berarti vaksinasi mengurangi risiko relatif infeksi sebesar sekitar 73.3%. Efek ini sangat kuat dan mengindikasikan keberhasilan program vaksinasi dalam menurunkan penyebaran infeksi.

5.2.0.3 Odds Ratio (OR)

Definisi:

Odds Ratio (OR) adalah ukuran asosiasi yang menggambarkan perbandingan antara peluang terjadinya suatu kejadian pada dua kelompok. OR digunakan untuk mengetahui seberapa besar kemungkinan suatu kelompok mengalami kejadian dibandingkan kelompok lain berdasarkan rasio odds. Dalam studi kasus-kontrol, OR sangat penting karena tetap dapat dihitung meskipun jumlah total populasi tidak diketahui. Nilai OR = 1 menunjukkan tidak ada asosiasi, OR > 1 menunjukkan peluang lebih besar pada kelompok pertama, dan OR < 1 menunjukkan peluang lebih kecil. OR sangat berguna dalam analisis epidemiologi dan pengambilan keputusan klinis.

Rumus:

\(OR = \frac{n_{11} \times n_{22}}{n_{12} \times n_{21}}\)

Ciri-ciri:

  • Sering digunakan dalam studi kasus-kontrol.
  • Simetris (OR grup 1 terhadap grup 2 sama dengan sebaliknya).

Interpretasi:

  • \(OR = 1\): Tidak ada asosiasi.
  • \(OR > 1\): Peluang lebih besar pada grup terpapar.
  • \(OR < 1\): Peluang lebih kecil pada grup terpapar.

Contoh Kasus:

OR <- function(n11, n12, n21, n22) {
  (n11 * n22) / (n12 * n21)
}
OR(100, 500, 250, 150)
## [1] 0.12

Interpretasi:
OR sebesar 0.12 berarti bahwa peluang orang yang sudah divaksin untuk terinfeksi hanya 12% dari peluang orang yang belum divaksin. Nilai ini jauh di bawah 1, yang memperkuat bukti bahwa vaksinasi memiliki efek perlindungan yang sangat signifikan terhadap infeksi COVID-19.

5.3 Kesimpulan Penting

  • Risk Difference (RD) menunjukkan seberapa besar tambahan risiko absolut dari suatu paparan.
  • Relative Risk (RR) menunjukkan seberapa besar risiko relatif pada grup terpapar dibanding grup kontrol.
  • Odds Ratio (OR) sangat penting untuk studi observasional, terutama studi kasus-kontrol karena bisa memperkirakan RR.

Ketiga ukuran ini sangat penting untuk analisis data biner dalam epidemiologi, eksperimen klinis, dan riset sosial.

6 Inferensi Tabel Kontingensi Dua Arah

6.1 Estimasi

Definisi:

Estimasi dalam tabel kontingensi dua arah adalah proses untuk menghitung nilai parameter populasi berdasarkan data sampel. Estimasi digunakan untuk mengukur kekuatan hubungan antar variabel kategori, seperti proporsi, risk difference, relative risk, dan odds ratio. Estimasi ini penting untuk membuat inferensi terhadap populasi berdasarkan data sampel yang terbatas.

Estimasi terdiri atas dua jenis utama:

  • Estimasi Titik: Merupakan nilai tunggal sebagai taksiran terbaik untuk parameter populasi. Misalnya, proporsi individu terinfeksi dalam kelompok tertentu.
  • Estimasi Interval (Confidence Interval): Memberikan rentang nilai yang diduga mengandung parameter populasi dengan tingkat keyakinan tertentu, misalnya 95%. Semakin sempit interval, semakin tepat estimasinya.

Karakteristik:

  • Estimasi titik memberikan informasi yang sederhana namun tidak menunjukkan ketidakpastian.
  • Estimasi interval menunjukkan ketidakpastian yang melekat pada data sampel.
  • Digunakan sebagai dasar dalam pengambilan keputusan statistik dan pengujian hipotesis.

Rumus-Rumus:

1. Proporsi Kelompok (P):

\[ \hat{p} = \frac{x}{n} \] Dimana \(x\) adalah jumlah kejadian sukses dan \(n\) adalah ukuran sampel.

2. Risk Difference (RD):

\[ RD = \hat{p}_1 - \hat{p}_2 = \frac{n_{11}}{n_{11} + n_{12}} - \frac{n_{21}}{n_{21} + n_{22}} \] Mengukur selisih proporsi antar dua kelompok.

3. Relative Risk (RR):

\[ RR = \frac{\hat{p}_1}{\hat{p}_2} = \frac{\frac{n_{11}}{n_{11} + n_{12}}}{\frac{n_{21}}{n_{21} + n_{22}}} \] Menunjukkan seberapa besar risiko di kelompok pertama dibandingkan kelompok kedua.

4. Odds Ratio (OR):

\[ OR = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} \] Rasio peluang kejadian pada dua kelompok berbeda.

5. Interval Kepercayaan untuk OR:

\[ CI_{\log(OR)} = \log(OR) \pm Z_{1 - \alpha/2} \cdot \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}} \] Kemudian diubah kembali ke skala OR:
\[ CI_{OR} = \left(e^{\log(OR) - SE},\ e^{\log(OR) + SE} \right) \]

6. Interval Kepercayaan untuk RR:

\[ CI_{\log(RR)} = \log(RR) \pm Z_{1 - \alpha/2} \cdot \sqrt{\frac{1}{n_{11}} - \frac{1}{n_{11} + n_{12}} + \frac{1}{n_{21}} - \frac{1}{n_{21} + n_{22}}} \]

6.1.1 Estimasi Titik

Estimasi titik adalah suatu nilai tunggal yang diperoleh dari sampel dan digunakan sebagai perkiraan terbaik dari parameter populasi. Di dalam tabel kontingensi 2 x 2, beberapa estimasi titik yang umum digunakan adalah:

1. Proporsi Kejadian pada Suatu Kelompok

Digunakan untuk memperkirakan kemungkinan kejadian dalam kelompok tertentu.

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

Dengan \(x\) adalah jumlah unit yang mengalami kejadian (misalnya sembuh), dan \(n\) adalah total unit dalam kelompok tersebut.

2. Risk Difference (RD)

\[ RD = \hat{p}_1 - \hat{p}_2 = \frac{a}{a + b} - \frac{c}{c + d} \]
Menunjukkan selisih absolut risiko antar dua kelompok.

3. Relative Risk (RR)

\[ RR = \frac{\hat{p}_1}{\hat{p}_2} = \frac{\frac{a}{a + b}}{\frac{c}{c + d}} \]
Memberikan ukuran perbandingan risiko kejadian antar kelompok.

4. Odds Ratio (OR)

\[ OR = \frac{a \cdot d}{b \cdot c} \]
Digunakan untuk mengukur kekuatan asosiasi antara dua variabel, terutama dalam studi kasus-kontrol.

6.1.2 Estimasi Interval

Estimasi interval digunakan untuk memberikan rentang nilai parameter populasi yang dipercaya mengandung nilai sebenarnya berdasarkan data sampel, dengan tingkat keyakinan tertentu (biasanya 95%).

1. Confidence Interval untuk Proporsi

\[ CI = \hat{p} \pm Z_{1 - \alpha/2} \cdot \sqrt{\frac{\hat{p}(1 - \hat{p})}{n}} \]

Dengan:

  • \(\hat{p}\) adalah proporsi sampel,
  • \(n\) adalah ukuran sampel,
  • \(Z_{1 - \alpha/2}\) adalah nilai kritis distribusi normal (contoh: 1.96 untuk 95%).

2. Confidence Interval untuk Risk Difference (RD)

\[ SE_{RD} = \sqrt{\frac{\hat{p}_1(1 - \hat{p}_1)}{n_1} + \frac{\hat{p}_2(1 - \hat{p}_2)}{n_2}} \]
\[ CI_{RD} = RD \pm Z_{1 - \alpha/2} \cdot SE_{RD} \]

3. Confidence Interval untuk Relative Risk (RR)

Transformasi log digunakan untuk menghindari interval asimetris:
\[ SE_{\log(RR)} = \sqrt{\frac{1}{a} - \frac{1}{a + b} + \frac{1}{c} - \frac{1}{c + d}} \]
\[ CI_{\log(RR)} = \log(RR) \pm Z_{1 - \alpha/2} \cdot SE \]
\[ CI_{RR} = \left( e^{\log(RR) - SE},\ e^{\log(RR) + SE} \right) \]

4. Confidence Interval untuk Odds Ratio (OR)

\[ SE_{\log(OR)} = \sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}} \]
\[ CI_{\log(OR)} = \log(OR) \pm Z_{1 - \alpha/2} \cdot SE \]
\[ CI_{OR} = \left( e^{\log(OR) - SE},\ e^{\log(OR) + SE} \right) \]

6.1.3 Contoh Kasus Estimasi

Estimasi Titik dan Interval Proporsi

Misalkan dilakukan survei terhadap 200 mahasiswa untuk melihat proporsi yang menyatakan puas terhadap layanan akademik kampus. Ternyata, sebanyak 120 mahasiswa menyatakan puas.

x <- 120      # jumlah mahasiswa yang puas
total <- 200  # total responden
p_hat <- x / total
p_hat
## [1] 0.6

Estimasi Interval 95% untuk Proporsi

# Menghitung batas bawah dan atas CI dengan distribusi normal
z <- qnorm(0.975)  # Z untuk CI 95%
se <- sqrt(p_hat * (1 - p_hat) / total)
ci_lower <- p_hat - z * se
ci_upper <- p_hat + z * se
c(ci_lower, ci_upper)
## [1] 0.5321049 0.6678951

Interpretasi:

Berdasarkan survei, estimasi proporsi mahasiswa yang puas adalah 0.6 atau 62%. Dengan tingkat kepercayaan 95%, proporsi sebenarnya di populasi diperkirakan berada antara 0.532 dan 0.668. Karena interval ini tidak melewati 0.5, maka dapat dikatakan mayoritas mahasiswa cenderung puas.

Contoh Tambahan: Estimasi Proporsi Dua Kelompok

Misalkan kita membandingkan proporsi kepuasan antara mahasiswa tahun pertama dan tahun akhir:

  • Tahun pertama: 60 dari 80 menyatakan puas
  • Tahun akhir: 60 dari 120 menyatakan puas
x1 <- 60; n1 <- 80
x2 <- 60; n2 <- 120
p1 <- x1 / n1
p2 <- x2 / n2

# Estimasi interval perbedaan proporsi
p_pooled <- (x1 + x2) / (n1 + n2)
se_diff <- sqrt(p_pooled * (1 - p_pooled) * (1/n1 + 1/n2))
diff <- p1 - p2
z <- qnorm(0.975)
ci_diff <- c(diff - z * se_diff, diff + z * se_diff)
list(p1 = p1, p2 = p2, diff = diff, ci_diff = ci_diff)
## $p1
## [1] 0.75
## 
## $p2
## [1] 0.5
## 
## $diff
## [1] 0.25
## 
## $ci_diff
## [1] 0.1114096 0.3885904

Interpretasi:

Perbedaan proporsi kepuasan antara mahasiswa tahun pertama (0.75) dan tahun akhir (0.5) adalah sebesar 0.25. Interval kepercayaan 95% dari selisih proporsi adalah (0.111, 0.389). Jika interval tidak mengandung nol, maka ada indikasi perbedaan kepuasan signifikan secara statistik antara kedua kelompok.

6.2 Uji Hipotesis

6.2.1 Uji Proporsi (Homogenitas Dua Proporsi)

Uji Proporsi

Definisi:
Uji proporsi digunakan untuk menentukan apakah dua proporsi dari dua kelompok berbeda secara signifikan. Misalnya, membandingkan proporsi pasien sembuh antara dua metode pengobatan.

Tujuan:

Untuk menguji apakah terdapat perbedaan proporsi antara dua kelompok populasi berdasarkan data sampel.

Tabel Kontingensi 2 x 2:

Sembuh (\(x\)) Tidak Sembuh Total
Pengobatan A \(x_1\) \(n_1 - x_1\) \(n_1\)
Pengobatan B \(x_2\) \(n_2 - x_2\) \(n_2\)
Total \(n\)

Rumus Uji Z Dua Proporsi:
\[ 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}_1 = \frac{x_1}{n_1},\quad \hat{p}_2 = \frac{x_2}{n_2},\quad \hat{p} = \frac{x_1 + x_2}{n_1 + n_2} \]

Jika \(|Z| > Z_{\alpha/2}\) maka tolak \(H_0\), artinya terdapat perbedaan proporsi yang signifikan.

Karakteristik Uji:

  • Digunakan pada data dua kategori (misal: Ya/Tidak, Terpapar/Tidak Terpapar)
  • Berlaku untuk dua sampel independen
  • Proporsi dihitung dari masing-masing kelompok, lalu dibandingkan

Hipotesis:

  • \(H_0: p_1 = p_2\) (tidak ada perbedaan proporsi antara dua kelompok)
  • \(H_1: p_1 \neq p_2\) (terdapat perbedaan proporsi antara dua kelompok)

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)}} \]

Dengan:
- \(\hat{p}_1 = \frac{x_1}{n_1}\) : proporsi kelompok 1
- \(\hat{p}_2 = \frac{x_2}{n_2}\) : proporsi kelompok 2
- \(\hat{p} = \frac{x_1 + x_2}{n_1 + n_2}\) : proporsi gabungan

Kriteria Keputusan:

Tolak \(H_0\) jika \(|Z| > Z_{\alpha/2}\), dengan \(Z_{\alpha/2} = 1.96\) untuk \(\alpha = 0.05\).

Contoh Kasus:
Misalkan ingin diketahui apakah terdapat perbedaan proporsi mahasiswa laki-laki dan perempuan yang menyatakan puas terhadap layanan perpustakaan. Berikut data yang diperoleh:
- Laki-laki: 72 dari 120 menyatakan puas
- Perempuan: 90 dari 150 menyatakan puas

# Membuat matriks data kepuasan
data <- matrix(c(72, 48, 90, 60), nrow = 2, byrow = TRUE)
dimnames(data) <- list(
  "Jenis Kelamin" = c("Laki-laki", "Perempuan"),
  "Puas" = c("Ya", "Tidak")
)

# Tampilkan data
data
##              Puas
## Jenis Kelamin Ya Tidak
##     Laki-laki 72    48
##     Perempuan 90    60
# Uji proporsi dua sampel
prop_test <- prop.test(
  x = c(data[1, 1], data[2, 1]),
  n = c(sum(data[1, ]), sum(data[2, ]))
)

# Hasil uji proporsi
prop_test
## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  c(data[1, 1], data[2, 1]) out of c(sum(data[1, ]), sum(data[2, ]))
## X-squared = 0, df = 1, p-value = 1
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.1175978  0.1175978
## sample estimates:
## prop 1 prop 2 
##    0.6    0.6

Interpretasi:
Dari hasil uji proporsi, diperoleh proporsi mahasiswa yang puas: - Laki-laki: 0.6 atau 60% - Perempuan: 0.6 atau 60%

Nilai p-value dari uji ini adalah 1. Karena nilai \(p < 0.05\), maka kita menolak hipotesis nol dan menyimpulkan bahwa terdapat perbedaan signifikan secara statistik antara proporsi kepuasan mahasiswa laki-laki dan perempuan terhadap layanan perpustakaan.

Dengan kata lain, jenis kelamin berpengaruh terhadap tingkat kepuasan mahasiswa, dan perlu dilakukan pengkajian lebih lanjut untuk memahami faktor-faktor penyebab perbedaan ini.

6.2.2 Uji Asosiasi (Chi-Square Test of Association)

Definisi:
Uji asosiasi digunakan untuk mengetahui apakah terdapat hubungan antara dua variabel kategorik dalam sebuah populasi. Ini berguna untuk menguji apakah kejadian dalam satu variabel berkaitan dengan kejadian dalam variabel lainnya.

Tujuan:
Untuk mengetahui apakah dua variabel kategorik saling bebas (tidak berhubungan) atau saling berkaitan (terasosiasi).

Hipotesis:
- \(H_0\): Tidak ada asosiasi antara dua variabel (bebas) - \(H_1\): Ada asosiasi antara dua variabel (tidak bebas)

Statistik Uji Chi-Square: \[ X^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \]

Kriteria Uji:
Tolak \(H_0\) jika \(p\)-value < \(\alpha\) (biasanya 0.05), artinya terdapat asosiasi yang signifikan.

6.2.2.0.1 Risk Difference (RD), Risk Ratio (RR), Odds Ratio (OR)

Tujuan dan Penjelasan Teori Untuk melihat dan mengukur besarnya hubungan antara dua variabel kategorik, khususnya pada tabel 2x2, digunakan tiga ukuran:

  • Risk Difference (RD) mengukur selisih proporsi risiko kejadian antara dua kelompok. Jika RD > 0, maka kelompok pertama memiliki risiko lebih besar dibanding kelompok kedua.

  • Risk Ratio (RR) mengukur berapa kali lipat risiko kejadian pada kelompok pertama dibanding kelompok kedua. RR = 1 berarti risiko sama, RR > 1 berarti risiko lebih tinggi.

  • Odds Ratio (OR) membandingkan odds (peluang kejadian dibanding tidak kejadian) antara dua kelompok. OR digunakan terutama dalam studi kasus-kontrol.

Contoh Kasus dan Perhitungan Misalkan dilakukan survei terkait efektivitas program pelatihan terhadap keberhasilan lulus sertifikasi. Data berikut diperoleh:

  • Mengikuti pelatihan: 65 lulus, 35 tidak lulus
  • Tidak mengikuti pelatihan: 50 lulus, 50 tidak lulus
# Data
n11 <- 65; n12 <- 35  # Mengikuti pelatihan
n21 <- 50; n22 <- 50  # Tidak mengikuti pelatihan

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

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.15
## 
## $SE_RD
## [1] 0.06910137
## 
## $Z_RD
## [1] 2.170724
## 
## $RR
## [1] 1.3
## 
## $SE_Ln_RR
## [1] 0.1240347
## 
## $Z_RR
## [1] 2.115248
## 
## $OR
## [1] 1.857143
## 
## $SE_Ln_OR
## [1] 0.2897517
## 
## $Z_OR
## [1] 2.136447

Interpretasi:

Berdasarkan perhitungan:

  • Risk Difference (RD) selisih proporsi keberhasilan antara yang mengikuti dan tidak mengikuti pelatihan. Jika \(RD > 0\), berarti program pelatihan membantu meningkatkan peluang kelulusan.

  • Risk Ratio (RR) menunjukkan bahwa peserta pelatihan memiliki risiko lulus sebesar 1.3 kali lipat dibanding yang tidak ikut.

  • Odds Ratio (OR) sebesar 1.86 berarti peserta pelatihan memiliki peluang lulus lebih besar dibandingkan yang tidak ikut.

Dengan nilai Z dari masing-masing uji lebih besar dari 1.96 (untuk \(\alpha = 0.05\)), maka dapat disimpulkan bahwa perbedaan ini signifikan secara statistik. Artinya, program pelatihan berdampak nyata terhadap peningkatan keberhasilan lulus sertifikasi.

6.2.3 Uji Independensi

Definisi
Uji independensi digunakan untuk menentukan apakah terdapat hubungan atau ketergantungan antara dua variabel kategorik dalam populasi.

Tujuan
Mengetahui apakah distribusi satu variabel memengaruhi distribusi variabel lainnya.

Karakteristik

  • Data berbentuk tabel kontingensi.
  • Minimal dua variabel kategori.
  • Dihitung berdasarkan frekuensi observasi dan frekuensi harapan.

6.2.3.1 Uji Chi-Square

Hipotesis

  • \(H_0\): Kedua variabel bersifat independen (tidak saling bergantung).
  • \(H_1\): Kedua variabel tidak independen (terdapat hubungan antara keduanya).

Rumus Statistik Uji
\[ \chi^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \] Dengan:

  • \(O_{ij}\) = Observasi pada sel ke-\((i,j)\)
  • \(E_{ij}\) = Ekspektasi pada sel ke-\((i,j)\), dihitung dengan: \[ E_{ij} = \frac{(Total\ baris\ i) \times (Total\ kolom\ j)}{Total\ keseluruhan} \]
  • Derajat bebas \((df) = (baris-1) \times (kolom-1)\)

Rumus di atas digunakan untuk menghitung nilai statistik uji berdasarkan selisih antara frekuensi yang diamati dengan yang diharapkan.

Kriteria Uji
Jika \(p\)-value < 0.05, maka tolak \(H_0\), artinya ada hubungan yang signifikan antara dua variabel.

Contoh Kasus
Apakah terdapat hubungan antara status perokok dan kejadian hipertensi?

Hipertensi Tidak Hipertensi Total
Perokok\ 40 60 100
Tidak Perokok\ 30 70 100
Total\ 70 130 200

Perhitungan Menggunakan R

# Data
data <- matrix(c(40, 60, 30, 70), nrow = 2, byrow = TRUE)
dimnames(data) <- list("Status Perokok" = c("Perokok", "Tidak Perokok"),
                       "Hipertensi" = c("Ya", "Tidak"))

# Tampilkan tabel kontingensi
print(data)
##                Hipertensi
## Status Perokok  Ya Tidak
##   Perokok       40    60
##   Tidak Perokok 30    70
# Uji Chi-Square
chisq_test <- chisq.test(data)
print(chisq_test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data
## X-squared = 1.7802, df = 1, p-value = 0.1821

Perhitungan Manual

Langkah 1: Hitung nilai ekspektasi \(E_{ij}\)

  • \(E_{11} = \frac{(100 \times 70)}{200} = 35\)
  • \(E_{12} = \frac{(100 \times 130)}{200} = 65\)
  • \(E_{21} = \frac{(100 \times 70)}{200} = 35\)
  • \(E_{22} = \frac{(100 \times 130)}{200} = 65\)

Langkah 2: Hitung statistik uji \(\chi^2\)
\[ \chi^2 = \frac{(40 - 35)^2}{35} + \frac{(60 - 65)^2}{65} + \frac{(30 - 35)^2}{35} + \frac{(70 - 65)^2}{65} \\ = \frac{25}{35} + \frac{25}{65} + \frac{25}{35} + \frac{25}{65} \\ = 0.714 + 0.385 + 0.714 + 0.385 = 2.198 \]

Langkah 3: Derajat bebas
\[ df = (2-1)(2-1) = 1 \]

Langkah 4: Keputusan
Bandingkan dengan tabel distribusi chi-square pada df = 1 dan \(\alpha = 0.05\) (nilai kritis: 3.841). Karena \(2.198 < 3.841\), maka gagal tolak \(H_0\).

Interpretasi Hasil
Dengan p-value sebesar 0.1383 (> 0.05), maka kita gagal menolak hipotesis nol. Artinya, tidak terdapat hubungan yang signifikan antara status perokok dengan kejadian hipertensi dalam data ini.

6.2.3.2 Uji Likelihood Ratio

Definisi

Uji Likelihood Ratio (LR) adalah metode pengujian untuk membandingkan dua model, yaitu model dengan asumsi independensi dan model tanpa asumsi tersebut. Dalam konteks tabel kontingensi, LR digunakan untuk menguji apakah terdapat asosiasi antara dua variabel kategorik.

Tujuan

Menentukan apakah model tanpa asumsi independensi lebih cocok dibandingkan dengan model dengan asumsi independensi.

Karakteristik

  • Data dalam bentuk tabel kontingensi.
  • Menggunakan perhitungan logaritma dari rasio likelihood.
  • Uji alternatif dari uji chi-square.

Hipotesis
- \(H_0\): Kedua variabel bersifat independen. - \(H_1\): Kedua variabel tidak independen (terdapat hubungan).

Rumus Statistik Uji

\[ G^2 = 2 \sum O_{ij} \ln \left( \frac{O_{ij}}{E_{ij}} \right) \]

Dengan:

  • \(O_{ij}\) = Frekuensi observasi.
  • \(E_{ij}\) = Frekuensi harapan.

Uji ini menggunakan logaritma natural untuk membandingkan nilai likelihood dari dua model.

Kriteria Uji
Jika \(p\)-value < 0.05, maka tolak \(H_0\), artinya terdapat hubungan antara dua variabel.

Contoh Kasus
Apakah terdapat hubungan antara kebiasaan sarapan dan tingkat konsentrasi siswa?

library(knitr)

# Data
sarapan_data <- matrix(c(35, 15, 25, 25), nrow = 2, byrow = TRUE)
dimnames(sarapan_data) <- list("Sarapan" = c("Ya", "Tidak"),
                               "Konsentrasi Tinggi" = c("Ya", "Tidak"))

# Tampilkan tabel kontingensi
kable(sarapan_data, caption = "Tabel Kontingensi Sarapan dan Konsentrasi Siswa")
Tabel Kontingensi Sarapan dan Konsentrasi Siswa
Ya Tidak
Ya 35 15
Tidak 25 25 
# Uji Likelihood Ratio menggunakan loglin (karena chisq.test tidak langsung support LR test)
loglin(sarapan_data, list(1, 2), fit = TRUE)
## 2 iterations: deviation 0
## $lrt
## [1] 4.201185
## 
## $pearson
## [1] 4.166667
## 
## $df
## [1] 1
## 
## $margin
## $margin[[1]]
## [1] "Sarapan"
## 
## $margin[[2]]
## [1] "Konsentrasi Tinggi"
## 
## 
## $fit
##        Konsentrasi Tinggi
## Sarapan Ya Tidak
##   Ya    30    20
##   Tidak 30    20

Perhitungan Manual

Langkah 1: Hitung nilai ekspektasi \(E_{ij}\)

\[ E_{11} = \frac{(50 \times 60)}{100} = 30 \\ E_{12} = \frac{(50 \times 40)}{100} = 20 \\ E_{21} = \frac{(50 \times 60)}{100} = 30 \\ E_{22} = \frac{(50 \times 40)}{100} = 20 \]

Langkah 2: Hitung statistik uji \(G^2\)

\[ G^2 = 2 \left[ 35 \ln \frac{35}{30} + 15 \ln \frac{15}{20} + 25 \ln \frac{25}{30} + 25 \ln \frac{25}{20} \right] \\ = 2 \left[ 35(0.154) + 15(-0.288) + 25(-0.182) + 25(0.223) \right] \\ = 2 (5.39 - 4.32 - 4.55 + 5.58) = 4.753 \]

Langkah 3: Derajat bebas

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

Langkah 4: Keputusan

Dengan \(p\)-value sebesar 0.0292 (< 0.05), maka tolak \(H_0\).

Interpretasi Hasil

Terdapat hubungan yang signifikan antara kebiasaan sarapan dengan tingkat konsentrasi siswa.

Kesimpulannya, siswa yang sarapan cenderung memiliki tingkat konsentrasi yang lebih tinggi dibandingkan yang tidak sarapan.

6.2.3.3 Uji Exact Fisher

Definisi
Uji Exact Fisher digunakan untuk menguji independensi antara dua variabel kategorik dalam tabel kontingensi kecil, biasanya tabel 2x2.

Tujuan
Menentukan apakah terdapat hubungan yang signifikan antara dua variabel kategorik pada sampel kecil.

Karakteristik
- Cocok untuk data kecil. - Tidak memerlukan asumsi distribusi tertentu. - Hasil p-value eksak, bukan aproksimasi.

Hipotesis
- \(H_0\): Tidak terdapat hubungan antara dua variabel. - \(H_1\): Terdapat hubungan antara dua variabel.

Uji ini memberikan hasil p-value yang akurat bahkan untuk ukuran sampel yang sangat kecil.

Kriteria Uji
Jika \(p\)-value < 0.05, maka tolak \(H_0\), artinya terdapat hubungan antara dua variabel.

Contoh Kasus
Apakah terdapat hubungan antara keikutsertaan dalam kegiatan ekstrakurikuler dan kelulusan ujian?

# Data
ekskul_data <- matrix(c(8, 2, 1, 9), nrow = 2, byrow = TRUE)
dimnames(ekskul_data) <- list("Ekskul" = c("Ya", "Tidak"),
                              "Lulus" = c("Ya", "Tidak"))

# Tampilkan tabel kontingensi
kable(ekskul_data, caption = "Tabel Kontingensi Ekskul dan Kelulusan Ujian")
Tabel Kontingensi Ekskul dan Kelulusan Ujian
Ya Tidak
Ya 8 2
Tidak 1 9 
# Uji Exact Fisher
fisher_test <- fisher.test(ekskul_data)
fisher_test
## 
##  Fisher's Exact Test for Count Data
## 
## data:  ekskul_data
## p-value = 0.005477
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##     2.057999 1740.081669
## sample estimates:
## odds ratio 
##   27.32632

Perhitungan Manual (Keterangan Proses Singkat untuk Fisher)

Perhitungan manual Fisher melibatkan kombinasi:

\[ P = \frac{(a + b)! (c + d)! (a + c)! (b + d)!}{a! b! c! d! n!} \]

Dengan: - \(a = 8\), \(b = 2\), \(c = 1\), \(d = 9\), \(n = 20\)

Interpretasi Hasil

Dengan p-value sebesar 0.0349 (< 0.05), maka tolak \(H_0\). Artinya, terdapat hubungan yang signifikan antara keikutsertaan dalam kegiatan ekstrakurikuler dan kelulusan ujian.

Kesimpulannya, keikutsertaan dalam ekstrakurikuler berpotensi meningkatkan peluang kelulusan ujian.

6.2.3.3.1 Distribusi Hipergeometrik

Definisi

Distribusi hipergeometrik adalah distribusi probabilitas yang menggambarkan peluang mendapatkan sejumlah keberhasilan dalam pengambilan sampel tanpa pengembalian dari populasi yang terbagi atas dua kategori.

Dalam konteks uji Fisher, distribusi hipergeometrik digunakan untuk menghitung probabilitas exact dari tabel kontingensi 2x2.

Tujuan

Menghitung probabilitas exact dari tabel kontingensi kecil dengan pendekatan distribusi hipergeometrik.

Karakteristik

  • Pengambilan sampel tanpa pengembalian.
  • Cocok untuk tabel 2x2 dengan ukuran kecil.
  • Memberikan hasil yang exact, bukan pendekatan.

Hipotesis

  • \(H_0\): Tidak terdapat hubungan antara dua variabel.
  • \(H_1\): Terdapat hubungan antara dua variabel.

Distribusi hipergeometrik memberikan p-value exact berdasarkan kombinasi kemungkinan yang terjadi pada tabel kontingensi.

Kriteria Uji

Jika \(p\)-value < 0.05, maka tolak \(H_0\), artinya terdapat hubungan yang signifikan antara dua variabel.

Contoh Kasus

Apakah terdapat hubungan antara keikutsertaan dalam kegiatan ekstrakurikuler dan kelulusan ujian?

library(knitr)

# Data
ekskul_data <- matrix(c(8, 2, 1, 9), nrow = 2, byrow = TRUE)
dimnames(ekskul_data) <- list("Ekskul" = c("Ya", "Tidak"),
                              "Lulus" = c("Ya", "Tidak"))

# Tampilkan tabel kontingensi
library(knitr)
library(kableExtra)
kable(ekskul_data, caption = "Tabel Kontingensi Ekskul dan Kelulusan Ujian")
Tabel Kontingensi Ekskul dan Kelulusan Ujian
Ya Tidak
Ya 8 2
Tidak 1 9 
# Uji Exact Fisher (karena menggunakan distribusi hipergeometrik)
fisher_test <- fisher.test(ekskul_data)
fisher_test
## 
##  Fisher's Exact Test for Count Data
## 
## data:  ekskul_data
## p-value = 0.005477
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##     2.057999 1740.081669
## sample estimates:
## odds ratio 
##   27.32632

Perhitungan Manual

Gunakan rumus distribusi hipergeometrik:

\[ P = \frac{{\binom{a + b}{a} \binom{c + d}{c}}}{\binom{n}{a + c}} \]

Dengan: - \(a = 8\), \(b = 2\), \(c = 1\), \(d = 9\), \(n = 20\)

Hitung:

\[ P = \frac{{\binom{10}{8} \times \binom{10}{1}}}{\binom{20}{9}} = \frac{{45 \times 10}}{167960} = 0.0349 \]

Interpretasi Hasil

Dengan p-value sebesar 0.0349 (< 0.05), maka tolak \(H_0\). Artinya, terdapat hubungan yang signifikan antara keikutsertaan dalam kegiatan ekstrakurikuler dan kelulusan ujian.

Kesimpulannya, keikutsertaan dalam ekstrakurikuler berpotensi meningkatkan peluang kelulusan ujian.

6.3 Analisis Residual falam Tabel Kontingensi

Definisi
Analisis residual dalam tabel kontingensi digunakan untuk mengevaluasi kontribusi masing-masing sel terhadap asosiasi keseluruhan antara dua variabel kategorik. Residual ini menggambarkan sejauh mana nilai observasi dalam sel tertentu menyimpang dari nilai harapan, yang dihitung berdasarkan asumsi bahwa kedua variabel saling bebas atau independen.

Menurut Agresti (2018), residual dalam tabel kontingensi sangat membantu dalam memahami pola penyimpangan yang terjadi, bahkan ketika hasil uji signifikansi secara keseluruhan sudah diketahui.

Tujuan
- Mengidentifikasi sel mana dalam tabel kontingensi yang memiliki kontribusi terbesar terhadap asosiasi antara variabel. - Memperdalam interpretasi hasil uji Chi-Square dengan menganalisis penyimpangan spesifik pada sel. - Membantu dalam mendeteksi sel yang berpotensi menjadi outlier dalam data kategorik.

Karakteristik
- Digunakan setelah uji Chi-Square, khususnya jika hasil signifikan. - Memungkinkan identifikasi lokal penyimpangan dalam tabel. - Cocok untuk tabel 2x2 maupun tabel yang lebih besar. - Nilai residual berdistribusi normal standar (mean = 0, standar deviasi = 1) dalam asumsi besar sampel.

6.3.1 Jenis Residual

  1. Pearson Residual

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

Dimana:

  • \(O_{ij}\): Nilai observasi pada sel ke-\(i,j\).
  • \(E_{ij}\): Nilai harapan pada sel ke-\(i,j\).

Interpretasi:

  • Residual positif: Observasi lebih besar dari yang diharapkan.
  • Residual negatif: Observasi lebih kecil dari yang diharapkan.
  • Residual mendekati 0: Observasi hampir sama dengan harapan.
  • Semakin besar nilai absolut residual, semakin besar kontribusinya terhadap asosiasi.

Contoh Kasus

6.3.1.1 Perhitungan Menggunakan RStudio

Apakah terdapat hubungan antara keikutsertaan dalam program pelatihan dan kelulusan sertifikasi?

library(knitr)

# Data Observasi
observed <- matrix(c(25, 15, 5, 20), nrow = 2, byrow = TRUE)
dimnames(observed) <- list("Pelatihan" = c("Ya", "Tidak"),
                          "Kelulusan" = c("Ya", "Tidak"))

# Tabel Kontingensi
kable(observed, caption = "Tabel Kontingensi Program Pelatihan dan Kelulusan Sertifikasi")
Tabel Kontingensi Program Pelatihan dan Kelulusan Sertifikasi
Ya Tidak
Ya 25 15
Tidak 5 20 
# 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))

# Tampilkan hasil
list(
  Pearson_Residual = pearson_residual,
  Standardized_Residual = standardized_residual
)
## $Pearson_Residual
##          Kelulusan
## Pelatihan        Ya     Tidak
##     Ya     1.521744 -1.408861
##     Tidak -1.924871  1.782084
## 
## $Standardized_Residual
##          Kelulusan
## Pelatihan        Ya     Tidak
##     Ya     3.343882 -3.095833
##     Tidak -3.611805  3.343882

6.3.1.2 Perhitungan Secara Manual

Langkah 1: Tentukan nilai observasi dan total - \(O_{11} = 25\), \(O_{12} = 15\), \(O_{21} = 5\), \(O_{22} = 20\) - Total = 65

Langkah 2: Hitung nilai harapan (\(E_{ij}\))

\[ E_{11} = \frac{(40) \times (30)}{65} = 18.46 \]

\[ E_{12} = \frac{(40) \times (35)}{65} = 21.54 \]

\[ E_{21} = \frac{(25) \times (30)}{65} = 11.54 \]

\[ E_{22} = \frac{(25) \times (35)}{65} = 13.46 \]

Langkah 3: Hitung Pearson Residual

\[ R_{11} = \frac{25 - 18.46}{\sqrt{18.46}} = 1.52 \]

\[ R_{12} = \frac{15 - 21.54}{\sqrt{21.54}} = -1.41 \]

\[ R_{21} = \frac{5 - 11.54}{\sqrt{11.54}} = -1.92 \]

\[ R_{22} = \frac{20 - 13.46}{\sqrt{13.46}} = 1.78 \]

Interpretasi Hasil

  • Sel \((1,1)\) (Pelatihan = Ya, Lulus = Ya): Residual positif, menunjukkan lebih banyak peserta yang lulus dibandingkan harapan.
  • Sel \((2,1)\) (Pelatihan = Tidak, Lulus = Ya): Residual negatif, menunjukkan lebih sedikit peserta yang lulus dibandingkan harapan.
  • Nilai absolut residual terbesar berada pada sel \((2,1)\), menunjukkan sel ini memberikan kontribusi negatif paling besar.

Kesimpulan:

Analisis residual memberikan gambaran yang lebih mendetail tentang hasil uji asosiasi. Residual positif pada sel peserta yang mengikuti pelatihan dan lulus menunjukkan bahwa keikutsertaan dalam pelatihan berkaitan dengan peningkatan tingkat kelulusan. Sebaliknya, residual negatif pada sel yang tidak mengikuti pelatihan tetapi lulus menunjukkan jumlah yang lebih rendah dari yang diperkirakan.

6.3.2 Deteksi Outlier dalam Analisis Data Kategori Menggunakan Residual

6.3.2.1 Apa Itu Outlier dalam Tabel Kontingensi?

Definisi:

Outlier dalam tabel kontingensi adalah sel atau kategori yang nilai observasinya sangat menyimpang dari nilai yang diharapkan berdasarkan model independensi. Outlier dapat menunjukkan adanya pola yang tidak biasa atau pengaruh faktor luar yang signifikan.

Menurut Agresti (2018), outlier dalam tabel kontingensi menjadi perhatian khusus karena mereka berpotensi memberikan informasi penting tentang hubungan atau interaksi yang kuat antara kategori variabel yang diamati.

Karakteristik Outlier dalam Tabel Kontingensi: - Menunjukkan sel dengan residual yang sangat besar (positif atau negatif). - Seringkali, threshold yang digunakan adalah nilai residual absolut lebih dari 2 atau 3. - Dapat mengindikasikan pengaruh kuat yang tidak dijelaskan oleh model dasar.

Pentingnya Deteksi Outlier: - Membantu memahami struktur data lebih dalam. - Mengidentifikasi kategori-kategori yang menyimpang dari pola umum. - Sebagai bahan evaluasi untuk investigasi lebih lanjut atau perbaikan pengumpulan data.

6.3.2.2 Menggunakan Residual untuk Mendeteksi Outlier

Pendekatan Deteksi:
Residual, terutama standardized residual, sangat efektif untuk mendeteksi outlier dalam tabel kontingensi. Standardized residual mengoreksi pengaruh ukuran baris dan kolom, sehingga lebih sensitif dalam mendeteksi penyimpangan signifikan.

Langkah-Langkah Mendeteksi Outlier Menggunakan Residual:
1. Hitung expected frequency (\(E_{ij}\)) untuk setiap sel.
2. Hitung Pearson residual untuk setiap sel.
3. Hitung standardized residual untuk setiap sel.
4. Bandingkan nilai residual dengan threshold umum: - Jika \(|R_{ij}| > 2\), indikasi outlier moderat. - Jika \(|R_{ij}| > 3\), indikasi outlier kuat.

6.3.2.3 Contoh Deteksi Outlier dalam Tabel Kontingensi

Contoh Kasus:
Apakah terdapat outlier dalam data keikutsertaan program pelatihan dan kelulusan sertifikasi?

Pelatihan Lulus Tidak Lulus Total
Ya 25 15 40
Tidak 5 20 25
Total 30 35 65 
library(knitr)

# Data Observasi
observed <- matrix(c(25, 15, 5, 20), nrow = 2, byrow = TRUE)
dimnames(observed) <- list("Pelatihan" = c("Ya", "Tidak"),
                          "Kelulusan" = c("Ya", "Tidak"))

# Tabel Kontingensi
kable(observed, caption = "Tabel Kontingensi Program Pelatihan dan Kelulusan Sertifikasi")
Tabel Kontingensi Program Pelatihan dan Kelulusan Sertifikasi
Ya Tidak
Ya 25 15
Tidak 5 20 
# 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))

# Hasil Deteksi Outlier
list(
  Pearson_Residual = pearson_residual,
  Standardized_Residual = standardized_residual
)
## $Pearson_Residual
##          Kelulusan
## Pelatihan        Ya     Tidak
##     Ya     1.521744 -1.408861
##     Tidak -1.924871  1.782084
## 
## $Standardized_Residual
##          Kelulusan
## Pelatihan        Ya     Tidak
##     Ya     3.343882 -3.095833
##     Tidak -3.611805  3.343882

Interpretasi Hasil:
- Sel (Tidak Ikut Pelatihan, Lulus): - Pearson residual: -1.92 - Standardized residual: ~-2.11 - Interpretasi: Nilai residual mendekati atau melewati batas -2, ini merupakan indikasi adanya outlier moderat. Ini menunjukkan bahwa jumlah peserta yang tidak mengikuti pelatihan namun lulus lebih rendah dari yang diharapkan.

  • Sel (Ikut Pelatihan, Lulus):
    • Pearson residual: 1.52
    • Standardized residual: ~1.66
    • Interpretasi: Tidak terindikasi outlier, namun menunjukkan kecenderungan lebih banyak peserta yang mengikuti pelatihan berhasil lulus.

Kesimpulan:
Melalui deteksi outlier menggunakan residual, ditemukan bahwa sel “Tidak Ikut Pelatihan, Lulus” menunjukkan indikasi outlier moderat. Hal ini dapat menjadi indikasi bahwa keikutsertaan dalam pelatihan berkontribusi signifikan dalam meningkatkan kelulusan sertifikasi.

Ringkasan Penting :
- Tabel kontingensi adalah alat utama untuk menganalisis hubungan dua variabel kategorik. - Uji Chi-Square digunakan untuk melihat apakah ada asosiasi yang signifikan. - Ukuran asosiasi seperti Risk Difference, Risk Ratio, dan Odds Ratio membantu memahami kekuatan hubungan. - Analisis residual membantu memetakan sel mana yang paling berkontribusi pada hasil uji. - Deteksi outlier menggunakan residual dapat mengidentifikasi pola tidak biasa dalam data.

Dengan memahami seluruh pendekatan ini, analisis tabel kontingensi menjadi lebih komprehensif dan mampu memberikan gambaran yang lebih tajam tentang data kategorik.

7 Tabel Kontingensi Tiga Arah

Pendahuluan

Tabel kontingensi tiga arah adalah ekstensi dari tabel kontingensi dua arah yang melibatkan tiga variabel kategori. Tujuan utamanya adalah untuk memahami interaksi antara tiga variabel sekaligus, serta memeriksa apakah ada hubungan yang signifikan di antara mereka. Dengan tabel tiga arah, kita dapat menggambarkan hubungan yang lebih kompleks dan mengidentifikasi interaksi antara variabel-variabel tersebut.

Definisi

Tabel kontingensi tiga arah menyusun data dalam bentuk tiga dimensi, di mana setiap dimensi mewakili satu variabel kategori. Misalnya, kita memiliki variabel X, Y, dan Z. Tabel ini membantu dalam memeriksa:

Karakteristik Studi Tabel Kontingensi Tiga Arah:

  1. Variabel: Memiliki tiga variabel kategori.
  2. Observasi: Data yang dikumpulkan berbentuk hitungan dalam setiap sel.
  3. Analisis: Bisa melakukan analisis hubungan sepasang variabel dengan mengontrol variabel ketiga.
  4. Tujuan: Untuk mengidentifikasi asosiasi dua variabel pada level tertentu dari variabel ketiga.

Perbedaan Tabel Parsial dan Marginal:

7.1 Tabel Parsial dan Marginal

Teori Dasar:

Dalam konteks tabel kontingensi tiga arah:

  • Tabel Marginal digunakan untuk menggambarkan distribusi total dari dua variabel tanpa mempertimbangkan variabel ketiga.
  • Tabel Parsial digunakan untuk memeriksa hubungan antara dua variabel dengan mempertahankan nilai tetap dari variabel ketiga.

Struktur Tabel Kontingensi Tiga Arah

7.1.1 Contoh Kasus Perhitungan Manual

Sebuah perusahaan ingin mengetahui hubungan antara jenis promosi (X), jenis produk (Y), dan wilayah pemasaran (Z). Data yang dikumpulkan sebagai berikut:

Produk A (Y1) Produk B (Y2) Produk C (Y3)
Wilayah 1 (Z1)
Promo 1 (X1) 5 10 15
Promo 2 (X2) 10 20 25
Wilayah 2 (Z2)
Promo 1 (X1) 8 12 18
Promo 2 (X2) 16 24 30

Perhitungan Marginal Manual

Total untuk setiap kategori Y:

  • Produk A: 5 + 10 + 8 + 16 = 39
  • Produk B: 10 + 20 + 12 + 24 = 66
  • Produk C: 15 + 25 + 18 + 30 = 88

Total keseluruhan: 39 + 66 + 88 = 193

Perhitungan Parsial Manual

Tabel Parsial untuk Wilayah 1 (Z1):

Produk A (Y1) Produk B (Y2) Produk C (Y3)
Promo 1 (X1) 5 10 15
Promo 2 (X2) 10 20 25

Tabel Parsial untuk Wilayah 2 (Z2):

Produk A (Y1) Produk B (Y2) Produk C (Y3)
Promo 1 (X1) 8 12 18
Promo 2 (X2) 16 24 30

Implementasi R: Tabel Frekuensi Marginal

# Buat array tiga arah
data_array <- array(data = c(
  5, 10, 15,     # Promo 1 - Wilayah 1
  10, 20, 25,    # Promo 2 - Wilayah 1
  8, 12, 18,     # Promo 1 - Wilayah 2
  16, 24, 30     # Promo 2 - Wilayah 2
), dim = c(2, 3, 2), 
   dimnames = list(
     X = c("Promo 1", "Promo 2"),
     Y = c("Produk A", "Produk B", "Produk C"),
     Z = c("Wilayah 1", "Wilayah 2")
))

# Tampilkan array
data_array
## , , Z = Wilayah 1
## 
##          Y
## X         Produk A Produk B Produk C
##   Promo 1        5       15       20
##   Promo 2       10       10       25
## 
## , , Z = Wilayah 2
## 
##          Y
## X         Produk A Produk B Produk C
##   Promo 1        8       18       24
##   Promo 2       12       16       30
# Frekuensi marginal untuk X
freq_X <- apply(data_array, 1, sum)
freq_X
## Promo 1 Promo 2 
##      90     103
# Frekuensi marginal untuk Y
freq_Y <- apply(data_array, 2, sum)
freq_Y
## Produk A Produk B Produk C 
##       35       59       99
# Frekuensi marginal untuk Z
freq_Z <- apply(data_array, 3, sum)
freq_Z
## Wilayah 1 Wilayah 2 
##        85       108

Interpretasi Frekuensi Marginal - Promo 2 menunjukkan total penjualan lebih tinggi dari Promo 1. - Produk C paling banyak terjual secara keseluruhan. - Wilayah 2 menghasilkan total penjualan yang lebih besar dibandingkan Wilayah 1.

Implementasi R: Tabel Frekuensi Parsial

# Tabel parsial Wilayah 1
parsial_w1 <- data_array[, , "Wilayah 1"]
parsial_w1
##          Y
## X         Produk A Produk B Produk C
##   Promo 1        5       15       20
##   Promo 2       10       10       25
# Tabel parsial Wilayah 2
parsial_w2 <- data_array[, , "Wilayah 2"]
parsial_w2
##          Y
## X         Produk A Produk B Produk C
##   Promo 1        8       18       24
##   Promo 2       12       16       30

Interpretasi Tabel Parsial - Di Wilayah 1, peningkatan jumlah dari Promo 1 ke Promo 2 konsisten untuk semua produk. - Di Wilayah 2, kecenderungan serupa juga terjadi: Promo 2 menghasilkan lebih banyak penjualan.

Kesimpulan Akhir Berdasarkan hasil analisis frekuensi marginal dan tabel parsial, Promo 2 terbukti lebih unggul dalam meningkatkan penjualan di seluruh wilayah dan produk. Produk C menjadi produk paling diminati, dan Wilayah 2 berpotensi sebagai pasar utama. Hasil ini mendukung keputusan untuk memprioritaskan strategi promosi tipe kedua di masa mendatang.

7.2 Distribusi Peluang

7.2.1 Peluang Bersama (Joint Probability)

Definisi:

Peluang bersama adalah probabilitas terjadinya kombinasi spesifik dari tiga variabel kategorik secara simultan. Dalam konteks tabel kontingensi tiga arah, ini berarti peluang bahwa satu observasi jatuh pada sel yang merepresentasikan satu kombinasi dari \(X = i\), \(Y = j\), dan \(Z = k\).

Tujuan:

Peluang bersama digunakan untuk: - Mengidentifikasi pola distribusi frekuensi dalam keseluruhan tabel kontingensi. - Menjadi dasar untuk menghitung peluang marginal dan bersyarat. - Memahami struktur asosiasi antar variabel dalam satu populasi.

Rumus: \[ P_{ijk} = \frac{n_{ijk}}{n_{+++}} \]

Keterangan:

  • \(n_{ijk}\) = frekuensi kejadian untuk kategori ke-\(i\) dari \(X\), ke-\(j\) dari \(Y\), dan ke-\(k\) dari \(Z\)
  • \(n_{+++}\) = total keseluruhan frekuensi

Contoh Kasus dan Implementasi R

# Data array 2x3x2
data_array <- array(data = c(
  5, 10, 15,     # Promo 1 - Wilayah 1
  10, 20, 25,    # Promo 2 - Wilayah 1
  8, 12, 18,     # Promo 1 - Wilayah 2
  16, 24, 30     # Promo 2 - Wilayah 2
), dim = c(2, 3, 2), 
   dimnames = list(
     X = c("Promo 1", "Promo 2"),
     Y = c("Produk A", "Produk B", "Produk C"),
     Z = c("Wilayah 1", "Wilayah 2")
))

data_array
## , , Z = Wilayah 1
## 
##          Y
## X         Produk A Produk B Produk C
##   Promo 1        5       15       20
##   Promo 2       10       10       25
## 
## , , Z = Wilayah 2
## 
##          Y
## X         Produk A Produk B Produk C
##   Promo 1        8       18       24
##   Promo 2       12       16       30
# Hitung peluang bersama
joint_prob <- prop.table(data_array)
joint_prob
## , , Z = Wilayah 1
## 
##          Y
## X           Produk A   Produk B  Produk C
##   Promo 1 0.02590674 0.07772021 0.1036269
##   Promo 2 0.05181347 0.05181347 0.1295337
## 
## , , Z = Wilayah 2
## 
##          Y
## X           Produk A   Produk B  Produk C
##   Promo 1 0.04145078 0.09326425 0.1243523
##   Promo 2 0.06217617 0.08290155 0.1554404

Interpretasi:

Misal \(P_{111} = \frac{5}{213} \approx 0.0235\) → berarti sekitar 2.35% observasi berada di kombinasi (Promo 1, Produk A, Wilayah 1).

7.2.2 Peluang Marginal (Marginal Probability)

Definisi:

Peluang marginal adalah probabilitas dari kejadian satu variabel kategorik tanpa memperhitungkan nilai dari dua variabel lainnya dalam tabel kontingensi tiga arah.

Tujuan:

  • Memberikan distribusi individual untuk masing-masing variabel.
  • Dasar untuk membandingkan frekuensi antar kategori.
  • Digunakan untuk menghitung peluang bersyarat dan untuk menguji independensi variabel.

Rumus:

  • Peluang marginal \(X\): \[ P_{i++} = \sum_j \sum_k \frac{n_{ijk}}{n_{+++}} \]
  • Peluang marginal \(Y\): \[ P_{+j+} = \sum_i \sum_k \frac{n_{ijk}}{n_{+++}} \]
  • Peluang marginal \(Z\): \[ P_{++k} = \sum_i \sum_j \frac{n_{ijk}}{n_{+++}} \]

Implementasi R

# Peluang marginal
margin_X <- apply(joint_prob, 1, sum)
margin_Y <- apply(joint_prob, 2, sum)
margin_Z <- apply(joint_prob, 3, sum)

margin_X
##   Promo 1   Promo 2 
## 0.4663212 0.5336788
margin_Y
##  Produk A  Produk B  Produk C 
## 0.1813472 0.3056995 0.5129534
margin_Z
## Wilayah 1 Wilayah 2 
## 0.4404145 0.5595855

Interpretasi:

  • Jika \(P_{\text{Promo 1}++} = 0.2535\), maka 25.35% observasi berasal dari Promo 1 tanpa memperhatikan produk dan wilayah.

7.2.3 Peluang Bersyarat (Conditional Probability)

Definisi:

Peluang bersyarat adalah probabilitas dari kombinasi dua variabel kategori, dengan asumsi nilai variabel ketiga sudah diketahui. Hal ini memberikan wawasan mengenai hubungan dua variabel pada strata tertentu.

Tujuan:

  • Mengkaji pengaruh variabel ketiga terhadap hubungan antara dua variabel lainnya.
  • Menguji independensi bersyarat antar variabel.
  • Dasar untuk pengembangan model log-linear dan analisis lanjut.

Rumus (misal, diketahui \(Z = k\)): \[ P_{ij|k} = \frac{n_{ijk}}{n_{++k}} \]

Implementasi R

# Peluang bersyarat: P(X,Y | Z = "Wilayah 1")
cond_prob_w1 <- prop.table(data_array[, , "Wilayah 1"])
cond_prob_w2 <- prop.table(data_array[, , "Wilayah 2"])

cond_prob_w1
##          Y
## X           Produk A  Produk B  Produk C
##   Promo 1 0.05882353 0.1764706 0.2352941
##   Promo 2 0.11764706 0.1176471 0.2941176
cond_prob_w2
##          Y
## X           Produk A  Produk B  Produk C
##   Promo 1 0.07407407 0.1666667 0.2222222
##   Promo 2 0.11111111 0.1481481 0.2777778

Interpretasi:

  • \(P(\text{Promo 1}, \text{Produk A} | \text{Wilayah 1}) = \frac{5}{105} = 0.0476\) menunjukkan proporsi pengamatan untuk kombinasi tersebut di Wilayah 1.

7.3 Tabel Peluang Bersyarat

Definisi:

Tabel peluang bersyarat menyajikan distribusi peluang dari dua variabel kategorik yang dihitung dengan asumsi bahwa variabel ketiga memiliki nilai tertentu (dikondisikan). Dalam konteks tabel kontingensi tiga arah, hal ini penting untuk melihat bagaimana hubungan antara dua variabel berubah pada setiap level dari variabel ketiga.

Tujuan:

  • Untuk menganalisis apakah hubungan antara dua variabel tetap konsisten pada setiap nilai variabel ketiga.
  • Menjadi dasar analisis independensi bersyarat dan model log-linear.
  • Digunakan untuk mendeteksi interaksi antar variabel dan potensi efek perancu (confounding).

Landasan Teori:

Dalam statistika, probabilitas bersyarat didefinisikan sebagai: \[ P(A \mid B) = \frac{P(A \cap B)}{P(B)} \] Yang berarti peluang dari \(A\) terjadi dengan syarat bahwa \(B\) telah terjadi. Untuk tiga variabel kategorik (\(X\), \(Y\), \(Z\)), peluang bersyarat dari \(X\) dan \(Y\) terhadap nilai \(Z=k\) didefinisikan sebagai: \[ P_{ij|k} = \frac{n_{ijk}}{n_{++k}} \]

Contoh Kasus:

Gunakan data dari perhitungan sebelumnya.

# Peluang bersyarat: P(X,Y | Z = Wilayah 1)
cond_prob_w1 <- prop.table(data_array[, , "Wilayah 1"])
cond_prob_w2 <- prop.table(data_array[, , "Wilayah 2"])

knitr::kable(cond_prob_w1, caption = "Tabel Peluang Bersyarat Wilayah 1")
Tabel Peluang Bersyarat Wilayah 1
Produk A Produk B Produk C
Promo 1 0.0588235 0.1764706 0.2352941
Promo 2 0.1176471 0.1176471 0.2941176 
knitr::kable(cond_prob_w2, caption = "Tabel Peluang Bersyarat Wilayah 2")
Tabel Peluang Bersyarat Wilayah 2
Produk A Produk B Produk C
Promo 1 0.0740741 0.1666667 0.2222222
Promo 2 0.1111111 0.1481481 0.2777778

Interpretasi:

  • Di Wilayah 1, produk C lebih banyak dikaitkan dengan Promo 2 (0.2381) dibanding Promo 1 (0.1429).
  • Di Wilayah 2, tren serupa terlihat, tetapi perbedaannya semakin besar: Promo 2 mendominasi semua kategori produk.
  • Tabel ini menunjukkan bahwa hubungan antara Promo dan Produk bervariasi tergantung wilayah → indikasi adanya interaksi atau perbedaan pola distribusi antar strata.

7.4 Ukuran Asosiasi

Definisi dan Landasan Teori:

Ukuran asosiasi merupakan indikator statistik yang digunakan untuk mengetahui kekuatan dan arah hubungan antara dua variabel kategorik. Dalam konteks tabel kontingensi tiga arah, ukuran asosiasi dihitung secara terpisah pada setiap strata dari variabel ketiga. Hal ini memungkinkan untuk mengevaluasi apakah hubungan antara dua variabel tetap konsisten atau berubah tergantung strata yang ditinjau.

Ukuran asosiasi umum yang digunakan meliputi Risk Difference (RD), Risk Ratio (RR), dan Odds Ratio (OR), yang masing-masing memiliki interpretasi dan kegunaan tersendiri, terutama dalam konteks epidemiologi, pemasaran, dan ilmu sosial.

Tujuan:

  • Mengukur kekuatan dan arah hubungan antara dua variabel kategorik dalam setiap strata.
  • Membandingkan hubungan antar kelompok untuk mendeteksi efek perancu atau interaksi.
  • Memberikan dasar kuantitatif dalam analisis hubungan tiga arah sebelum melakukan uji statistik lanjutan atau model log-linear.

Ukuran asosiasi seperti Risk Difference (RD), Risk Ratio (RR), dan Odds Ratio (OR) digunakan untuk mengukur kekuatan dan arah hubungan antara dua variabel kategorik pada setiap nilai dari variabel ketiga (strata). Setiap ukuran memiliki interpretasi dan kegunaan tersendiri tergantung pada konteks data.

7.4.1 Risk Difference (RD)

Definisi:

Risk Difference adalah selisih antara dua probabilitas risiko, yaitu probabilitas terjadinya kejadian pada dua kelompok yang berbeda. Dalam konteks \(2 \times 2 \times K\):

\[ RD_k = p_{1\cdot k} - p_{2\cdot k} \]

Keterangan:

  • \(p_{1\cdot k}\): proporsi kejadian di kelompok 1 pada strata ke-\(k\)
  • \(p_{2\cdot k}\): proporsi kejadian di kelompok 2 pada strata ke-\(k\)

Contoh Perhitungan:

Data yang digunakan ialah data yang mengacu pada perhitungan sebelumnya

rd_w1 <- (5 / (5 + 10)) - (10 / (10 + 20))  # Wilayah 1
rd_w2 <- (8 / (8 + 12)) - (16 / (16 + 24))  # Wilayah 2
rd_w1
## [1] 0
rd_w2
## [1] 0

Interpretasi:

  • \(RD_1 = -0.0556\): risiko ProdukA lebih rendah pada Promo1 dibanding Promo2 di Wilayah 1
  • \(RD_2 = -0.0667\): hasil serupa di Wilayah 2

7.4.2 Risk Ratio (RR)

Definisi:

Risk Ratio adalah perbandingan dua risiko kejadian dari dua kelompok. Dalam bentuk:

\[ RR_k = \frac{p_{1\cdot k}}{p_{2\cdot k}} \]

Contoh Perhitungan:

rr_w1 <- (5 / (5 + 10)) / (10 / (10 + 20))  # Wilayah 1
rr_w2 <- (8 / (8 + 12)) / (16 / (16 + 24))  # Wilayah 2
rr_w1
## [1] 1
rr_w2
## [1] 1

Interpretasi:

  • \(RR_1 = 0.75\): risiko Promo1 terhadap ProdukA lebih rendah 25% dibanding Promo2 di Wilayah 1
  • \(RR_2 = 0.75\): hasil konsisten di Wilayah 2

7.4.3 Odds Ratio (OR)

Definisi:

Odds Ratio adalah rasio antara odds dari kejadian di dua kelompok. Dalam tabel \(2 \times 2 \times K\):

\[ OR_k = \frac{n_{11k} \cdot n_{22k}}{n_{12k} \cdot n_{21k}} \]

Contoh Perhitungan:

or_w1 <- (5 * 20) / (10 * 10)  # Wilayah 1
or_w2 <- (8 * 24) / (12 * 16)  # Wilayah 2
or_w1
## [1] 1
or_w2
## [1] 1

Interpretasi:

  • \(OR_1 = 1\): tidak ada perbedaan odds antara Promo dan Produk di Wilayah 1
  • \(OR_2 = 1\): hasil konsisten di Wilayah 2

Struktur Data Tabel Kontingensi dan Tabel Marginal:

# Buat array data tiga arah
array_data <- array(data = c(
  5, 10, 15,     # Promo 1 - Wilayah 1
  10, 20, 25,    # Promo 2 - Wilayah 1
  8, 12, 18,     # Promo 1 - Wilayah 2
  16, 24, 30     # Promo 2 - Wilayah 2
), dim = c(2, 3, 2), 
   dimnames = list(
     X = c("Promo 1", "Promo 2"),
     Y = c("Produk A", "Produk B", "Produk C"),
     Z = c("Wilayah 1", "Wilayah 2")
))

# Tampilkan array tiga dimensi
array_data
## , , Z = Wilayah 1
## 
##          Y
## X         Produk A Produk B Produk C
##   Promo 1        5       15       20
##   Promo 2       10       10       25
## 
## , , Z = Wilayah 2
## 
##          Y
## X         Produk A Produk B Produk C
##   Promo 1        8       18       24
##   Promo 2       12       16       30
library(kableExtra)
# Tabel Marginal: Menjumlahkan seluruh Wilayah (Z)
marginal <- apply(array_data, c(1, 2), sum)
df_marginal <- as.data.frame(marginal)
colnames(df_marginal) <- c("Produk A", "Produk B", "Produk C")
rownames(df_marginal) <- c("Promo 1", "Promo 2")

kable(df_marginal, format = "latex", booktabs = TRUE, align = "c",
      caption = "Tabel Marginal (Total untuk Semua Wilayah)") %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  column_spec(1:4, latex_valign = "m", width = "3cm")

Tabel marginal \(X imes Y\) memberikan total frekuensi dari setiap kombinasi kategori Promo dan Produk setelah seluruh Wilayah dijumlahkan.

Ukuran asosiasi seperti Risk Difference (RD), Risk Ratio (RR), dan Odds Ratio (OR) digunakan untuk mengukur kekuatan dan arah hubungan antara dua variabel kategorik pada setiap nilai dari variabel ketiga (strata). Setiap ukuran memiliki interpretasi dan kegunaan tersendiri tergantung pada konteks data.

Kesimpulan Umum:

Nilai RD negatif, RR < 1, dan OR = 1 di kedua wilayah menunjukkan tidak ada asosiasi positif antara Promo1 dan ProdukA. Promo2 justru cenderung lebih efektif, namun perbedaannya relatif konsisten antar wilayah → tidak ada interaksi besar.

7.5 Tabel Kontingensi Parsial dan Ukuran Asosiasi (RR, RD, OR)

Tujuan:

Subbab ini bertujuan menyajikan tabel kontingensi parsial untuk masing-masing level variabel ketiga (Z = Wilayah), serta menghitung dan menginterpretasi ukuran asosiasi yaitu Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR) baik secara manual maupun menggunakan R.

Contoh Data:

Gunakan data dari array array_data sebelumnya, lalu fokus pada dua kategori produk pertama: Produk A dan Produk B. Struktur tabel 2x2 akan:

  • Baris: Promo 1 dan Promo 2
  • Kolom: Produk A dan Produk B

7.5.1 Wilayah 1: Tabel Parsial

tab_w1 <- array_data[, 1:2, "Wilayah 1"]
kable(tab_w1, format = "latex", booktabs = TRUE, align = "c", caption = "Tabel Kontingensi Parsial Wilayah 1") %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  column_spec(1:3, latex_valign = "m", width = "3cm")

Perhitungan Manual:

Diketahui: - \(n_{11} = 5\), \(n_{12} = 10\) (Promo 1) - \(n_{21} = 10\), \(n_{22} = 20\) (Promo 2)

  • \(RD = \frac{5}{15} - \frac{10}{30} = 0.333 - 0.333 = 0\)
  • \(RR = \frac{5/15}{10/30} = 1\)
  • \(OR = \frac{5 \times 20}{10 \times 10} = 100/100 = 1\)
# Wilayah 1 - RR, RD, OR
n11 <- 5; n12 <- 10; n21 <- 10; n22 <- 20
n1. <- n11 + n12; n2. <- n21 + n22

rd1 <- (n11 / n1.) - (n21 / n2.)
rr1 <- (n11 / n1.) / (n21 / n2.)
or1 <- (n11 * n22) / (n12 * n21)

list(RD = rd1, RR = rr1, OR = or1)
## $RD
## [1] 0
## 
## $RR
## [1] 1
## 
## $OR
## [1] 1

7.5.2 Wilayah 2: Tabel Parsial

tab_w2 <- array_data[, 1:2, "Wilayah 2"]
kable(tab_w2, format = "latex", booktabs = TRUE, align = "c", caption = "Tabel Kontingensi Parsial Wilayah 2") %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  column_spec(1:3, latex_valign = "m", width = "3cm")

Perhitungan Manual:

  • \(n_{11} = 8\), \(n_{12} = 12\) (Promo 1)

  • \(n_{21} = 16\), \(n_{22} = 24\) (Promo 2)

  • \(RD = \frac{8}{20} - \frac{16}{40} = 0.4 - 0.4 = 0\)

  • \(RR = 1\)

  • \(OR = \frac{8 \times 24}{12 \times 16} = 192/192 = 1\)

# Wilayah 2 - RR, RD, OR
n11 <- 8; n12 <- 12; n21 <- 16; n22 <- 24
n1. <- n11 + n12; n2. <- n21 + n22

rd2 <- (n11 / n1.) - (n21 / n2.)
rr2 <- (n11 / n1.) / (n21 / n2.)
or2 <- (n11 * n22) / (n12 * n21)

list(RD = rd2, RR = rr2, OR = or2)
## $RD
## [1] 0
## 
## $RR
## [1] 1
## 
## $OR
## [1] 1

Kesimpulan Interpretatif:

  • RD dan RR = 0 & 1 di kedua wilayah, mengindikasikan tidak ada perbedaan risiko maupun asosiasi antara jenis promo dan preferensi Produk A/B.
  • Nilai OR = 1 juga memperkuat bahwa tidak ada kecondongan odds pada promo tertentu.
  • Analisis ini memperlihatkan bahwa pada masing-masing strata (Wilayah 1 & 2), hubungan antara promo dan preferensi produk stabil dan tidak signifikan secara statistik.

Jika nanti ditemukan perbedaan nilai OR antar wilayah, maka itu bisa ditindaklanjuti ke subbab independensi bersyarat.

7.6 Conditional Independence

Definisi

Conditional independence (kemandirian bersyarat) dalam tabel kontingensi terjadi ketika dua variabel menjadi independen setelah dikendalikan oleh variabel ketiga.

Secara matematis, dua variabel X dan Y dikatakan independen secara kondisional terhadap variabel Z jika:

\[ P(X, Y | Z) = P(X | Z) \cdot P(Y | Z) \]

atau dalam bentuk frekuensi:

\[ \frac{n_{ijk}}{n_{++k}} = \frac{n_{i+k}}{n_{++k}} \cdot \frac{n_{+jk}}{n_{++k}} \]

Contoh Kasus Baru

Misalkan kita memiliki data tentang Gaya Hidup (X), Kondisi Tidur (Y), dan Kelompok Umur (Z).

  • X (Gaya Hidup): Sehat, Tidak Sehat
  • Y (Tidur): Nyenyak, Tidak Nyenyak
  • Z (Umur): Remaja, Dewasa
library(knitr)
library(kableExtra)

# Data array 3 dimensi: Gaya Hidup, Tidur, Umur
lifestyle_sleep <- array(
  c(25, 15, 30, 10, 20, 30, 25, 35),
  dim = c(2, 2, 2),
  dimnames = list(
    "Gaya Hidup" = c("Sehat", "Tidak Sehat"),
    "Tidur" = c("Nyenyak", "Tidak Nyenyak"),
    "Umur" = c("Remaja", "Dewasa")
  )
)

# Tampilkan tabel 2x2 untuk masing-masing Umur
for (umur in dimnames(lifestyle_sleep)$Umur) {
  df <- as.data.frame.matrix(lifestyle_sleep[,,umur])
  colnames(df) <- c("Nyenyak", "Tidak Nyenyak")
  rownames(df) <- c("Sehat", "Tidak Sehat")
  print(
    kable(df, format = "latex", booktabs = TRUE,
          caption = paste("Tabel Gaya Hidup vs Tidur untuk Umur =", umur)) %>%
      kable_styling(full_width = FALSE, position = "center")
  )
}
## \begin{table}
## \centering
## \caption{\label{tab:unnamed-chunk-48}Tabel Gaya Hidup vs Tidur untuk Umur = Remaja}
## \centering
## \begin{tabular}[t]{lrr}
## \toprule
##   & Nyenyak & Tidak Nyenyak\\
## \midrule
## Sehat & 25 & 30\\
## Tidak Sehat & 15 & 10\\
## \bottomrule
## \end{tabular}
## \end{table}
## \begin{table}
## \centering
## \caption{\label{tab:unnamed-chunk-48}Tabel Gaya Hidup vs Tidur untuk Umur = Dewasa}
## \centering
## \begin{tabular}[t]{lrr}
## \toprule
##   & Nyenyak & Tidak Nyenyak\\
## \midrule
## Sehat & 20 & 25\\
## Tidak Sehat & 30 & 35\\
## \bottomrule
## \end{tabular}
## \end{table}

Perhitungan Manual Probabilitas Bersyarat

Contoh: Untuk kelompok umur Remaja:

  • \(P(Tidur Nyenyak | Sehat, Remaja) = 25 / (25+15) = 0.625\)
  • \(P(Tidur Nyenyak | Tidak Sehat, Remaja) = 30 / (30+10) = 0.75\)
  • \(P(Tidur Nyenyak | Remaja) = (25+30) / (25+15+30+10) = 55 / 80 = 0.6875\)

Jika nilai probabilitas bersyarat tidak berbeda jauh dari probabilitas marjinal, maka X dan Y bersifat conditional independent terhadap Z.

Perhitungan dengan R

# Uji chi-square untuk masing-masing strata umur
chisq_remaja <- chisq.test(lifestyle_sleep[, , "Remaja"])
chisq_dewasa <- chisq.test(lifestyle_sleep[, , "Dewasa"])

chisq_remaja
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  lifestyle_sleep[, , "Remaja"]
## X-squared = 0.93091, df = 1, p-value = 0.3346
chisq_dewasa
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  lifestyle_sleep[, , "Dewasa"]
## X-squared = 0, df = 1, p-value = 1

Interpretasi Hasil

  • Jika p-value > 0.05 → tidak ada hubungan signifikan antara Gaya Hidup dan Tidur setelah dikendalikan oleh Umur (conditional independence).
  • Jika p-value < 0.05 → masih ada ketergantungan antara Gaya Hidup dan Tidur meskipun Umur dikendalikan.

Tabel Marginal

# Marginal XY
margin_XY <- apply(lifestyle_sleep, c(1, 2), sum)
cat("\n## Tabel Marginal: Gaya Hidup x Tidur\n")
## 
## ## Tabel Marginal: Gaya Hidup x Tidur
print(kable(margin_XY, format = "latex", booktabs = TRUE) %>%
        kable_styling(full_width = FALSE, position = "center"))
## \begin{table}
## \centering
## \begin{tabular}{lrr}
## \toprule
##   & Nyenyak & Tidak Nyenyak\\
## \midrule
## Sehat & 45 & 55\\
## Tidak Sehat & 45 & 45\\
## \bottomrule
## \end{tabular}
## \end{table}
# Marginal XZ
margin_XZ <- apply(lifestyle_sleep, c(1, 3), sum)
cat("\n## Tabel Marginal: Gaya Hidup x Umur\n")
## 
## ## Tabel Marginal: Gaya Hidup x Umur
print(kable(margin_XZ, format = "latex", booktabs = TRUE) %>%
        kable_styling(full_width = FALSE, position = "center"))
## \begin{table}
## \centering
## \begin{tabular}{lrr}
## \toprule
##   & Remaja & Dewasa\\
## \midrule
## Sehat & 55 & 45\\
## Tidak Sehat & 25 & 65\\
## \bottomrule
## \end{tabular}
## \end{table}
# Marginal YZ
margin_YZ <- apply(lifestyle_sleep, c(2, 3), sum)
cat("\n## Tabel Marginal: Tidur x Umur\n")
## 
## ## Tabel Marginal: Tidur x Umur
print(kable(margin_YZ, format = "latex", booktabs = TRUE) %>%
        kable_styling(full_width = FALSE, position = "center"))
## \begin{table}
## \centering
## \begin{tabular}{lrr}
## \toprule
##   & Remaja & Dewasa\\
## \midrule
## Nyenyak & 40 & 50\\
## Tidak Nyenyak & 40 & 60\\
## \bottomrule
## \end{tabular}
## \end{table}

Analisis lebih lanjut dapat dilakukan menggunakan odds ratio atau model log-linear untuk melihat konsistensi hubungan antara variabel setelah mengendalikan variabel ketiga.

7.7 Ukuran Asosiasi Berdasarkan Tabel Marginal Y dan X

Teori Dasar

Ukuran asosiasi dalam tabel kontingensi digunakan untuk mengevaluasi kekuatan dan arah hubungan antara dua variabel kategorik. Dalam konteks tabel tiga arah, ukuran asosiasi seperti Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR) dihitung pada setiap strata dari variabel ketiga (Z) untuk mengevaluasi adanya pengaruh atau efek pembaur (confounding) dari variabel tersebut.

Tujuan:

  • Mengukur seberapa besar pengaruh dari satu variabel terhadap variabel lain.
  • Membandingkan kekuatan asosiasi antar strata.
  • Mendeteksi kemungkinan efek pembaur atau interaksi.

Karakteristik:

  • RD (Risk Difference): Mengukur selisih risiko antara dua kelompok.
  • RR (Risk Ratio): Mengukur rasio antara dua risiko.
  • OR (Odds Ratio): Mengukur perbandingan odds (peluang) antara dua kelompok.

Ukuran ini sangat penting dalam epidemiologi, ilmu sosial, dan penelitian kesehatan untuk mengevaluasi efek perlakuan atau faktor risiko terhadap suatu outcome.

7.8 Marginal Y dan X

Definisi

Analisis marginal antara variabel Y (Tidur) dan X (Gaya Hidup), tanpa mempertimbangkan variabel Z (Umur).

# Marginal antara Tidur dan Gaya Hidup
margin_YX <- apply(lifestyle_sleep, c(2, 1), sum)

# Tampilkan tabel marginal Y dan X
kable(margin_YX, format = "latex", booktabs = TRUE,
      caption = "Tabel Marginal antara Tidur dan Gaya Hidup") %>%
  kable_styling(full_width = FALSE, position = "center")

Perhitungan Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR)

Kita menghitung RD, RR, dan OR untuk setiap tingkat Z (Remaja dan Dewasa) berdasarkan data lifestyle_sleep.

Remaja

# Nilai frekuensi
ar <- lifestyle_sleep["Sehat", "Nyenyak", "Remaja"]
br <- lifestyle_sleep["Sehat", "Tidak Nyenyak", "Remaja"]
cr <- lifestyle_sleep["Tidak Sehat", "Nyenyak", "Remaja"]
dr <- lifestyle_sleep["Tidak Sehat", "Tidak Nyenyak", "Remaja"]

# Perhitungan
p1r <- ar / (ar + br)
p2r <- cr / (cr + dr)
RD_remaja <- p1r - p2r
RR_remaja <- p1r / p2r
OR_remaja <- (ar * dr) / (br * cr)

# Tabel hasil
data_RD_RR_OR_remaja <- data.frame(RD = RD_remaja, RR = RR_remaja, OR = OR_remaja)
kable(data_RD_RR_OR_remaja, format = "latex", booktabs = TRUE,
      caption = "Risk Difference, Relative Risk, dan Odds Ratio untuk Remaja") %>%
  kable_styling(full_width = FALSE, position = "center")

Dewasa

# Nilai frekuensi
ad <- lifestyle_sleep["Sehat", "Nyenyak", "Dewasa"]
bd <- lifestyle_sleep["Sehat", "Tidak Nyenyak", "Dewasa"]
cd <- lifestyle_sleep["Tidak Sehat", "Nyenyak", "Dewasa"]
dd <- lifestyle_sleep["Tidak Sehat", "Tidak Nyenyak", "Dewasa"]

# Perhitungan
p1d <- ad / (ad + bd)
p2d <- cd / (cd + dd)
RD_dewasa <- p1d - p2d
RR_dewasa <- p1d / p2d
OR_dewasa <- (ad * dd) / (bd * cd)

# Tabel hasil
data_RD_RR_OR_dewasa <- data.frame(RD = RD_dewasa, RR = RR_dewasa, OR = OR_dewasa)
kable(data_RD_RR_OR_dewasa, format = "latex", booktabs = TRUE,
      caption = "Risk Difference, Relative Risk, dan Odds Ratio untuk Dewasa") %>%
  kable_styling(full_width = FALSE, position = "center")

Kesimpulan

Analisis tabel kontingensi tiga arah ini menunjukkan bahwa Umur (Z) mempengaruhi hubungan antara Gaya Hidup (X) dan Kondisi Tidur (Y). Efek positif dari gaya hidup sehat terhadap kualitas tidur lebih besar pada kelompok usia remaja dibandingkan dewasa. Hasil ini dapat digunakan untuk menyusun rekomendasi kesehatan yang lebih spesifik berdasarkan kelompok umur.

7.9 Inferensi Tabel Kontingensi Tiga Arah

Tabel kontingensi tiga arah merupakan alat analisis yang digunakan untuk mengevaluasi hubungan antara dua variabel kategorik (\(X\) dan \(Y\)) dengan mempertimbangkan pengaruh variabel ketiga sebagai variabel kontrol (\(Z\)). Dalam praktiknya, pendekatan ini membantu mendeteksi apakah hubungan antara \(X\) dan \(Y\) tetap konsisten dalam setiap level \(Z\), atau apakah terdapat pengaruh pembaur (confounding) atau interaksi.

Tabel Parsial dan Marginal

Tabel kontingensi tiga arah tersusun atas beberapa tabel parsial \(2 \times 2\) berdasarkan setiap strata \(Z\) (misalnya: usia muda, dewasa, tua), dan satu tabel marginal yang mengabaikan nilai \(Z\).

Sebagai contoh, untuk mempelajari hubungan antara kebiasaan merokok (\(X\)) dan kanker paru-paru (\(Y\)) dengan kontrol usia (\(Z\)), kita bisa membuat dua tabel parsial berdasarkan strata Z: - Usia Muda - Usia Tua

Setiap tabel parsial memiliki struktur berikut:

Ukuran asosiasi yang umum digunakan di sini adalah Odds Ratio (OR):

\[ OR = \frac{a \cdot d}{b \cdot c} \]

Jika nilai OR dari seluruh tabel parsial konsisten (tidak berbeda jauh), maka bisa dilakukan penggabungan dengan menghitung Odds Ratio Mantel-Haenszel sebagai estimasi keseluruhan efek asosiasi.

7.9.1 Independensi Bersyarat

Definisi: Independensi bersyarat terjadi apabila dua variabel (\(X\) dan \(Y\)) tidak memiliki hubungan secara statistik dalam setiap strata dari variabel ketiga (\(Z\)). Artinya, setelah mengendalikan pengaruh \(Z\), tidak ada ketergantungan antara \(X\) dan \(Y\).

Secara matematis ditulis sebagai:

\[ OR(X, Y \mid Z) = 1 \]

Artinya, Odds Ratio di setiap level \(Z\) sama dengan 1, sehingga tidak ada asosiasi antara \(X\) dan \(Y\) setelah mengontrol \(Z\).

Catatan Penting: - Jika \(X\) dan \(Y\) independen bersyarat pada \(Z\), belum tentu mereka independen secara marginal. Hal ini dikenal sebagai Paradoks Simpson, di mana arah hubungan pada strata dapat berlawanan dengan arah hubungan pada tabel gabungan. - Independensi bersyarat sangat penting dalam analisis epidemiologi, sosial, serta machine learning (misalnya dalam struktur jaringan Bayesian).

7.9.2 Uji Cochran-Mantel-Haenszel (CMH)

Untuk menguji independensi bersyarat secara statistik, digunakan Uji Cochran-Mantel-Haenszel (CMH). Uji ini menghitung apakah terdapat asosiasi antara \(X\) dan \(Y\) dengan mengontrol variabel \(Z\). CMH sangat berguna ketika data terbagi dalam beberapa strata dan kita ingin menyimpulkan asosiasi keseluruhan tanpa efek pembaur dari \(Z\).

Uji CMH menguji apakah terdapat asosiasi antara dua variabel kategorik (\(X\) dan \(Y\)), dengan mengendalikan pengaruh variabel ke-3 (\(Z\)). Rumus statistik CMH adalah:

\[ CMH = \frac{\left[ \sum_k (n_{11k} - \mu_{11k}) \right]^2}{\sum_k \text{Var}(n_{11k})} \]

dengan:

  • \(n_{11k}\): Observasi pada sel baris 1 kolom 1 dari strata ke-\(k\)
  • \(\mu_{11k} = \frac{n_{1+k} n_{+1k}}{n_{++k}}\): Ekspektasi dari \(n_{11k}\)
  • \(\text{Var}(n_{11k}) = \frac{n_{1+k} n_{2+k} n_{+1k} n_{+2k}}{n_{++k}^2 (n_{++k}-1)}\): Variansi

Hipotesis:

  • \(H_0\): Tidak ada asosiasi antara \(X\) dan \(Y\) di setiap strata \(Z\) → \(\theta_{xy(k)} = 1\)
  • \(H_1\): Terdapat asosiasi di paling tidak satu strata → \(\theta_{xy(k)} \neq 1\)

Statistik uji Cochran-Mantel-Haenszel (CMH) digunakan untuk menentukan apakah terdapat hubungan antara dua variabel kategorik (\(X\) dan \(Y\)) setelah mengendalikan pengaruh variabel ketiga (\(Z\)).

Rumus umum CMH dinyatakan sebagai:

\[ CMH = \frac{\left[ \sum_k (n_{11k} - \mu_{11k}) \right]^2}{\sum_k \text{Var}(n_{11k})} \]

Keterangan:

  • \(n_{11k}\): Jumlah observasi pada sel baris 1 kolom 1 dari strata ke-\(k\)
  • \(\mu_{11k}\): Ekspektasi dari \(n_{11k}\), dihitung sebagai:

\[ \mu_{11k} = \frac{n_{1+k} \cdot n_{+1k}}{n_{++k}} \]

  • \(\text{Var}(n_{11k})\): Variansi dari \(n_{11k}\), diberikan oleh:

\[ \text{Var}(n_{11k}) = \frac{n_{1+k} \cdot n_{2+k} \cdot n_{+1k} \cdot n_{+2k}}{n_{++k}^2 (n_{++k} - 1)} \]

Distribusi:

Statistik CMH mengikuti distribusi chi-square dengan derajat kebebasan = 1.

Keputusan Uji:

  • Tolak \(H_0\) jika nilai CMH > \(\chi^2_{(1, \alpha)}\) atau jika p-value < \(\alpha\).

Uji ini sangat bermanfaat dalam studi epidemiologi dan eksperimen sosial yang membutuhkan pengendalian terhadap faktor pembaur (confounder), terutama saat menganalisis hubungan dua variabel dalam konteks data berlapis.

7.9.2.1 Contoh Kasus 1 CMH: Pola Konsumsi, Kualitas Tidur, dan Lingkungan

Pada contoh ini, kita ingin menganalisis apakah terdapat hubungan antara Pola Konsumsi dan Kualitas Tidur, dengan mempertimbangkan Lingkungan Tempat Tinggal sebagai variabel kontrol. Ini adalah kasus tipikal dalam studi gaya hidup dan kesehatan masyarakat.

Variabel yang digunakan:

  • \(X\) = Pola Konsumsi (Seimbang, Tidak Seimbang)
  • \(Y\) = Kualitas Tidur (Baik, Buruk)
  • \(Z\) = Lingkungan (Perkotaan, Pedesaan)
# Data array 2x2x2: Pola Konsumsi x Kualitas Tidur x Lingkungan
konsumsi_tidur <- array(
  c(40, 20, 30, 30, 25, 35, 20, 40),
  dim = c(2, 2, 2),
  dimnames = list(
    Konsumsi = c("Seimbang", "Tidak Seimbang"),
    Tidur = c("Baik", "Buruk"),
    Lingkungan = c("Perkotaan", "Pedesaan")
  )
)

Tabel Kontingensi Per Strata Lingkungan

Baik Buruk
Seimbang 40 20
Tidak Seimbang 25 35

Tabel Pola Konsumsi vs Tidur - Pedesaan

Baik Buruk
Seimbang 30 30
Tidak Seimbang 20 40

Perhitungan Manual Probabilitas Bersyarat

Langkah 1: Menghitung nilai ekspektasi \(\mu_{11k}\)

Rumus: \[ \mu_{11k} = \frac{n_{1+k} \cdot n_{+1k}}{n_{++k}} \]

  • Perkotaan:
    • \(n_{1+1} = 40 + 20 = 60\) (Total Seimbang)
    • \(n_{+1,1} = 40 + 25 = 65\) (Total Tidur Baik)
    • \(n_{++1} = 40 + 20 + 25 + 35 = 120\)
    • \(\mu_{111} = \frac{60 \cdot 65}{120} = 32.5\)
  • Pedesaan:
    • \(n_{1+2} = 30 + 30 = 60\)
    • \(n_{+1,2} = 30 + 20 = 50\)
    • \(n_{++2} = 30 + 30 + 20 + 40 = 120\)
    • \(\mu_{112} = \frac{60 \cdot 50}{120} = 25\)

Langkah 2: Menghitung varians \(\text{Var}(n_{11k})\)

Rumus: \[ \text{Var}(n_{11k}) = \frac{n_{1+k} n_{2+k} n_{+1k} n_{+2k}}{n_{++k}^2 (n_{++k} - 1)} \]

  • Perkotaan:
    • \(n_{1+1} = 60\), \(n_{2+1} = 60\)
    • \(n_{+1,1} = 65\), \(n_{+2,1} = 55\)
    • \(\text{Var}(n_{111}) = \frac{60 \cdot 60 \cdot 65 \cdot 55}{120^2 \cdot 119} \approx 8.97\)
  • Pedesaan:
    • \(n_{1+2} = 60\), \(n_{2+2} = 60\)
    • \(n_{+1,2} = 50\), \(n_{+2,2} = 70\)
    • \(\text{Var}(n_{112}) = \frac{60 \cdot 60 \cdot 50 \cdot 70}{120^2 \cdot 119} \approx 10.26\)

Langkah 3: Menghitung Statistik CMH

\[ X^2_{CMH} = \frac{(n_{111} - \mu_{111} + n_{112} - \mu_{112})^2}{\text{Var}(n_{111}) + \text{Var}(n_{112})} \]

  • \(n_{111} = 40\), \(\mu_{111} = 32.5\)
  • \(n_{112} = 30\), \(\mu_{112} = 25\)

\[ X^2_{CMH} = \frac{(40 - 32.5 + 30 - 25)^2}{8.97 + 10.26} = \frac{(12.5)^2}{19.23} = \frac{156.25}{19.23} \approx 8.13 \]

Interpretasi: Dengan \(X^2_{CMH} \approx 8.13\) dan derajat kebebasan 1, kita bandingkan dengan nilai kritis \(\chi^2_{0.05,1} = 3.841\). Karena 8.13 > 3.841, maka kita menolak hipotesis nol dan menyimpulkan bahwa pola konsumsi berpengaruh signifikan terhadap kualitas tidur setelah mengontrol lingkungan tempat tinggal.

Uji CMH dengan Base R

# Uji Cochran-Mantel-Haenszel tanpa continuity correction
cmh_konsumsi <- mantelhaen.test(konsumsi_tidur, correct = FALSE)
cmh_konsumsi
## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  konsumsi_tidur
## Mantel-Haenszel X-squared = 3.8945, df = 1, p-value = 0.04844
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.003627 2.853555
## sample estimates:
## common odds ratio 
##          1.692308

Interpretasi Hasil Uji CMH

  • Jika p-value < 0.05, maka terdapat hubungan signifikan antara Pola Konsumsi dan Tidur setelah mengendalikan Lingkungan → tolak \(H_0\)
  • Jika p-value > 0.05, maka tidak cukup bukti adanya asosiasi setelah mengontrol Z → dukung adanya independensi bersyarat

Kesimpulan

Analisis Cochran-Mantel-Haenszel memungkinkan kita mengevaluasi hubungan dua variabel kategorik dengan mempertimbangkan faktor pembaur. Dalam studi ini, uji CMH membantu kita memahami apakah pola konsumsi benar-benar memengaruhi kualitas tidur, setelah mengontrol kondisi lingkungan tempat tinggal.

7.9.2.2 Contoh Kasus 2 CMH: Pelatihan, Peningkatan Keterampilan, dan Industri

Dalam contoh ini, kita ingin mengetahui apakah partisipasi dalam program pelatihan berpengaruh terhadap peningkatan keterampilan, dengan mengontrol faktor industri tempat bekerja.

Spesifikasi Variabel:

Simbol Variabel Kategori Peran dalam Analisis
\(X\) Pelatihan Mengikuti, Tidak Mengikuti Prediktor / independen
\(Y\) Peningkatan Keterampilan Ya, Tidak Respons / dependen
\(Z\) Industri Stratum 1, Stratum 2, Stratum 3 Kontrol / confounder 
# Data array: Pelatihan x Peningkatan x Industri
pelatihan_data <- array(c(
  20, 30, 10, 40,   # Stratum 1
  15, 25, 5, 35,    # Stratum 2
  18, 32, 12, 38    # Stratum 3
),
  dim = c(2, 2, 3),
  dimnames = list(
    Pelatihan = c("Mengikuti", "Tidak Mengikuti"),
    Peningkatan = c("Ya", "Tidak"),
    Industri = c("Stratum 1", "Stratum 2", "Stratum 3")
  )
)
# Uji CMH tanpa continuity correction
cmh_pelatihan <- mantelhaen.test(pelatihan_data, correct = FALSE)
cmh_pelatihan
## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  pelatihan_data
## Mantel-Haenszel X-squared = 11.733, df = 1, p-value = 0.0006139
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.482202 4.377281
## sample estimates:
## common odds ratio 
##          2.547159

Interpretasi Hasil: - Statistik uji Mantel-Haenszel menghasilkan nilai \(\chi^2 \approx 11.73\) dengan \(p\)-value \(\approx 0.0006\) - Karena \(p < 0.05\), maka kita menolak hipotesis nol. - Artinya, terdapat hubungan signifikan antara pelatihan dan peningkatan keterampilan setelah mempertimbangkan pengaruh industri.

Kesimpulan: Partisipasi dalam program pelatihan terbukti secara statistik meningkatkan peluang peningkatan keterampilan pekerja, bahkan setelah mengontrol perbedaan antar industri.

7.9.2.3 Contoh Kasus 3 CMH: Pengaruh Intervensi Digital terhadap Kepatuhan Pasien dengan Provinsi sebagai Kontrol

Diketahui sebuah studi dilakukan untuk menilai apakah intervensi digital (misalnya aplikasi pengingat obat) berpengaruh terhadap kepatuhan pasien (\(Y\)), dengan mempertimbangkan pengaruh provinsi tempat tinggal sebagai variabel kontrol (\(Z\)). Analisis dilakukan melalui odds ratio per provinsi, dilanjutkan dengan uji CMH.

Langkah 1: Menyusun Data

library(knitr)
# Data dalam format data.frame
data <- data.frame(
  Provinsi = rep(c("Jakarta", "Bandung", "Surabaya", "Medan", "Yogyakarta", "Makassar", "Denpasar", "Palembang"), each = 2),
  Intervensi = rep(c("Ya", "Tidak"), 8),
  Patuh = c(130, 40, 110, 50, 100, 35, 95, 30, 85, 25, 80, 20, 70, 18, 65, 15),
  Tidak_Patuh = c(90, 70, 85, 95, 75, 70, 65, 65, 60, 60, 55, 50, 45, 42, 40, 38)
)
kable(data, caption = "Data Kepatuhan dan Intervensi Digital di 8 Provinsi")
Data Kepatuhan dan Intervensi Digital di 8 Provinsi
Provinsi Intervensi Patuh Tidak_Patuh
Jakarta Ya 130 90
Jakarta Tidak 40 70
Bandung Ya 110 85
Bandung Tidak 50 95
Surabaya Ya 100 75
Surabaya Tidak 35 70
Medan Ya 95 65
Medan Tidak 30 65
Yogyakarta Ya 85 60
Yogyakarta Tidak 25 60
Makassar Ya 80 55
Makassar Tidak 20 50
Denpasar Ya 70 45
Denpasar Tidak 18 42
Palembang Ya 65 40
Palembang Tidak 15 38

Langkah 2: Menyusun Array untuk Uji CMH

# Susun data menjadi array 2x2x8
array_data <- array(
  c(data$Patuh, data$Tidak_Patuh),
  dim = c(2, 2, 8),
  dimnames = list(
    Intervensi = c("Ya", "Tidak"),
    Kepatuhan = c("Patuh", "Tidak Patuh"),
    Provinsi = unique(data$Provinsi)
  )
)

Langkah 3: Perhitungan Odds Ratio per Provinsi

library(epitools)
or_results <- matrix(nrow = 8, ncol = 2)
colnames(or_results) <- c("Provinsi", "Odds Ratio")
for (i in 1:8) {
  subset_data <- data[(2*i-1):(2*i), ]
  table_matrix <- matrix(c(subset_data$Patuh, subset_data$Tidak_Patuh), nrow = 2, byrow = TRUE)
  or_value <- oddsratio.wald(table_matrix)$measure[2, 1]
  or_results[i, ] <- c(subset_data$Provinsi[1], round(or_value, 3))
}
or_results_df <- as.data.frame(or_results)
kable(or_results_df, caption = "Odds Ratio per Provinsi")
Odds Ratio per Provinsi
Provinsi Odds Ratio
Jakarta 2.528
Bandung 2.459
Surabaya 2.667
Medan 3.167
Yogyakarta 3.4
Makassar 3.636
Denpasar 3.63
Palembang 4.117

Langkah 4: Uji CMH untuk Independensi Bersyarat

cmh_test <- mantelhaen.test(array_data)
cmh_test
## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  array_data
## Mantel-Haenszel X-squared = 0.96306, df = 1, p-value = 0.3264
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.9139751 1.3357175
## sample estimates:
## common odds ratio 
##          1.104904

Interpretasi Hasil:

  • Jika nilai p-value < 0.05, maka kita menolak hipotesis nol dan menyimpulkan terdapat hubungan signifikan antara intervensi digital dan kepatuhan pasien, setelah mengontrol efek provinsi.
  • Jika p-value > 0.05, maka tidak terdapat cukup bukti untuk menyatakan adanya hubungan.

Kesimpulan:

Uji CMH merupakan alat penting untuk menguji hubungan dua variabel kategorik dengan kontrol stratifikasi. Dalam contoh ini, hasil OR dan CMH membantu menyimpulkan apakah intervensi teknologi berdampak nyata dalam berbagai wilayah administratif.

7.9.3 Odds Ratio Bersama

Pengantar dan Tujuan

Dalam tabel kontingensi 2 × 2 × K, odds ratio bersama (common odds ratio) digunakan untuk mengukur kekuatan asosiasi antara dua variabel utama (\(X\) dan \(Y\)), setelah mengontrol pengaruh dari variabel ketiga (\(Z\)). Odds ratio bersama sangat berguna ketika odds ratio pada masing-masing strata konsisten atau searah.

Rumus Odds Ratio Bersama (Estimator Mantel-Haenszel)

\[ \hat{\theta}_{MH} = \frac{\sum_k \frac{n_{11k} n_{22k}}{n_{..k}}}{\sum_k \frac{n_{12k} n_{21k}}{n_{..k}}} \]

Keterangan:

  • \(n_{11k}\): Frekuensi sel baris 1 kolom 1 pada strata ke-\(k\)
  • \(n_{12k}\): Frekuensi sel baris 1 kolom 2 pada strata ke-\(k\)
  • \(n_{21k}\): Frekuensi sel baris 2 kolom 1 pada strata ke-\(k\)
  • \(n_{22k}\): Frekuensi sel baris 2 kolom 2 pada strata ke-\(k\)
  • \(n_{..k}\): Total observasi pada strata ke-\(k\)

Contoh Perhitungan dengan data dari Contoh Kasus 3 CMH:

intervensi_data <- array(
  c(130, 90, 70, 40,    # Jakarta
    110, 85, 95, 50,    # Bandung
    100, 75, 70, 35,    # Surabaya
    95, 65, 65, 30,     # Medan
    85, 60, 60, 25,     # Yogyakarta
    80, 55, 50, 20,     # Makassar
    70, 45, 42, 18,     # Denpasar
    65, 40, 38, 15),    # Palembang
  dim = c(2, 2, 8),
  dimnames = list(
    Intervensi = c("Ya", "Tidak"),
    Kepatuhan = c("Patuh", "Tidak Patuh"),
    Provinsi = c("Jakarta", "Bandung", "Surabaya", "Medan", 
                 "Yogyakarta", "Makassar", "Denpasar", "Palembang")
  )
)

Contoh Perhitungan Manual

\[ \hat{\theta}_{MH} = \frac{(130 \cdot 40 + 110 \cdot 50 + 100 \cdot 35 + 95 \cdot 30 + 85 \cdot 25 + 80 \cdot 20 + 70 \cdot 18 + 65 \cdot 15)/n}{(90 \cdot 70 + 85 \cdot 95 + 75 \cdot 70 + 65 \cdot 65 + 60 \cdot 60 + 55 \cdot 50 + 45 \cdot 42 + 40 \cdot 38)/n} \approx \frac{19900/100}{12400/100} = \frac{199}{124} \approx 1.60 \]

7.9.3.1 Standard Error Log Odds Ratio

\[ SE(\log(\hat{\theta}_{MH})) \approx 0.19 \]

7.9.3.2 Interval Kepercayaan 95%

\[ \log(1.60) \pm 1.96 \cdot 0.19 \Rightarrow 0.470 \pm 0.372 \Rightarrow (0.098, 0.842) \] \[ CI\ 95\%\ OR = (e^{0.098}, e^{0.842}) = (1.10, 2.32) \]

7.9.3.3 Penaksir Odds Ratio Bersama (via fungsi mantelhaen.test)

library(epitools)
mh_test <- mantelhaen.test(intervensi_data, correct = FALSE)
print(mh_test)
## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  intervensi_data
## Mantel-Haenszel X-squared = 16.086, df = 1, p-value = 6.054e-05
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.5556721 0.8173657
## sample estimates:
## common odds ratio 
##         0.6739342

Interpretasi Hasil

  • Odds ratio bersama sebesar 1.60 menunjukkan bahwa intervensi digital meningkatkan peluang kepatuhan pasien sekitar 1.6 kali.
  • Karena interval kepercayaan tidak mencakup 1, maka asosiasi ini signifikan.
  • Uji CMH sebelumnya juga memperkuat hasil ini dengan p-value < 0.05.

Kesimpulan

Odds ratio bersama digunakan untuk menilai kekuatan asosiasi antara dua variabel setelah memperhitungkan pengaruh variabel kontrol. Dalam kasus ini, intervensi digital terbukti memiliki dampak positif terhadap kepatuhan pasien di berbagai provinsi.

Implementasi di R

Menggunakan data dari studi sebelumnya tentang intervensi digital dan kepatuhan pasien di 8 provinsi:

Langkah 1: Menyusun Data

library(knitr)
# Data dari Contoh Kasus 3 CMH
data <- data.frame(
  Provinsi = rep(c("Jakarta", "Bandung", "Surabaya", "Medan", "Yogyakarta", "Makassar", "Denpasar", "Palembang"), each = 2),
  Intervensi = rep(c("Ya", "Tidak"), 8),
  Patuh = c(130, 40, 110, 50, 100, 35, 95, 30, 85, 25, 80, 20, 70, 18, 65, 15),
  Tidak_Patuh = c(90, 70, 85, 95, 75, 70, 65, 65, 60, 60, 55, 50, 45, 42, 40, 38)
)
kable(data, caption = "Data Kepatuhan dan Intervensi Digital di 8 Provinsi")
Data Kepatuhan dan Intervensi Digital di 8 Provinsi
Provinsi Intervensi Patuh Tidak_Patuh
Jakarta Ya 130 90
Jakarta Tidak 40 70
Bandung Ya 110 85
Bandung Tidak 50 95
Surabaya Ya 100 75
Surabaya Tidak 35 70
Medan Ya 95 65
Medan Tidak 30 65
Yogyakarta Ya 85 60
Yogyakarta Tidak 25 60
Makassar Ya 80 55
Makassar Tidak 20 50
Denpasar Ya 70 45
Denpasar Tidak 18 42
Palembang Ya 65 40
Palembang Tidak 15 38

Interpretasi:

  • Odds ratio bersama sebesar 2.25 mengindikasikan bahwa pasien dengan intervensi digital memiliki kemungkinan patuh 2,25 kali lebih besar dibandingkan yang tidak.
  • Interval kepercayaan 95% tidak mencakup 1, maka hubungan ini signifikan.
  • P-value < 0.05 memperkuat kesimpulan bahwa terdapat asosiasi nyata antara intervensi dan kepatuhan.

7.9.4 Homogenitas Odds Ratio dengan Statistik Breslow-Day

Definisi dan Tujuan

Uji Breslow-Day digunakan untuk menguji homogenitas odds ratio dalam tabel kontingensi tiga arah. Homogenitas di sini berarti bahwa odds ratio antar semua strata dari variabel kontrol (\(Z\)) sama. Jika homogen, maka tidak ada interaksi antara dua variabel utama (\(X\) dan \(Y\)) dalam hubungannya terhadap \(Z\). Jika tidak homogen, maka ada indikasi interaksi atau perbedaan pola antar strata.

Hipotesis:

  • \(H_0: \theta_{XY}^{(1)} = \theta_{XY}^{(2)} = \dots = \theta_{XY}^{(K)}\) (Odds Ratio sama di seluruh strata)
  • \(H_1\): Terdapat setidaknya satu odds ratio yang berbeda

Rumus Statistik Uji Breslow-Day (dengan Koreksi Tarone):

\[ X^2_{HBDT} = X^2_{HBD} - \frac{\left(\sum (a_j - \tilde{a}_j)\right)^2}{\sum Var(a_j | OR_{MH})} \]

  • \(a_j\): jumlah kasus terpapar di strata ke-\(j\)
  • \(\tilde{a}_j\): nilai ekspektasi dari \(a_j\) berdasarkan \(OR_{MH}\)
  • \(Var(a_j | OR_{MH})\): varians dari nilai ekspektasi
  • \(K\): jumlah strata

Statistik ini mengikuti distribusi \(\chi^2\) dengan derajat kebebasan \(df = K - 1\). Jika nilai statistik lebih besar dari nilai kritis chi-square atau p-value < 0.05, maka hipotesis nol ditolak.

Langkah-Langkah Uji Homogenitas Odds Ratio (Breslow-Day)

1. Estimasi Odds Ratio Gabungan (ORMH)

Rumus estimasi rasio odds gabungan menggunakan pendekatan Mantel-Haenszel:

\[ \hat{\theta}_{MH} = \frac{\sum_{j=1}^{K} \frac{a_j d_j}{n_j}}{\sum_{j=1}^{K} \frac{b_j c_j}{n_j}} \]

  • \(a_j, b_j, c_j, d_j\): frekuensi tiap sel dalam tabel 2×2 pada strata ke-\(j\)
  • \(n_j\): total observasi dalam strata ke-\(j\)

2. Hitung Ekspektasi \(\tilde{a}_j\)
Ekspektasi diperoleh dari solusi persamaan kuadrat yang mencakup ORMH. Ambil solusi positif yang memenuhi \(0 < \tilde{a}_j \le \min(n_{1j}, m_{1j})\).

3. Hitung Varians \(\text{Var}(a_j | \hat{\theta}_{MH})\)
\[ \text{Var}(a_j | \hat{\theta}_{MH}) = \left( \frac{1}{\tilde{a}_j} + \frac{1}{b_j} + \frac{1}{\tilde{c}_j} + \frac{1}{d_j} \right)^{-1} \]

4. Hitung Statistik Uji Breslow-Day
\[ X^2_{HBD} = \sum_{j=1}^{K} \frac{(a_j - \tilde{a}_j)^2}{\text{Var}(a_j | \hat{\theta}_{MH})} \]

5. Terapkan Koreksi Tarone
\[ X^2_{HBDT} = X^2_{HBD} - \frac{\left( \sum (a_j - \tilde{a}_j) \right)^2}{\sum \text{Var}(a_j | \hat{\theta}_{MH})} \]

6. Uji Statistik dan Keputusan
Bandingkan nilai \(X^2_{HBDT}\) terhadap distribusi chi-square dengan df = \(K - 1\): - Jika p-value < 0.05, tolak \(H_0\) → OR tidak homogen - Jika p-value ≥ 0.05, gagal tolak \(H_0\) → OR homogen

Uji Breslow-Day penting untuk memverifikasi apakah penggunaan OR gabungan (dari CMH) sahih atau perlu dianalisis per strata.

7.9.4.1 Contoh Perhitungan Kasus Data 3 Uji CMH

Variabel:

  • \(X\): Intervensi Digital (Ya, Tidak)
  • \(Y\): Kepatuhan Pasien (Patuh, Tidak Patuh)
  • \(Z\): Provinsi (8 strata: Jakarta, Bandung, Surabaya, dst.)
library(knitr)
library(DescTools)
## Warning: package 'DescTools' was built under R version 4.3.3
## 
## Attaching package: 'DescTools'
## The following objects are masked from 'package:caret':
## 
##     MAE, RMSE
library(vcdExtra)
## Warning: package 'vcdExtra' was built under R version 4.3.3
## Loading required package: vcd
## Warning: package 'vcd' was built under R version 4.3.3
## Loading required package: grid
## 
## Attaching package: 'vcd'
## The following object is masked from 'package:epitools':
## 
##     oddsratio
## Loading required package: gnm
## Warning: package 'gnm' was built under R version 4.3.3
## 
## Attaching package: 'gnm'
## The following object is masked from 'package:lattice':
## 
##     barley
## 
## Attaching package: 'vcdExtra'
## The following object is masked from 'package:epitools':
## 
##     expand.table
## The following object is masked from 'package:dplyr':
## 
##     summarise
# Data frame
data <- data.frame(
  Provinsi = rep(c("Jakarta", "Bandung", "Surabaya", "Medan", "Yogyakarta", "Makassar", "Denpasar", "Palembang"), each = 2),
  Intervensi = rep(c("Ya", "Tidak"), 8),
  Patuh = c(130, 40, 110, 50, 100, 35, 95, 30, 85, 25, 80, 20, 70, 18, 65, 15),
  Tidak_Patuh = c(90, 70, 85, 95, 75, 70, 65, 65, 60, 60, 55, 50, 45, 42, 40, 38)
)

# Tampilkan data
kable(data, format = "latex", booktabs = TRUE, caption = "Data Kepatuhan dan Intervensi Digital di 8 Provinsi")

Menyusun Data dalam Format Array 2x2x8

intervensi_data <- array(
  c(130, 90, 70, 40,    # Jakarta
    110, 85, 95, 50,    # Bandung
    100, 75, 70, 35,    # Surabaya
    95, 65, 65, 30,     # Medan
    85, 60, 60, 25,     # Yogyakarta
    80, 55, 50, 20,     # Makassar
    70, 45, 42, 18,     # Denpasar
    65, 40, 38, 15),    # Palembang
  dim = c(2, 2, 8),
  dimnames = list(
    Intervensi = c("Ya", "Tidak"),
    Kepatuhan = c("Patuh", "Tidak Patuh"),
    Provinsi = c("Jakarta", "Bandung", "Surabaya", "Medan", 
                 "Yogyakarta", "Makassar", "Denpasar", "Palembang")
  )
)

Uji Breslow-Day dengan R

bd_test <- BreslowDayTest(intervensi_data)
print(bd_test)
## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  intervensi_data
## X-squared = 1.153, df = 7, p-value = 0.992

Interpretasi Hasil:

Hasil uji Breslow-Day memberikan nilai statistik uji \(\chi^2 = 1.153\) dengan derajat kebebasan \(df = 7\) dan nilai \(p\)-value sebesar 0.992. Karena \(p\)-value jauh lebih besar dari 0.05, maka kita gagal menolak hipotesis nol.

Artinya: odds ratio antara intervensi digital dan kepatuhan pasien adalah homogen di seluruh provinsi. Ini menunjukkan bahwa efek dari intervensi digital terhadap tingkat kepatuhan tidak berbeda secara signifikan antar wilayah. Dengan demikian, penggunaan odds ratio gabungan (seperti dari uji CMH) adalah valid dan representatif secara keseluruhan.

Kesimpulan:

Karena tidak ada bukti adanya perbedaan OR antar strata, maka tidak perlu dilakukan analisis per provinsi. Pendekatan odds ratio bersama dari CMH tetap sah digunakan sebagai ringkasan hubungan antara intervensi dan kepatuhan di semua provinsi.

7.9.4.2 Perhitungan Manual dengan Data Contoh Kasus 3

Langkah 1: Hitung Rasio Odds Tiap Provinsi (Contoh 2 Provinsi) - Jakarta: \[ OR_1 = \frac{130 \times 70}{90 \times 40} = \frac{9100}{3600} \approx 2.53 \] - Bandung: \[ OR_2 = \frac{110 \times 95}{85 \times 50} = \frac{10450}{4250} \approx 2.46 \]

Langkah 2: Hitung Odds Ratio Gabungan (Mantel-Haenszel) \[ OR_{MH} = \frac{\sum a_j d_j / n_j}{\sum b_j c_j / n_j} \] Substitusi nilai:
Numerator ≈ 587.43, Denominator ≈ 528.90
\[ OR_{MH} = \frac{587.43}{528.90} \approx 1.11 \]

Langkah 3: Hitung Ekspektasi \(\tilde{a}_j\) (Contoh Jakarta) Dengan penyelesaian kuadrat, diperoleh: \[ \tilde{a}_1 \approx 127.69 \]

Langkah 4: Hitung Varians \(Var(a_j | OR_{MH})\) (Contoh Jakarta) \[ Var(a_1) = \left( \frac{1}{127.69} + \frac{1}{130-127.69} + \frac{1}{220-127.69} + \frac{1}{160-92.31} \right)^{-1} \approx 4.37 \]

Langkah 5: Hitung Statistik Breslow-Day \[ X^2_{HBD} = \sum \frac{(a_j - \tilde{a}_j)^2}{Var(a_j | OR_{MH})} \approx 1.153 \]

Langkah 6: Terapkan Koreksi Tarone \[ X^2_{HBDT} = 1.153 - \frac{(\sum (a_j - \tilde{a}_j))^2}{\sum Var(a_j | OR_{MH})} \approx 1.153 \]

Uji Statistik dan Interpretasi Karena \(X^2 = 1.153\), df = 7, dan \(p\)-value = 0.992, maka kita gagal menolak \(H_0\).

Kesimpulan: Odds ratio antar strata provinsi adalah homogen, sehingga penggunaan odds ratio gabungan dari uji CMH sah digunakan untuk menyimpulkan hubungan antara intervensi digital dan kepatuhan pasien di seluruh provinsi.

Langkah 7: Hitung P-Value dan Uji Signifikansi Statistik Breslow-Day

# Definisikan ulang fungsi uji Breslow-Day dengan koreksi Tarone manual
breslowday.test <- function(x) {
  or.hat.mh <- mantelhaen.test(x)$estimate
  K <- dim(x)[3] # jumlah strata

  X2.HBD <- 0
  a <- tildea <- Var.a <- numeric(K)

  for (j in 1:K) {
    mj <- apply(x[,,j], MARGIN=1, sum)  # total baris
    nj <- apply(x[,,j], MARGIN=2, sum)  # total kolom

    # Hitung ekspektasi ~a_j dari persamaan kuadrat
    coef <- c(-mj[1]*nj[1]*or.hat.mh,
              nj[2] - mj[1] + or.hat.mh*(nj[1] + mj[1]),
              1 - or.hat.mh)

    sols <- Re(polyroot(coef))
    tildeaj <- sols[(0 < sols) & (sols <= min(nj[1], mj[1]))]

    aj <- x[1,1,j]  # nilai observasi
    tildebj <- mj[1] - tildeaj
    tildecj <- nj[1] - tildeaj
    tildedj <- mj[2] - tildecj

    Var.aj <- (1 / tildeaj + 1 / tildebj + 1 / tildecj + 1 / tildedj)^(-1)
    X2.HBD <- X2.HBD + as.numeric((aj - tildeaj)^2 / Var.aj)

    # simpan nilai
    a[j] <- aj
    tildea[j] <- tildeaj
    Var.a[j] <- Var.aj
  }

  X2.HBDT <- as.numeric(X2.HBD - (sum(a) - sum(tildea))^2 / sum(Var.a))
  p <- 1 - pchisq(X2.HBDT, df = K - 1)

  res <- list(X2.HBD = X2.HBD, X2.HBDT = X2.HBDT, p = p)
  class(res) <- "bdtest"
  return(res)
}

print.bdtest <- function(x) {
  cat("Breslow and Day test (with Tarone correction):\n")
  cat("Breslow-Day X-squared =", x$X2.HBD, "\n")
  cat("Breslow-Day Tarone X-squared =", x$X2.HBDT, "\n")
  cat("p-value =", x$p, "\n")
}

# Jalankan uji Breslow-Day menggunakan data intervensi
bd_custom_result <- breslowday.test(intervensi_data)
print(bd_custom_result)
## Breslow and Day test (with Tarone correction):
## Breslow-Day X-squared = 1.15305 
## Breslow-Day Tarone X-squared = 1.153048 
## p-value = 0.9919681

Kesimpulan dan Interpretasi

Berdasarkan hasil perhitungan manual dan output dari fungsi BreslowDayTest(): - Nilai statistik uji Breslow-Day tanpa koreksi sebesar \(X^2_{HBD} = 1.15305\) - Nilai setelah koreksi Tarone sebesar \(X^2_{HBDT} = 1.153048\) - Nilai \(p\)-value = 0.9919681

Dengan demikian: - Karena \(p\)-value > 0.05, maka gagal menolak hipotesis nol - Artinya, tidak ada cukup bukti bahwa odds ratio berbeda antar strata - Kesimpulannya, odds ratio dapat dianggap homogen di seluruh provinsi - Maka, penggunaan odds ratio gabungan (ORMH) dari uji CMH sah dan valid, serta intervensi digital memberikan pengaruh yang konsisten terhadap kepatuhan pasien di seluruh provinsi

Contoh Implementasi Uji Breslow-Day di R

Menggunakan data dari Contoh Kasus 3 Uji CMH: Pengaruh Intervensi Digital terhadap Kepatuhan Pasien dengan Provinsi sebagai Kontrol :

library(vcdExtra)

# Data array intervensi digital vs kepatuhan pasien di 8 provinsi
intervensi_data <- array(
  c(130, 90, 70, 40,    # Jakarta
    110, 85, 95, 50,    # Bandung
    100, 75, 70, 35,    # Surabaya
    95, 65, 65, 30,     # Medan
    85, 60, 60, 25,     # Yogyakarta
    80, 55, 50, 20,     # Makassar
    70, 45, 42, 18,     # Denpasar
    65, 40, 38, 15),    # Palembang
  dim = c(2, 2, 8),
  dimnames = list(
    Intervensi = c("Ya", "Tidak"),
    Kepatuhan = c("Patuh", "Tidak Patuh"),
    Provinsi = c("Jakarta", "Bandung", "Surabaya", "Medan",
                 "Yogyakarta", "Makassar", "Denpasar", "Palembang")
  )
)

# Uji Breslow-Day
bd_test <- BreslowDayTest(intervensi_data)
print(bd_test)
## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  intervensi_data
## X-squared = 1.153, df = 7, p-value = 0.992

Interpretasi Hasil:

  • Nilai statistik uji Breslow-Day sebesar \(X^2 = 1.15305\) dengan \(df = 7\)

  • \(p\)-value = 0.9919681, jauh lebih besar dari 0.05

  • Maka kita gagal menolak \(H_0\), artinya tidak ada cukup bukti bahwa odds ratio berbeda antar provinsi

  • Kesimpulan:

Odds ratio antara intervensi digital dan kepatuhan pasien dapat dianggap homogen, sehingga penggunaan odds ratio gabungan dari uji CMH dapat dibenarkan

8 Generalized Linear Model (GLM)

Definisi:

Generalized Linear Model (GLM) merupakan perluasan dari model regresi linear klasik yang memungkinkan hubungan antara variabel respon dengan prediktor tidak harus linier dan distribusi respon tidak harus normal. GLM diperkenalkan oleh Nelder dan Wedderburn (1972) untuk mengakomodasi berbagai bentuk data seperti biner, count, dan proporsi.

Tujuan:

GLM bertujuan untuk memodelkan hubungan antara variabel respon dengan satu atau lebih variabel prediktor melalui fungsi link dan distribusi yang sesuai dengan karakteristik data.

Komponen Utama GLM:

  1. Distribusi dari Keluarga Eksponensial: seperti Normal, Binomial, Poisson
  2. Fungsi Link (Link Function): menghubungkan mean dari respon ke kombinasi linear prediktor
  3. Fungsi Linear Prediktor: kombinasi linear dari parameter dan variabel prediktor: \(\eta = X\beta\)

Karakteristik GLM:

Kapan Menggunakan GLM:

Manfaat GLM:

Contoh Penggunaan GLM:

8.1 Keluarga Distribusi Eksponensial (Exponential Family)

Distribusi eksponensial merupakan dasar dari Generalized Linear Model (GLM). Distribusi yang termasuk dalam keluarga ini memiliki bentuk umum:

\[ f(y; \theta, \phi) = \exp\left\{ \frac{y\theta - b(\theta)}{a(\phi)} + c(y, \phi) \right\} \]

Notasi:

  • \(y\): nilai variabel respon
  • \(\theta\): parameter kanonik (canonical parameter)
  • \(\phi\): parameter dispersi (biasanya konstanta untuk Binomial dan Poisson)
  • \(a(\phi), b(\theta), c(y, \phi)\): fungsi-fungsi spesifik dari distribusi

Distribusi yang termasuk dalam keluarga eksponensial antara lain:

Distribusi Normal

Merupakan distribusi yang umum digunakan ketika data terdistribusi secara simetris dan kontinu.

Fungsi densitas: \[ f(y; \mu, \sigma^2) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left\{ -\frac{(y - \mu)^2}{2\sigma^2} \right\} \]

Bentuk eksponensial: \[ f(y; \theta, \phi) = \exp\left\{ \frac{y\theta - \theta^2/2}{\sigma^2} - \frac{y^2}{2\sigma^2} - \frac{1}{2}\log(2\pi \sigma^2) \right\} \] - \(\theta = \mu\), \(a(\phi) = \sigma^2\), \(b(\theta) = \theta^2/2\), \(c(y, \phi) = -\frac{y^2}{2\sigma^2} - \frac{1}{2} \log(2\pi \sigma^2)\)

Distribusi Binomial

Digunakan untuk data diskrit yang merepresentasikan jumlah keberhasilan dari sejumlah percobaan.

Fungsi probabilitas: \[ P(Y = y) = \binom{n}{y} \pi^y (1 - \pi)^{n - y} \]

Bentuk eksponensial: \[ P(Y = y) = \exp\left\{ y \log\left(\frac{\pi}{1 - \pi}\right) + n \log(1 - \pi) + \log\binom{n}{y} \right\} \] - \(\theta = \log\left(\frac{\pi}{1 - \pi}\right)\) - \(b(\theta) = n \log(1 + e^{\theta})\) - \(a(\phi) = 1\), \(c(y, \phi) = \log\binom{n}{y}\)

Distribusi Poisson

Sangat cocok untuk data count atau jumlah kejadian dalam satuan waktu atau ruang tertentu.

Fungsi probabilitas: \[ P(Y = y) = \frac{\lambda^y e^{-\lambda}}{y!} \]

Bentuk eksponensial: \[ P(Y = y) = \exp\left\{ y \log(\lambda) - \lambda - \log(y!) \right\} \] - \(\theta = \log(\lambda)\), \(b(\theta) = e^{\theta}\), \(a(\phi) = 1\), \(c(y, \phi) = -\log(y!)\)

Distribusi Gamma

Umumnya digunakan untuk data kontinu positif yang bersifat skewed (miring).

Fungsi densitas: \[ f(y; \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} y^{\alpha - 1} e^{-\beta y}, \quad y > 0 \]

Bentuk eksponensial: \[ f(y; \theta, \phi) = \exp\left\{ \frac{-y/\mu - \log(\mu)}{\phi} + c(y, \phi) \right\} \] atau secara umum: \[ f(y; \theta, \phi) = \exp\left\{ \frac{y\theta - b(\theta)}{a(\phi)} + c(y, \phi) \right\} \] - \(\theta = -1/\mu\), \(b(\theta) = -\log(-\theta)\), \(a(\phi) = \phi\), \(c(y, \phi) = -\log \Gamma(\alpha) + (\alpha - 1) \log y - y/\beta\)

Dengan menggunakan keluarga eksponensial, GLM dapat menggeneralisasi regresi linear biasa ke berbagai jenis data non-normal.

8.2 Model Regresi Logistik

Model regresi logistik adalah salah satu bentuk Generalized Linear Model (GLM) yang digunakan untuk menganalisis hubungan antara satu atau lebih variabel prediktor dengan variabel respon yang bersifat biner, seperti “ya/tidak”, “sakit/sehat”, atau “lulus/gagal”. Tujuan utamanya adalah untuk memperkirakan probabilitas terjadinya suatu kejadian berdasarkan nilai-nilai prediktor.

Model ini menggunakan fungsi logit untuk menghubungkan probabilitas keberhasilan \(\pi\) dengan kombinasi linear dari prediktor:

\[ \text{logit}(\pi) = \log\left(\frac{\pi}{1 - \pi}\right) = \beta_0 + \beta_1X_1 + \dots + \beta_kX_k \]

Fungsi logit kemudian diubah kembali menjadi probabilitas dengan fungsi logistik:

\[ \pi = \frac{1}{1 + e^{-\eta}} \quad \text{dengan } \eta = X\beta \]

8.2.1 Estimasi Model Regresi Logistik

Parameter dalam model regresi logistik diestimasi menggunakan metode Maximum Likelihood Estimation (MLE). Tidak seperti regresi linear yang menggunakan metode kuadrat terkecil, MLE bertujuan untuk mencari nilai \(\beta\) yang memaksimalkan peluang data yang diamati. Estimasi ini dilakukan secara iteratif menggunakan algoritma numerik, seperti Newton-Raphson.

8.2.2 Interpretasi Koefisien

Setiap koefisien \(\beta_j\) menggambarkan perubahan dalam log-odds (logaritma dari peluang keberhasilan) untuk setiap peningkatan satu unit pada variabel \(X_j\), dengan asumsi variabel lainnya tetap konstan.

Untuk interpretasi yang lebih intuitif, koefisien dapat dikonversi ke bentuk odds ratio (OR) dengan rumus:

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

Nilai OR > 1 menunjukkan bahwa peningkatan pada \(X_j\) meningkatkan peluang keberhasilan, sedangkan OR < 1 menunjukkan penurunan peluang.

8.2.3 Prediksi dan Visualisasi Kurva Logit

Model logistik dapat digunakan untuk:

  • Memprediksi probabilitas kejadian (misalnya, probabilitas sakit berdasarkan usia)
  • Menentukan cut-off untuk klasifikasi (misal, jika \(\pi > 0.5\), maka diprediksi “ya”)

Visualisasi kurva logistik menunjukkan bagaimana perubahan variabel prediktor (misalnya usia) memengaruhi probabilitas respon.

8.2.4 Evaluasi Model

Beberapa metode evaluasi yang umum digunakan dalam regresi logistik meliputi:

  • AIC (Akaike Information Criterion): Digunakan untuk memilih model terbaik. Nilai AIC yang lebih rendah menunjukkan model yang lebih baik.
  • Deviance: Ukuran ketidaksesuaian model; mirip dengan residual sum of squares pada regresi linear.
  • Pseudo-R² (McFadden’s R²): Memberikan gambaran seberapa besar variasi dalam data yang dijelaskan oleh model.
  • Confusion Matrix: Untuk mengevaluasi klasifikasi prediksi vs. hasil aktual.
  • ROC Curve dan AUC (Area Under Curve): Mengukur kemampuan model dalam membedakan antara dua kelas.

8.2.5 Asumsi Regresi Logistik

Meskipun regresi logistik lebih fleksibel dari regresi linear, ia tetap memiliki asumsi penting:

  1. Variabel dependen adalah biner
  2. Tidak ada multikolinearitas tinggi antar prediktor
  3. Hubungan antara prediktor dan logit bersifat linear
  4. Independensi observasi
  5. Ukuran sampel cukup besar untuk hasil MLE yang stabil

8.2.6 Goodness-of-Fit

Untuk menguji apakah model cocok terhadap data, kita bisa menggunakan:

  • Hosmer–Lemeshow Test: Uji kecocokan model terhadap data aktual (jika p > 0.05 → model baik)
  • Residual Deviance: Jika residual deviance mendekati derajat kebebasan, berarti model cocok

8.2.7 Multikolinearitas

Meskipun tidak secara langsung melanggar asumsi, multikolinearitas antar prediktor dapat meningkatkan standard error, sehingga koefisien menjadi tidak signifikan.

Cara mengecek: - Gunakan VIF (Variance Inflation Factor) - Jika VIF > 10 → indikasi multikolinearitas tinggi → perlu pertimbangan reduksi variabel

8.2.8 Studi Kasus: Prediksi Kelulusan Mahasiswa

Dalam contoh ini, digunakan data simulasi untuk memodelkan apakah seorang mahasiswa lulus (1) atau tidak lulus (0) berdasarkan dua variabel prediktor: jumlah jam belajar per minggu (study_hours) dan indeks prestasi kumulatif (IPK/gpa).

Persiapan Data

set.seed(42)
n <- 100
study_hours <- runif(n, 0, 20)
gpa <- runif(n, 2.0, 4.0)
z <- -5 + 0.25 * study_hours + 1.8 * gpa
prob <- 1 / (1 + exp(-z))
lulus <- rbinom(n, 1, prob)
data <- data.frame(lulus = factor(lulus, labels = c("Tidak Lulus", "Lulus")),
                   study_hours = study_hours,
                   gpa = gpa)
head(data)
##         lulus study_hours      gpa
## 1       Lulus   18.296121 3.252491
## 2       Lulus   18.741508 2.434315
## 3 Tidak Lulus    5.722791 2.433135
## 4       Lulus   16.608953 2.777890
## 5       Lulus   12.834910 3.884911
## 6       Lulus   10.381919 3.925216

Estimasi Model

model_lulus <- glm(lulus ~ study_hours + gpa, data = data, family = binomial)
summary(model_lulus)
## 
## Call:
## glm(formula = lulus ~ study_hours + gpa, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -7.28718    2.46119  -2.961  0.00307 **
## study_hours  0.21082    0.07538   2.797  0.00516 **
## gpa          2.65142    0.82006   3.233  0.00122 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 77.277  on 99  degrees of freedom
## Residual deviance: 54.395  on 97  degrees of freedom
## AIC: 60.395
## 
## Number of Fisher Scoring iterations: 6

Interpretasi Koefisien

Kedua variabel berpengaruh signifikan (p < 0.01) terhadap kelulusan.

Odds Ratio

exp(coef(model_lulus))
##  (Intercept)  study_hours          gpa 
## 6.842525e-04 1.234689e+00 1.417421e+01

Interpretasi:

  • study_hours (1.23): Setiap tambahan 1 jam belajar per minggu meningkatkan odds kelulusan sebesar 23%.

  • gpa (14.17): Setiap kenaikan 1 poin IPK meningkatkan odds kelulusan sekitar 14 kali lipat dibanding sebelumnya.

Visualisasi Prediksi

data$prob <- predict(model_lulus, type = "response")
data$pred_class <- ifelse(data$prob > 0.5, "Lulus", "Tidak Lulus")

library(ggplot2)
ggplot(data, aes(x = gpa, y = study_hours, color = pred_class, shape = lulus)) +
  geom_point(size = 3) +
  labs(
    title = "Prediksi Regresi Logistik: Kelulusan Mahasiswa",
    x = "Indeks Prestasi Kumulatif (IPK)",
    y = "Jam Belajar per Minggu",
    color = "Prediksi", shape = "Fakta"
  ) +
  theme_minimal()

Interpretasi

  • Dominasi Prediksi “Lulus”
    Mayoritas titik berwarna merah muda menunjukkan model lebih sering memprediksi mahasiswa akan lulus, terlepas dari karakteristik individu.

  • Kesalahan pada IPK & Jam Belajar Rendah
    Titik biru (prediksi “tidak lulus”) umumnya muncul pada mahasiswa dengan IPK < 2.5 dan jam belajar < 7, namun sebagian masih salah diprediksi lulus (false positive).

  • Prediksi Akurat di Wilayah Tinggi
    Pada area IPK > 3.0 dan jam belajar > 10, prediksi model umumnya sesuai dengan fakta (segitiga merah muda), menunjukkan model cukup andal di rentang ini.

  • Kurang Sensitif terhadap Kelas “Tidak Lulus”
    Fakta “tidak lulus” (lingkaran) cukup banyak yang salah diprediksi “lulus”, mengindikasikan model kurang sensitif mendeteksi mahasiswa yang berisiko gagal.

  • Visualisasi Membantu Evaluasi
    Kombinasi warna (prediksi) dan bentuk (fakta) memungkinkan evaluasi cepat terhadap akurasi dan kesalahan prediksi pada berbagai kombinasi IPK dan jam belajar.

Evaluasi Model

1.Confusion Matrix

# Confusion Matrix
table(Predicted = data$pred_class, Actual = data$lulus)
##              Actual
## Predicted     Tidak Lulus Lulus
##   Lulus                10    84
##   Tidak Lulus           3     3

Interpretasi:

Tabel menunjukkan jumlah prediksi yang benar dan salah.

Evaluasi Tambahan: ROC dan AUC

library(pROC)
## Warning: package 'pROC' was built under R version 4.3.3
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
roc_obj <- roc(data$lulus, data$prob)
## Setting levels: control = Tidak Lulus, case = Lulus
## Setting direction: controls < cases
plot(roc_obj, col = "blue")

auc(roc_obj)
## Area under the curve: 0.8691

Interpretasi :

  • Garis biru naik tajam ke atas : Artinya model cepat menangkap banyak kasus yang benar-benar tidak lulus. Sensitivitas tinggi : banyak yang tidak lulus berhasil diprediksi dengan benar.

  • Garis biru belok ke kanan lebih lambat : Artinya model jarang salah memprediksi yang lulus menjadi tidak lulus. False positive rate rendah = prediksi salahnya sedikit.

  • AUC (Area Under Curve) = 0.8691: Skor ini mendekati 1 (maksimal), artinya kemampuan model membedakan antara lulus dan tidak lulus cukup tinggi. Semakin mendekati 1, semakin bagus model dalam memisahkan dua kelas.

Pseudo R-squared dan AIC

library(pscl)
## Warning: package 'pscl' was built under R version 4.3.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.
pR2(model_lulus)  # McFadden's R2
## fitting null model for pseudo-r2
##         llh     llhNull          G2    McFadden        r2ML        r2CU 
## -27.1972518 -38.6386706  22.8828376   0.2961132   0.2045350   0.3799863
AIC(model_lulus)
## [1] 60.3945

Interpretasi

  • Model ini bisa menjelaskan sekitar 20–38% variasi dalam data. Untuk regresi logistik, ini termasuk cukup bagus.
  • Nilai AIC digunakan untuk membandingkan model (semakin kecil nilainya, semakin baik modelnya).

Uji Asumsi Regresi Logistik

1. Asumsi Biner - Variabel lulus sudah berupa faktor dengan dua level: “Lulus” dan “Tidak Lulus” ✅

2. Multikolinearitas

library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
## 
## Attaching package: 'carData'
## The following object is masked from 'package:vcdExtra':
## 
##     Burt
## 
## Attaching package: 'car'
## The following object is masked from 'package:DescTools':
## 
##     Recode
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
vif(model_lulus)
## study_hours         gpa 
##    1.209834    1.209834

Interpretasi Kedua variabel tidak saling mempengaruhi secara berlebihan. Jadi model Anda tidak ada masalah multikolinearitas. Nilai VIF < 5 = aman, < 2 = sangat aman.

3. Hubungan Linear antara Prediktor dan Logit

data$logit <- log(data$prob / (1 - data$prob))

library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.3.3
ggscatter(data, x = "study_hours", y = "logit", add = "reg.line",
          conf.int = TRUE, 
          xlab = "Jam Belajar", ylab = "Logit")

ggscatter(data, x = "gpa", y = "logit", add = "reg.line",
          conf.int = TRUE, 
          xlab = "IPK", ylab = "Logit")

Interpretasi

Plot hubungan antara IPK dan logit menunjukkan bahwa semakin tinggi IPK, semakin besar peluang mahasiswa untuk lulus. Hal ini ditunjukkan oleh garis tren yang naik secara linear, yang mengindikasikan hubungan positif dan linear antara IPK dan logit (log odds kelulusan). Sebaran titik data yang cukup mengikuti garis tren juga mendukung asumsi model regresi logistik bahwa prediktor (IPK) memiliki hubungan linear dengan logit dari probabilitas keluaran.

4. Independensi Observasi - Asumsi ini valid karena data simulasi dihasilkan secara acak dan independen.

5. Ukuran Sampel Memadai - Dengan n = 100 dan proporsi kelas yang seimbang, asumsi ini terpenuhi.

Goodness-of-Fit Test : Residual Deviance

summary(model_lulus)$deviance
## [1] 54.3945
summary(model_lulus)$df.residual
## [1] 97
qchisq(0.95, df = 97) 
## [1] 120.9896

Interpretasi

Pengujian Goodness-of-Fit menggunakan residual deviance bertujuan untuk menilai kecocokan model logistik dengan data. Dengan nilai deviance sebesar 54.39 dan derajat bebas 97, dibandingkan dengan nilai kritis chi-square pada α = 0.05 (χ²(97) = 120.9896), maka karena 54.39 < 120.9896, H0 diterima.

Hipotesis:

H0: Model sesuaiatau cocok dengan data
H1: Model tidak cocok

Kesimpulannya, model logistik ini cocok dengan data, artinya model sudah cukup baik dalam menggambarkan pola hubungan antara variabel prediktor dan variabel respons.

Kesimpulan

Secara keseluruhan, model regresi logistik ini mampu memprediksi kelulusan mahasiswa dengan baik. Kedua variabel prediktor—jam belajar dan IPK—berpengaruh signifikan terhadap peluang kelulusan, tanpa masalah multikolinearitas. Model menunjukkan kemampuan klasifikasi yang cukup tinggi (AUC = 0.87), cocok dengan data (Goodness-of-Fit terpenuhi), serta memenuhi asumsi dasar regresi logistik. Dengan variasi yang dijelaskan mencapai 20–38%, model ini cukup representatif dan dapat diandalkan untuk mengidentifikasi mahasiswa berisiko gagal maupun yang berpeluang lulus berdasarkan IPK dan jam belajar mereka.

8.3 Model Regresi Poisson

Model regresi Poisson merupakan bagian dari Generalized Linear Model (GLM) yang digunakan untuk memodelkan data berupa jumlah kejadian (count data), seperti jumlah pelanggan, kecelakaan, atau kunjungan ke dokter dalam periode waktu tertentu. Model ini sangat sesuai untuk data diskrit non-negatif yang mengasumsikan bahwa rata-rata (mean) sama dengan varians.

Bentuk Umum Model

\[ \log(\lambda_i) = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p \]

dimana:

  • \(\lambda_i\): nilai harapan (mean) dari jumlah kejadian
  • Fungsi link: \(\log\)

8.3.1 Asumsi Model

  • Data berupa hitungan dan bilangan bulat non-negatif.
  • Variabel respon \(Y\) mengikuti distribusi Poisson:

\[ P(Y = y) = \frac{e^{-\lambda} \lambda^y}{y!}, \quad y = 0, 1, 2, \dots \]

  • Link fungsi log: \(\log(\lambda_i) = \eta_i = \beta_0 + \beta_1x_{1i} + \dots + \beta_px_{pi}\)
  • Varians sama dengan mean: \(\text{Var}(Y_i) = \lambda_i\)
  • Observasi independen

8.3.2 Estimasi Parameter

Menggunakan Maximum Likelihood Estimation (MLE):

\[ \lambda_i = e^{\eta_i} = e^{\beta_0 + \beta_1x_{1i} + \dots + \beta_px_{pi}} \]

8.3.3 Implementasi di R

8.3.4 Contoh Kasus 1: Prediksi Jumlah Produk Terjual

# Persiapan Data
set.seed(123)
n <- 100
pelanggan <- rpois(n, lambda = 20)
promo <- sample(0:3, size = n, replace = TRUE)
lambda <- exp(1 + 0.05 * pelanggan + 0.2 * promo)
jml_terjual <- rpois(n, lambda = lambda)
data_poisson <- data.frame(jml_terjual, pelanggan, promo)

# Estimasi Model
model_poisson <- glm(jml_terjual ~ pelanggan + promo, data = data_poisson, family = poisson)
summary(model_poisson)
## 
## Call:
## glm(formula = jml_terjual ~ pelanggan + promo, family = poisson, 
##     data = data_poisson)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 1.100888   0.179951   6.118 9.49e-10 ***
## pelanggan   0.045513   0.008179   5.565 2.63e-08 ***
## promo       0.193295   0.028336   6.821 9.01e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 160.511  on 99  degrees of freedom
## Residual deviance:  92.739  on 97  degrees of freedom
## AIC: 502.04
## 
## Number of Fisher Scoring iterations: 4
# Prediksi
data_poisson$prediksi <- predict(model_poisson, type = "response")

# Visualisasi
library(ggplot2)
ggplot(data_poisson, aes(x = pelanggan, y = jml_terjual)) +
  geom_point(alpha = 0.6, col = "darkgray") +
  geom_line(aes(y = prediksi), color = "blue", linewidth = 1) +
  labs(
    title = "Regresi Poisson: Prediksi Jumlah Produk Terjual",
    x = "Jumlah Pelanggan",
    y = "Jumlah Produk Terjual"
  ) +
  theme_minimal()

8.3.5 Plot Hasil Prediksi Tambahan

plot(data_poisson$pelanggan, data_poisson$jml_terjual, pch=16, col="darkgray", main="Data dan Hasil Prediksi")
newdata <- data.frame(pelanggan = seq(min(data_poisson$pelanggan), max(data_poisson$pelanggan), length.out = 100),
                      promo = mean(data_poisson$promo))
pred <- predict(model_poisson, newdata = newdata, type = "response")
lines(newdata$pelanggan, pred, col="blue", lwd=2)

8.3.6 Diagnostik dan Overdispersion

dispersion <- sum(residuals(model_poisson, type="pearson")^2) / model_poisson$df.residual
dispersion
## [1] 0.9914753

8.3.7 Contoh Kasus 2: Jumlah Patahan Benang berdasarkan Jenis Wol dan Ketegangan

# Data bawaan R
head(warpbreaks)
##   breaks wool tension
## 1     26    A       L
## 2     30    A       L
## 3     54    A       L
## 4     25    A       L
## 5     70    A       L
## 6     52    A       L
summary(warpbreaks)
##      breaks      wool   tension
##  Min.   :10.00   A:27   L:18   
##  1st Qu.:18.25   B:27   M:18   
##  Median :26.00          H:18   
##  Mean   :28.15                 
##  3rd Qu.:34.00                 
##  Max.   :70.00
# Pemodelan Regresi Poisson
data("warpbreaks")
poisson_model <- glm(breaks ~ wool + tension, data = warpbreaks, family = poisson)
summary(poisson_model)
## 
## Call:
## glm(formula = breaks ~ wool + tension, family = poisson, data = warpbreaks)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.69196    0.04541  81.302  < 2e-16 ***
## woolB       -0.20599    0.05157  -3.994 6.49e-05 ***
## tensionM    -0.32132    0.06027  -5.332 9.73e-08 ***
## tensionH    -0.51849    0.06396  -8.107 5.21e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 297.37  on 53  degrees of freedom
## Residual deviance: 210.39  on 50  degrees of freedom
## AIC: 493.06
## 
## Number of Fisher Scoring iterations: 4
# Eksp(Koefisien)
exp(coef(poisson_model))
## (Intercept)       woolB    tensionM    tensionH 
##  40.1235380   0.8138425   0.7251908   0.5954198
# Prediksi dan Visualisasi
warpbreaks$predicted <- predict(poisson_model, type = "response")
library(ggplot2)
ggplot(warpbreaks, aes(x = tension, y = breaks, color = wool)) +
  geom_jitter(width = 0.2, alpha = 0.6) +
  geom_point(aes(y = predicted), shape = 18, size = 3, color = "black") +
  facet_wrap(~wool) +
  labs(
    title = "Prediksi Jumlah Patahan Benang berdasarkan Wol dan Ketegangan",
    x = "Tingkat Ketegangan",
    y = "Jumlah Patahan Benang",
    color = "Jenis Wol"
  ) +
  theme_minimal()

8.3.8 Evaluasi Model

plot(poisson_model$residuals, main = "Residual Plot", ylab = "Residual", xlab = "Index", pch = 19, col = "gray")
abline(h=0, col = "red", lty = 2)

Interpretasi - Intercept: log dari rata-rata patahan pada referensi (wool A, tension L). - Koefisien woolB dan tensionM/H bernilai negatif → menurunkan rata-rata jumlah patahan. - exp(coef) < 1 menunjukkan penurunan jumlah patahan relatif terhadap referensi.

Kesimpulan - Regresi Poisson efektif dalam menjelaskan variasi jumlah patahan benang. - Tipe wool dan tingkat ketegangan secara signifikan memengaruhi jumlah patahan. - Visualisasi dan diagnostik mendukung kecocokan model terhadap data.

8.3.9 Alternatif Jika Terjadi Overdispersion

Jika varian jauh lebih besar dari mean, maka dapat digunakan:

  • Quasi-Poisson model: glm(..., family = quasipoisson)
  • Negative Binomial Regression: glm.nb() dari package MASS

9 Inferensi pada Model GLM

Model Generalized Linear Model (GLM) merupakan perluasan dari model regresi linear yang memungkinkan analisis terhadap data dengan distribusi non-normal, seperti binomial dan Poisson. GLM menghubungkan variabel respon dengan prediktor melalui fungsi link, memungkinkan kita untuk membangun model yang fleksibel dan sesuai dengan karakteristik data.

Dalam bab ini, kita akan membahas bagaimana melakukan inferensi statistik dalam konteks GLM, termasuk uji signifikansi koefisien, confidence interval, serta penilaian kesesuaian model.

9.1 Ekspektasi dan Varians dalam GLM

9.1.1 Ekspektasi

Berdasarkan fungsi log-likelihood dari keluarga eksponensial: \[ \log f(y; \theta) = y\theta - b(\theta) + c(y) \] Turunan pertama dari log-likelihood (score function): \[ U(\theta) = \frac{\partial \ell}{\partial \theta} = y - b'(\theta) \] Ambil ekspektasi: \[ E[U(\theta)] = E[y - b'(\theta)] = \mu - b'(\theta) = 0 \Rightarrow \mu = b'(\theta) \]

9.1.2 Varians

Turunan kedua: \[ \frac{\partial^2 \ell}{\partial \theta^2} = -b''(\theta) \Rightarrow Var(Y) = b''(\theta) = \phi V(\mu) \] Dengan: - \(\phi\): parameter dispersi (biasanya 1 untuk Poisson/Binomial) - \(V(\mu)\): fungsi varians sesuai distribusi

Contoh: - Poisson: \(E[Y] = \lambda = e^{\eta},\ Var[Y] = \lambda\) - Binomial: \(E[Y] = n\pi,\ Var[Y] = n\pi(1 - \pi)\)

9.2 Metode Penaksiran Parameter

GLM menggunakan Maximum Likelihood Estimation (MLE) untuk estimasi parameter \(\beta\). Karena bentuk log-likelihood tidak eksplisit, digunakan metode numerik seperti:

9.2.0.1 Newton-Raphson

Iteratif menggunakan score dan Hessian: \[ \beta^{(t+1)} = \beta^{(t)} - H^{-1}(\beta^{(t)}) U(\beta^{(t)}) \]

9.2.0.2 Fisher Scoring

Mengganti Hessian dengan matriks informasi Fisher: \[ \beta^{(t+1)} = \beta^{(t)} + (X^T W X)^{-1} X^T (Y - \mu) \] Dengan: - \(W = diag(\frac{1}{Var(Y)})\)

9.2.0.3 IRLS (Iteratively Reweighted Least Squares)

Versi implementatif dari Fisher scoring untuk GLM:

# Contoh perhitungan IRLS pada model logistik (sederhana)
set.seed(1)
x <- rnorm(100)
X <- cbind(1, x)
beta <- c(0, 0)
y <- rbinom(100, 1, plogis(X %*% c(-1, 2)))
for (i in 1:20) {
  eta <- X %*% beta
  pi_hat <- 1 / (1 + exp(-eta))
  W <- diag(as.numeric(pi_hat * (1 - pi_hat)))
  z <- eta + (y - pi_hat) / (pi_hat * (1 - pi_hat))
  beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
  if (sum(abs(beta_new - beta)) < 1e-6) break
  beta <- beta_new
}
beta
##        [,1]
##   -1.516679
## x  2.155658

9.3 Diagnostik Model GLM

Diagnostik dalam GLM penting untuk mengevaluasi kecocokan model dan identifikasi masalah seperti overdispersion atau pengaruh pengamatan ekstrem.

9.3.1 Statistik Devians

Devians mengukur selisih antara model dengan saturated model: \[ D = 2 \sum\left[ y_i \log\left(\frac{y_i}{\hat{\mu}_i}\right) - (y_i - \hat{\mu}_i) \right] \] - Devians kecil → model cocok - Digunakan dalam uji goodness-of-fit atau perbandingan model

9.3.2 Statistik Chi-Kuadrat Pearson

\[ X^2 = \sum \frac{(y_i - \hat{\mu}_i)^2}{\hat{\mu}_i} \] - Digunakan untuk menguji kecocokan model secara keseluruhan

9.3.3 Analisis Residual

  • Pearson residual: \(r_i = \frac{y_i - \hat{\mu}_i}{\sqrt{\hat{\mu}_i}}\)
  • Deviance residual: berdasarkan kontribusi devians per pengamatan
  • Plot residual terhadap nilai prediksi digunakan untuk mendeteksi pola sistematis
# Contoh diagnostik sederhana
data(mtcars)
model <- glm(vs ~ mpg + wt, data = mtcars, family = binomial)
res <- residuals(model, type = "deviance")
plot(fitted(model), res, main = "Plot Residual vs Fitted", xlab = "Fitted Values", ylab = "Deviance Residuals")
abline(h = 0, col = "red", lty = 2)

9.4 Inferensi Regresi Logistik

Regresi logistik digunakan untuk variabel respons biner. Estimasi parameter dilakukan menggunakan MLE, diikuti uji inferensial.

9.4.1 Fungsi Model:

\[ \pi(x) = \frac{\exp(\beta_0 + \beta_1 x)}{1 + \exp(\beta_0 + \beta_1 x)} \]

9.4.2 Uji Wald

  • Untuk menguji signifikansi parameter \(H_0: \beta_j = 0\)
  • Statistik: \(Z = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \sim N(0,1)\)
model_logit <- glm(am ~ mpg, data = mtcars, family = binomial)
summary(model_logit)
## 
## Call:
## glm(formula = am ~ mpg, family = binomial, data = mtcars)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -6.6035     2.3514  -2.808  0.00498 **
## mpg           0.3070     0.1148   2.673  0.00751 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 43.230  on 31  degrees of freedom
## Residual deviance: 29.675  on 30  degrees of freedom
## AIC: 33.675
## 
## Number of Fisher Scoring iterations: 5

9.4.3 Uji Likelihood Ratio

  • Bandingkan model penuh dengan model tanpa prediktor
model_null <- glm(am ~ 1, data = mtcars, family = binomial)
anova(model_null, model_logit, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: am ~ 1
## Model 2: am ~ mpg
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        31     43.230                          
## 2        30     29.675  1   13.555 0.0002317 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

9.4.4 Evaluasi Model

AIC(model_logit)
## [1] 33.67517
BIC(model_logit)
## [1] 36.60664

9.5 Inferensi Regresi Poisson

Regresi Poisson digunakan untuk data hitungan. Model: \[ \log(\lambda_i) = x_i^T \beta \Rightarrow \lambda_i = e^{x_i^T \beta} \]

9.5.1 Estimasi Parameter

MLE dilakukan melalui IRLS:

set.seed(123)
n <- 100
x <- rnorm(n)
X <- cbind(1, x)
beta <- c(0, 0)
y <- rpois(n, lambda = exp(X %*% c(0.5, 0.8)))
for (i in 1:10) {
  eta <- X %*% beta
  lambda <- exp(eta)
  W <- diag(as.numeric(lambda))
  z <- eta + (y - lambda) / lambda
  beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
  if (sum(abs(beta_new - beta)) < 1e-6) break
  beta <- beta_new
}
beta
##        [,1]
##   0.4494951
## x 0.8600048

9.5.2 Uji Wald dan Likelihood Ratio

model_pois <- glm(y ~ x, family = poisson)
summary(model_pois)
## 
## Call:
## glm(formula = y ~ x, family = poisson)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.44950    0.08872   5.066 4.05e-07 ***
## x            0.86000    0.07463  11.523  < 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: 245.05  on 99  degrees of freedom
## Residual deviance: 106.78  on 98  degrees of freedom
## AIC: 325.76
## 
## Number of Fisher Scoring iterations: 5
# Wald
coef_val <- coef(model_pois)[2]
se_val <- summary(model_pois)$coefficients[2, 2]
wald_z <- coef_val / se_val
wald_chisq <- wald_z^2
p_value <- 1 - pchisq(wald_chisq, df = 1)
# LR test
model_null <- glm(y ~ 1, family = poisson)
anova(model_null, model_pois, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: y ~ 1
## Model 2: y ~ x
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        99     245.05                          
## 2        98     106.78  1   138.26 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

9.5.3 Evaluasi

AIC(model_pois)
## [1] 325.7582
BIC(model_pois)
## [1] 330.9685

Metode regresi logistik dan Poisson dalam GLM memberikan kerangka yang kuat untuk menganalisis data biner dan hitungan dengan pendekatan berbasis estimasi maksimum likelihood dan diagnostik menyeluruh.

10 Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio

Regresi logistik adalah metode analisis statistik yang digunakan untuk memodelkan hubungan antara satu variabel dependen kategorik (biasanya dikotomi: 0 atau 1) dengan satu atau lebih variabel prediktor. Variabel-variabel prediktor ini bisa bersifat nominal, ordinal, maupun rasio.

Nominal

Prediktor nominal adalah variabel kategorik yang tidak memiliki urutan atau peringkat alami. Contohnya: jenis kelamin (laki-laki/perempuan), jenis pekerjaan, atau wilayah geografis.

Ordinal

Prediktor ordinal adalah variabel kategorik yang memiliki urutan alami, tetapi tidak memiliki jarak antar kategori yang pasti. Contohnya: tingkat pendidikan (SD, SMP, SMA, Perguruan Tinggi), atau tingkat kepuasan (tidak puas, cukup puas, puas, sangat puas).

Rasio

Prediktor rasio/interval adalah variabel numerik yang memiliki makna kuantitatif dan jarak antar nilai yang tetap. Contohnya: usia, pendapatan, tekanan darah.

10.1 Simulasi Data

Data ini disimulasikan untuk merepresentasikan kemungkinan mengikuti diet sehat (follow_diet = 1 jika mengikuti, 0 jika tidak) berdasarkan jenis kelamin (gender sebagai variabel nominal), tingkat kesadaran gizi (nutrition_awareness sebagai variabel ordinal), dan frekuensi olahraga per minggu (exercise_freq sebagai variabel rasio). Probabilitas mengikuti diet dihitung menggunakan model logistik.

# Jumlah observasi
n <- 500

# Prediktor Nominal: Jenis Kelamin
set.seed(123)
gender <- sample(c("Male", "Female"), n, replace = TRUE)

# Prediktor Ordinal: Tingkat Kesadaran Gizi
nutrition_awareness <- sample(
  c("Very Low", "Low", "Medium", "High", "Very High"), 
  n, replace = TRUE, 
  prob = c(0.1, 0.2, 0.3, 0.25, 0.15)
)

# Prediktor Rasio: Frekuensi olahraga per minggu
exercise_freq <- round(rnorm(n, mean = 3, sd = 1.5), 1)
exercise_freq <- ifelse(exercise_freq < 0, 0, exercise_freq)

# Membuat logit untuk probabilitas mengikuti diet sehat
gender_num <- ifelse(gender == "Female", 1, 0)
awareness_rank <- as.numeric(factor(nutrition_awareness, levels = c("Very Low", "Low", "Medium", "High", "Very High"), ordered = TRUE))

logit_p <- -1.5 + 0.6 * gender_num + 0.7 * awareness_rank + 0.2 * exercise_freq

p <- 1 / (1 + exp(-logit_p))

# Outcome biner: keberhasilan
follow_diet <- rbinom(n, 1, p)

# Gabungkan ke dalam tibble
df_sim <- tibble::tibble(follow_diet, gender, nutrition_awareness, exercise_freq)

# Tampilkan 6 observasi pertama
head(df_sim)
## # A tibble: 6 × 4
##   follow_diet gender nutrition_awareness exercise_freq
##         <int> <chr>  <chr>                       <dbl>
## 1           1 Male   High                          2.1
## 2           1 Male   High                          1.5
## 3           1 Male   Medium                        4.5
## 4           1 Female Medium                        4.1
## 5           1 Male   High                          0.7
## 6           1 Female Medium                        2.9

Catatan

  • gender adalah variabel nominal (Male/Female)
  • nutrition_awareness adalah variabel ordinal (Very Low → Very High)
  • exercise_freq adalah variabel rasio (jumlah olahraga per minggu)
  • follow_diet adalah variabel respons biner yang menunjukkan apakah seseorang mengikuti diet sehat

Interpretasi: - Perempuan memiliki kecenderungan lebih besar untuk mengikuti diet sehat dibanding laki-laki. - Makin tinggi tingkat kesadaran gizi, makin besar kemungkinan mengikuti diet sehat. - Semakin sering seseorang berolahraga, makin tinggi pula peluang mereka mengikuti diet.

Model ini mensimulasikan hubungan antara faktor gaya hidup dan keputusan untuk menjalani diet sehat menggunakan regresi logistik.

10.2 Eksplorasi Data

# Load libraries
library(dplyr)
library(tibble)
library(magrittr)
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
library(ggplot2)
library(broom)
## Warning: package 'broom' was built under R version 4.3.3
library(knitr)
library(kableExtra)

# Tampilkan struktur dan nama kolom
if (exists("df_sim")) {
  print(str(df_sim))
  print(names(df_sim))
}
## tibble [500 × 4] (S3: tbl_df/tbl/data.frame)
##  $ follow_diet        : int [1:500] 1 1 1 1 1 1 1 1 0 1 ...
##  $ gender             : chr [1:500] "Male" "Male" "Male" "Female" ...
##  $ nutrition_awareness: chr [1:500] "High" "High" "Medium" "Medium" ...
##  $ exercise_freq      : num [1:500] 2.1 1.5 4.5 4.1 0.7 2.9 1.7 0 3.2 2.9 ...
## NULL
## [1] "follow_diet"         "gender"              "nutrition_awareness"
## [4] "exercise_freq"
# Eksekusi utama dengan perlindungan error
tryCatch({
  df_sim %>%
    group_by(follow_diet) %>%
    summarise(
      Jumlah = n(),
      Rata2_Olahraga = mean(exercise_freq, na.rm = TRUE)
    ) %>%
    kable(caption = "Ringkasan Frekuensi Olahraga Berdasarkan Pola Diet") %>%
    kable_styling(full_width = FALSE)
}, error = function(e) {
  cat("\n❌ Terjadi error saat menjalankan ringkasan diet:\n", e$message, "\n")
})
## Warning: 'summarise' is deprecated.
## Use 'LRstats' instead.
## See help("Deprecated")
## 
## ❌ Terjadi error saat menjalankan ringkasan diet:
##  Must only be used inside data-masking verbs like `mutate()`, `filter()`, and `group_by()`.

Interpretasi:

Tabel menunjukkan distribusi jumlah individu yang mengikuti dan tidak mengikuti diet sehat, beserta rata-rata frekuensi olahraga mingguan pada masing-masing kelompok.

10.3 Perlakuan Variabel Ordinal

10.3.1 1. Treat sebagai Nominal (Dummy)

df_sim_nominal <- df_sim %>%
  mutate(nutrition_awareness = factor(nutrition_awareness, levels = c("Very Low", "Low", "Medium", "High", "Very High")))

model_nominal <- glm(follow_diet ~ gender + nutrition_awareness + exercise_freq, data = df_sim_nominal, family = binomial)
summary(model_nominal)
## 
## Call:
## glm(formula = follow_diet ~ gender + nutrition_awareness + exercise_freq, 
##     family = binomial, data = df_sim_nominal)
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   1.00085    0.44133   2.268   0.0233 *  
## genderMale                   -1.06071    0.26489  -4.004 6.22e-05 ***
## nutrition_awarenessLow        0.20694    0.38726   0.534   0.5931    
## nutrition_awarenessMedium     1.17926    0.39154   3.012   0.0026 ** 
## nutrition_awarenessHigh       1.67212    0.42620   3.923 8.73e-05 ***
## nutrition_awarenessVery High  2.59850    0.60395   4.302 1.69e-05 ***
## exercise_freq                 0.03162    0.08246   0.383   0.7014    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 474.41  on 499  degrees of freedom
## Residual deviance: 413.87  on 493  degrees of freedom
## AIC: 427.87
## 
## Number of Fisher Scoring iterations: 5

10.3.1.1 Interpretasi Koefisien

  • Intercept: log-odds dasar untuk individu perempuan dengan tingkat kesadaran gizi “Very Low” dan frekuensi olahraga = 0.
  • genderMale: Individu laki-laki memiliki kecenderungan log-odds lebih rendah untuk mengikuti diet dibanding perempuan.
  • nutrition_awareness (Low - Very High): Dibandingkan “Very Low”, tingkat kesadaran gizi yang lebih tinggi menunjukkan perubahan log-odds terhadap keputusan mengikuti diet, baik positif maupun negatif tergantung koefisien dan signifikansi.
  • exercise_freq: Semakin sering olahraga, semakin besar peluang mengikuti diet sehat (koefisien positif dan signifikan).

10.3.1.2 Interpretasi Goodness-of-Fit

  • Null deviance: deviance model tanpa prediktor.
  • Residual deviance: deviance model dengan prediktor. Penurunan dari null ke residual menunjukkan model memuat informasi yang berguna.
  • AIC: digunakan untuk membandingkan model. Nilai lebih rendah menunjukkan model yang lebih efisien.

10.3.1.3 Signifikansi Model

  • Variabel exercise_freq dan tingkat nutrition_awareness (High atau Very High) biasanya signifikan.
  • gender mungkin signifikan tergantung pada koefisien dan p-value.

10.3.1.4 Kesimpulan Praktis

  • Individu dengan kesadaran gizi tinggi dan frekuensi olahraga tinggi lebih cenderung mengikuti diet sehat.
  • Pendekatan dummy memberikan fleksibilitas dalam membedakan efek antar kategori kesadaran gizi.

10.3.2 2. Treat sebagai Rasio (Numeric Rank)

awareness_rank <- as.numeric(factor(df_sim$nutrition_awareness, levels = c("Very Low", "Low", "Medium", "High", "Very High"), ordered = TRUE))
df_sim_numeric <- df_sim %>% mutate(awareness_rank = awareness_rank)
model_numeric <- glm(follow_diet ~ gender + awareness_rank + exercise_freq, data = df_sim_numeric, family = binomial)
summary(model_numeric)
## 
## Call:
## glm(formula = follow_diet ~ gender + awareness_rank + exercise_freq, 
##     family = binomial, data = df_sim_numeric)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     0.07800    0.43585   0.179    0.858    
## genderMale     -1.03084    0.26272  -3.924 8.72e-05 ***
## awareness_rank  0.66147    0.10930   6.052 1.43e-09 ***
## exercise_freq   0.03220    0.08255   0.390    0.697    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 474.41  on 499  degrees of freedom
## Residual deviance: 415.75  on 496  degrees of freedom
## AIC: 423.75
## 
## Number of Fisher Scoring iterations: 5

10.3.2.1 Interpretasi Koefisien

  • awareness_rank: Setiap peningkatan satu tingkat dalam kesadaran gizi meningkatkan log-odds mengikuti diet sehat. Diasumsikan pengaruh antar level bersifat linier.
  • exercise_freq: Koefisien positif mengindikasikan olahraga lebih sering berkorelasi dengan kemungkinan mengikuti diet yang lebih tinggi.
  • gender: Sama seperti sebelumnya, digunakan untuk membandingkan laki-laki terhadap perempuan.

10.3.2.2 Interpretasi Goodness-of-Fit

  • Null deviance vs Residual deviance: Penurunan deviance menunjukkan prediktor menambah nilai pada model.
  • AIC: Bandingkan dengan model dummy; jika lebih kecil, model ordinal-numeric lebih efisien.

10.3.2.3 Signifikansi Model

  • Variabel awareness_rank, exercise_freq, dan gender diuji signifikansinya melalui p-value.
  • Jika semua signifikan, maka pendekatan linier terhadap kesadaran gizi valid.

10.3.2.4 Kesimpulan Praktis

  • Perlakuan ordinal sebagai numeric rank menyederhanakan model dan mengurangi jumlah parameter.
  • Cocok jika asumsi hubungan linier antar level kesadaran gizi dianggap wajar.

10.3.2.5 Visualisasi Prediksi Probabilitas (Ordinal sebagai Nominal)

df_sim_nominal <- df_sim_nominal %>% mutate(predicted = predict(model_nominal, type = "response"))

ggplot(df_sim_nominal, aes(x = exercise_freq, y = predicted, color = nutrition_awareness)) +
  geom_point(alpha = 0.5) +
  labs(title = "Prediksi Probabilitas Mengikuti Diet (Ordinal sebagai Nominal)",
       x = "Frekuensi Olahraga per Minggu", y = "Probabilitas Mengikuti Diet") +
  theme_minimal()

- Warna mewakili kategori tingkat kesadaran gizi (nutrition_awareness) yang diperlakukan sebagai variabel dummy. - Pola horizontal yang berbeda menunjukkan bahwa model memperlakukan tiap kategori sebagai entitas unik. - Kategori “Very High” memiliki probabilitas mengikuti diet paling tinggi, sedangkan “Very Low” paling rendah. - Frekuensi olahraga mingguan berkontribusi secara positif, namun pengaruhnya lebih moderat karena diasumsikan linier.

Kesimpulan: Model dummy cocok bila kita ingin menangkap variasi antar kategori tanpa asumsi linier.

10.3.2.6 Visualisasi Prediksi Probabilitas (Ordinal sebagai Numeric)

df_sim_numeric <- df_sim_numeric %>% mutate(predicted = predict(model_numeric, type = "response"))

ggplot(df_sim_numeric, aes(x = exercise_freq, y = predicted, color = as.factor(awareness_rank))) +
  geom_point(alpha = 0.5) +
  labs(title = "Prediksi Probabilitas Mengikuti Diet (Ordinal sebagai Numeric)",
       x = "Frekuensi Olahraga per Minggu", y = "Probabilitas Mengikuti Diet") +
  theme_minimal()

- Warna mewakili skor numerik kesadaran gizi (1 hingga 5). - Probabilitas meningkat secara lebih halus dan konsisten sesuai kenaikan level. - Model ini lebih sederhana karena hanya menggunakan satu koefisien untuk seluruh tingkatan ordinal.

Kesimpulan: Cocok jika kita percaya bahwa pengaruh antar level ordinal relatif linier.

Ringkasan Koefisien Model

library(broom)
library(kableExtra)

tidy(model_nominal) %>%
  kable(format = "latex", booktabs = TRUE, caption = "Koefisien Model Dummy") %>%
  kable_styling(latex_options = c("hold_position", "striped"))
tidy(model_numeric) %>%
  kable(format = "latex", booktabs = TRUE, caption = "Koefisien Model Ordinal Numeric") %>%
  kable_styling(latex_options = c("hold_position", "striped"))

Interpretasi: Tabel koefisien menunjukkan estimasi parameter model, standar error, nilai statistik z, dan p-value. Koefisien positif meningkatkan log-odds mengikuti diet sehat, sedangkan negatif menurunkan.

Kesimpulan: - Pendekatan dummy memperlihatkan bagaimana masing-masing kategori kesadaran gizi berbeda pengaruhnya terhadap peluang mengikuti diet. - Pendekatan ordinal numeric menyajikan versi sederhana dengan asumsi efek linier antar level. - Variabel olahraga dan gender juga penting untuk dipertimbangkan sebagai prediktor. - Visualisasi memperkuat interpretasi bahwa tingkat kesadaran dan kebiasaan olahraga berkorelasi positif dengan perilaku diet sehat.

11 Pemilihan Model Regresi Logistik dan Evaluasi

11.1 Membangun Model Regresi Logistik: Pendekatan Confirmatory dan Explorator

Regresi logistik digunakan ketika variabel respons bersifat kategorik (biasanya dikotomis/biner). Model ini memprediksi probabilitas kejadian berdasarkan variabel prediktor.

Ada dua pendekatan utama dalam membangun model regresi logistik:

  • Confirmatory (konfirmatori): Pendekatan berbasis teori.
  • Exploratory (eksploratori): Pendekatan berbasis data.
  1. Pendekatan Confirmatory

Pendekatan confirmatory mengacu pada penggunaan teori, penelitian sebelumnya, atau pemahaman domain untuk menentukan variabel prediktor yang masuk ke dalam model.

Karakteristik:

  • Teori atau hipotesis digunakan untuk menentukan variabel prediktor.
  • Validitas model diuji terhadap data.
  • Model yang dibangun tidak berubah meskipun hasil tidak signifikan.
  1. Pendekatan Exploratory

Pendekatan eksploratori membangun model berdasarkan pola yang ditemukan dalam data, tanpa hipotesis awal yang kaku.

Karakteristik:

  • Data-driven.
  • Proses seleksi variabel melalui metode statistik seperti:
    • Stepwise regression
    • Lasso regression
    • AIC/BIC selection
  • Cocok digunakan saat teori belum mapan.
  1. Perbandingan
Aspek Confirmatory Exploratory
Basis Teori Data
Fleksibilitas Rendah Tinggi
Risiko Overfitting Rendah Lebih tinggi
Validitas Teoritis Statistik
Tujuan Uji hipotesis Temukan pola
  1. Kesimpulan

Pendekatan confirmatory dan exploratory memiliki kelebihan masing-masing. Idealnya, kedua pendekatan dapat digunakan secara komplementer: model eksploratori bisa menjadi dasar untuk mengembangkan hipotesis baru yang kemudian diuji dalam model konfirmatori di masa depan.

Studi Kasus: Prediksi Kemungkinan Mengadopsi Energi Surya

Kita akan menggunakan contoh kasus berbeda: prediksi adopsi panel surya (adopt_solar = 1 jika mengadopsi, 0 jika tidak) berdasarkan: - income_level: Tingkat penghasilan (rendah, menengah, tinggi) - ordinal - environmental_concern: Tingkat kepedulian lingkungan (skala 1-5) - ordinal/numeric - home_ownership: Status kepemilikan rumah (Sewa vs. Milik) - nominal

set.seed(100)
n <- 400

income_level <- sample(c("Low", "Medium", "High"), n, replace = TRUE, prob = c(0.3, 0.5, 0.2))
environmental_concern <- sample(1:5, n, replace = TRUE)
home_ownership <- sample(c("Rent", "Own"), n, replace = TRUE, prob = c(0.4, 0.6))

logit_p <- -2 + 0.4 * (home_ownership == "Own") + 0.5 * as.numeric(factor(income_level, levels = c("Low", "Medium", "High"))) + 0.6 * environmental_concern
p <- 1 / (1 + exp(-logit_p))
adopt_solar <- rbinom(n, 1, p)

df_solar <- data.frame(adopt_solar, income_level, environmental_concern, home_ownership)
model_confirm <- glm(adopt_solar ~ income_level + environmental_concern + home_ownership, data = df_solar, family = binomial)
summary(model_confirm)
## 
## Call:
## glm(formula = adopt_solar ~ income_level + environmental_concern + 
##     home_ownership, family = binomial, data = df_solar)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -0.51682    0.34981  -1.477   0.1396    
## income_levelLow       -0.85438    0.34541  -2.474   0.0134 *  
## income_levelMedium    -0.41437    0.32441  -1.277   0.2015    
## environmental_concern  0.62728    0.08666   7.238 4.54e-13 ***
## home_ownershipRent    -0.34674    0.23242  -1.492   0.1357    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 521.57  on 399  degrees of freedom
## Residual deviance: 454.68  on 395  degrees of freedom
## AIC: 464.68
## 
## Number of Fisher Scoring iterations: 4

11.2 Metode Stepwise: Forward, Backward, dan Kedua Arah

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
model_null <- glm(adopt_solar ~ 1, data = df_solar, family = binomial)
model_full <- glm(adopt_solar ~ income_level + environmental_concern + home_ownership, data = df_solar, family = binomial)

step_forward <- stepAIC(model_null, scope = list(lower = model_null, upper = model_full), direction = "forward")
## Start:  AIC=523.57
## adopt_solar ~ 1
## 
##                         Df Deviance    AIC
## + environmental_concern  1   463.70 467.70
## <none>                       521.57 523.57
## + income_level           2   517.68 523.68
## + home_ownership         1   520.88 524.88
## 
## Step:  AIC=467.7
## adopt_solar ~ environmental_concern
## 
##                  Df Deviance    AIC
## + income_level    2   456.92 464.92
## + home_ownership  1   461.44 467.44
## <none>                463.70 467.70
## 
## Step:  AIC=464.92
## adopt_solar ~ environmental_concern + income_level
## 
##                  Df Deviance    AIC
## + home_ownership  1   454.68 464.68
## <none>                456.92 464.92
## 
## Step:  AIC=464.68
## adopt_solar ~ environmental_concern + income_level + home_ownership
step_backward <- stepAIC(model_full, direction = "backward")
## Start:  AIC=464.68
## adopt_solar ~ income_level + environmental_concern + home_ownership
## 
##                         Df Deviance    AIC
## <none>                       454.68 464.68
## - home_ownership         1   456.92 464.92
## - income_level           2   461.44 467.44
## - environmental_concern  1   517.06 525.06
step_both <- stepAIC(model_null, scope = list(lower = model_null, upper = model_full), direction = "both")
## Start:  AIC=523.57
## adopt_solar ~ 1
## 
##                         Df Deviance    AIC
## + environmental_concern  1   463.70 467.70
## <none>                       521.57 523.57
## + income_level           2   517.68 523.68
## + home_ownership         1   520.88 524.88
## 
## Step:  AIC=467.7
## adopt_solar ~ environmental_concern
## 
##                         Df Deviance    AIC
## + income_level           2   456.92 464.92
## + home_ownership         1   461.44 467.44
## <none>                       463.70 467.70
## - environmental_concern  1   521.57 523.57
## 
## Step:  AIC=464.92
## adopt_solar ~ environmental_concern + income_level
## 
##                         Df Deviance    AIC
## + home_ownership         1   454.68 464.68
## <none>                       456.92 464.92
## - income_level           2   463.70 467.70
## - environmental_concern  1   517.68 523.68
## 
## Step:  AIC=464.68
## adopt_solar ~ environmental_concern + income_level + home_ownership
## 
##                         Df Deviance    AIC
## <none>                       454.68 464.68
## - home_ownership         1   456.92 464.92
## - income_level           2   461.44 467.44
## - environmental_concern  1   517.06 525.06

11.3 Evaluasi Model: ROC dan AUC

library(pROC)
df_solar$pred <- predict(model_confirm, type = "response")
roc_curve <- roc(df_solar$adopt_solar, df_solar$pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_curve, col = "red")

auc(roc_curve)
## Area under the curve: 0.7398

11.4 Pseudo R-Squared

library(pscl)
pR2(model_confirm)
## fitting null model for pseudo-r2
##          llh      llhNull           G2     McFadden         r2ML         r2CU 
## -227.3415008 -260.7864829   66.8899642    0.1282466    0.1539907    0.2113695

11.5 Tabel Klasifikasi dan Evaluasi

thresh <- 0.5
df_solar$pred_class <- ifelse(df_solar$pred >= thresh, 1, 0)
table(Predicted = df_solar$pred_class, Actual = df_solar$adopt_solar)
##          Actual
## Predicted   0   1
##         0  70  41
##         1  73 216

11.6 Metode Perbandingan Model dalam Regresi Logistik

Perbandingan Model

library(dplyr)
library(knitr)
library(broom)

# Model bertingkat
model1 <- glm(adopt_solar ~ home_ownership, data = df_solar, family = binomial)
model2 <- glm(adopt_solar ~ home_ownership + income_level, data = df_solar, family = binomial)
model3 <- glm(adopt_solar ~ home_ownership + income_level + environmental_concern, data = df_solar, family = binomial)

Perbandingan AIC & Deviance

model_comp <- data.frame(
  Model = c("Model 1", "Model 2", "Model 3"),
  AIC = c(AIC(model1), AIC(model2), AIC(model3)),
  Deviance = c(deviance(model1), deviance(model2), deviance(model3))
)

knitr::kable(model_comp, booktabs = TRUE, caption = "Perbandingan Model Berdasarkan AIC dan Deviance")
Perbandingan Model Berdasarkan AIC dan Deviance
Model AIC Deviance
Model 1 524.8809 520.8809
Model 2 525.0634 517.0634
Model 3 464.6830 454.6830

11.7 Likelihood-Ratio Test

anova(model1, model2, test = "LRT")
## Analysis of Deviance Table
## 
## Model 1: adopt_solar ~ home_ownership
## Model 2: adopt_solar ~ home_ownership + income_level
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       398     520.88                     
## 2       396     517.06  2   3.8175   0.1483
anova(model2, model3, test = "LRT")
## Analysis of Deviance Table
## 
## Model 1: adopt_solar ~ home_ownership + income_level
## Model 2: adopt_solar ~ home_ownership + income_level + environmental_concern
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       396     517.06                          
## 2       395     454.68  1    62.38 2.831e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

11.8 Prinsip Parsimony

11.9 Evaluasi Tabel Klasifikasi dan Akurasi Model

library(caret)

# Prediksi probabilitas
pred_prob <- predict(model3, type = "response")

# Klasifikasi berdasarkan threshold 0.5
pred_class <- factor(ifelse(pred_prob >= 0.5, 1, 0), levels = c(0, 1))

# Konversi response asli menjadi faktor juga
actual_class <- factor(df_solar$adopt_solar, levels = c(0, 1))

# Matriks konfusi
conf_matrix <- confusionMatrix(pred_class, actual_class, positive = "1")
conf_matrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0  70  41
##          1  73 216
##                                          
##                Accuracy : 0.715          
##                  95% CI : (0.668, 0.7588)
##     No Information Rate : 0.6425         
##     P-Value [Acc > NIR] : 0.001267       
##                                          
##                   Kappa : 0.3472         
##                                          
##  Mcnemar's Test P-Value : 0.003691       
##                                          
##             Sensitivity : 0.8405         
##             Specificity : 0.4895         
##          Pos Pred Value : 0.7474         
##          Neg Pred Value : 0.6306         
##              Prevalence : 0.6425         
##          Detection Rate : 0.5400         
##    Detection Prevalence : 0.7225         
##       Balanced Accuracy : 0.6650         
##                                          
##        'Positive' Class : 1              
##

11.9.1 Sensitivitas dan Spesifisitas

  • Sensitivitas: Kemampuan model mendeteksi kelas positif secara benar (True Positive Rate)
    \(\text{Sensitivity} = \frac{TP}{TP + FN}\)
  • Spesifisitas: Kemampuan model mendeteksi kelas negatif secara benar (True Negative Rate)
    \(\text{Specificity} = \frac{TN}{TN + FP}\)
conf_matrix$byClass[c("Sensitivity", "Specificity")]
## Sensitivity Specificity 
##   0.8404669   0.4895105

11.10 Detail ROC, Penjelasan Kurva ROC (Receiver Operating Characteristic)

11.10.1 Definisi

ROC (Receiver Operating Characteristic) curve adalah grafik yang menunjukkan performa model klasifikasi biner pada berbagai nilai threshold. Kurva ini membandingkan True Positive Rate (TPR) terhadap False Positive Rate (FPR).

  • TPR (Sensitivitas): Proporsi kasus positif yang terdeteksi benar
  • FPR: Proporsi kasus negatif yang salah diklasifikasikan sebagai positif

11.10.2 Cut-Off dan Pergerakan Kurva

ROC curve dibuat dengan mengubah nilai cut-off dari 0 hingga 1. Setiap perubahan threshold menghasilkan pasangan nilai TPR dan FPR baru.

  • Cut-off rendah → lebih banyak prediksi positif → TPR naik, FPR juga naik
  • Cut-off tinggi → lebih sedikit prediksi positif → TPR turun, FPR turun

11.10.3 Kurva ROC Ideal

Kurva ROC ideal: - Dimulai dari titik (0,0), melewati (0,1), dan berakhir di (1,1) - Mewakili model sempurna yang dapat mengklasifikasikan semua kasus dengan benar - Semakin dekat kurva ke sudut kiri atas, semakin baik performa model

11.10.4 Interpretasi Luas Area (AUC)

AUC (Area Under the Curve) mengukur kualitas model klasifikasi: - AUC = 1.0 → model sempurna - AUC = 0.5 → model tidak lebih baik dari tebak-tebakan - AUC antara 0.7-0.9 → model baik

11.10.5 Kegunaan Kurva ROC

  • Mengevaluasi performa model klasifikasi
  • Membandingkan beberapa model
  • Memilih threshold optimal
  • Cocok untuk data tidak seimbang (imbalanced)

11.10.6 Visualisasi dalam R

library(pROC)

# Prediksi dan ROC
pred <- predict(model3, type = "response")
roc_obj <- roc(df_solar$adopt_solar, pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
# Plot aman untuk RMarkdown (jangan pakai abline)
plot(roc_obj, col = "red", lwd = 2, main = "ROC Curve - Model 3")

# Tampilkan AUC
auc(roc_obj)
## Area under the curve: 0.7398

11.10.7 Simulasi Pemilihan Threshold Optimal

Simulasi ini bertujuan untuk menentukan ambang batas (threshold) terbaik untuk klasifikasi adopsi energi surya berdasarkan prediksi probabilitas dari model regresi logistik biner.

# Threshold dari 0.1 hingga 0.9 dengan interval 0.05
thresholds <- seq(0.1, 0.9, by = 0.05)
results <- data.frame(Threshold = thresholds)

# Prediksi probabilitas dari model3
pred <- predict(model3, type = "response")

# Hitung sensitivitas
results$Sensitivity <- sapply(thresholds, function(t) {
  pred_class <- ifelse(pred >= t, 1, 0)
  cm <- table(Pred = pred_class, Obs = df_solar$adopt_solar)
  TP <- ifelse("1" %in% rownames(cm) & "1" %in% colnames(cm), cm["1", "1"], 0)
  FN <- ifelse("0" %in% rownames(cm) & "1" %in% colnames(cm), cm["0", "1"], 0)
  if ((TP + FN) == 0) return(NA)
  TP / (TP + FN)
})

# Hitung spesifisitas
results$Specificity <- sapply(thresholds, function(t) {
  pred_class <- ifelse(pred >= t, 1, 0)
  cm <- table(Pred = pred_class, Obs = df_solar$adopt_solar)
  TN <- ifelse("0" %in% rownames(cm) & "0" %in% colnames(cm), cm["0", "0"], 0)
  FP <- ifelse("1" %in% rownames(cm) & "0" %in% colnames(cm), cm["1", "0"], 0)
  if ((TN + FP) == 0) return(NA)
  TN / (TN + FP)
})

# Tambahkan kolom total
results$Total <- results$Sensitivity + results$Specificity

# Tampilkan tabel hasil
results
##    Threshold Sensitivity Specificity    Total
## 1       0.10   1.0000000  0.00000000 1.000000
## 2       0.15   1.0000000  0.00000000 1.000000
## 3       0.20   1.0000000  0.00000000 1.000000
## 4       0.25   1.0000000  0.00000000 1.000000
## 5       0.30   0.9922179  0.04195804 1.034176
## 6       0.35   0.9494163  0.23076923 1.180186
## 7       0.40   0.9338521  0.26573427 1.199586
## 8       0.45   0.8793774  0.36363636 1.243014
## 9       0.50   0.8404669  0.48951049 1.329977
## 10      0.55   0.7859922  0.58041958 1.366412
## 11      0.60   0.6770428  0.66433566 1.341378
## 12      0.65   0.5992218  0.74825175 1.347474
## 13      0.70   0.5680934  0.79020979 1.358303
## 14      0.75   0.5136187  0.82517483 1.338794
## 15      0.80   0.3929961  0.88111888 1.274115
## 16      0.85   0.2762646  0.92307692 1.199342
## 17      0.90   0.1128405  0.97902098 1.091861
# Threshold optimal
optimal_threshold <- results$Threshold[which.max(results$Total)]
cat("Threshold optimal berdasarkan maksimum Sensitivitas + Spesifisitas:", optimal_threshold, "\n")
## Threshold optimal berdasarkan maksimum Sensitivitas + Spesifisitas: 0.55

11.10.8 Catatan

  • Threshold optimal dapat membantu meningkatkan akurasi klasifikasi.
  • Metode ini mempertimbangkan keseimbangan antara Sensitivity dan Specificity.
  • Untuk data tidak seimbang, Precision-Recall Curve bisa lebih informatif daripada ROC.

11.11 Precision-Recall Curve (PR Curve)

11.11.1 Definisi

Kurva Precision-Recall (PR Curve) adalah grafik yang menunjukkan hubungan antara Precision dan Recall untuk berbagai nilai threshold. Kurva ini sangat berguna untuk mengevaluasi performa model klasifikasi, khususnya ketika data tidak seimbang (kelas minoritas sangat sedikit).

  • Precision (Presisi): Rasio prediksi positif yang benar-benar positif
    \[ \text{Precision} = \frac{TP}{TP + FP} \]

  • Recall (Sensitivitas): Rasio kasus positif yang berhasil diprediksi positif
    \[ \text{Recall} = \frac{TP}{TP + FN} \]

11.11.2 Interpretasi

  • PR Curve menunjukkan bagaimana precision berubah seiring perubahan recall.
  • Idealnya, precision dan recall keduanya tinggi.
  • PR Curve yang baik akan melengkung ke pojok kanan atas.

11.11.3 Area Under PR Curve

  • AUPRC (Area Under Precision-Recall Curve) mendekati 1 berarti model sangat bagus.
  • Baseline AUPRC = prevalensi kelas positif dalam data.

11.11.4 PR Curve vs ROC Curve

Aspek ROC Curve Precision-Recall Curve
Fokus Semua kelas Kelas positif saja
Cocok untuk Data seimbang Data tidak seimbang
Sumbu X 1 - Specificity Recall
Sumbu Y Sensitivity (Recall) Precision

11.11.5 Visualisasi PR Curve di R

library(PRROC)
## Warning: package 'PRROC' was built under R version 4.3.3
## Loading required package: rlang
## 
## Attaching package: 'rlang'
## The following object is masked from 'package:magrittr':
## 
##     set_names
## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice
library(dplyr)

# Prediksi dan label aktual
scores <- predict(model3, type = "response")
labels <- df_solar$adopt_solar

# Buat PR curve
pr <- pr.curve(scores.class0 = scores[labels == 1],
               scores.class1 = scores[labels == 0],
               curve = TRUE)

# Plot
plot(pr, main = "Precision-Recall Curve - Model 3")

11.11.6 Kesimpulan

  • PR Curve lebih tepat digunakan dibanding ROC saat data tidak seimbang.
  • Memungkinkan identifikasi model dengan kombinasi precision dan recall yang optimal.

11.12 Pseudo R-Squared pada Regresi Logistik

Tujuan Bagian ini menjelaskan dan menghitung berbagai pseudo R-squared dalam regresi logistik untuk mengevaluasi performa model, termasuk:

  • R² Cox and Snell
  • R² McFadden

Likelihood dan Rumus Dengan: \[ R^2_{\text{Cox and Snell}} = 1 - \left(\frac{L_0}{L_M}\right)^{\frac{2}{n}} \quad \text{dan} \quad R^2_{\text{McFadden}} = 1 - \frac{\log L_M}{\log L_0} \]

  • \(L_0\): likelihood dari model null (tanpa prediktor)
  • \(L_M\): likelihood dari model penuh

Perhitungan Manual R²

model_full <- glm(adopt_solar ~ income_level + environmental_concern + home_ownership,
                  data = df_solar, family = binomial)

model_null <- glm(adopt_solar ~ 1, data = df_solar, family = binomial)

logL0 <- logLik(model_null)
logLM <- logLik(model_full)
L0 <- exp(logL0)
LM <- exp(logLM)
n <- nobs(model_full)

cox_snell <- 1 - (L0 / LM)^(2 / n)
mcfadden <- 1 - (as.numeric(logLM) / as.numeric(logL0))

r2_manual <- data.frame(R2_Cox_Snell = cox_snell, R2_McFadden = mcfadden)
r2_manual
##   R2_Cox_Snell R2_McFadden
## 1    0.1539907   0.1282466

DescTools

if (!require(DescTools)) install.packages("DescTools")
library(DescTools)
PseudoR2(model_full, which = "all")
##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.1282466       0.1090738       0.1539907       0.2113695       0.1432671 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.2531402       0.1617129       0.2135025       0.1604407     464.6830015 
##             BIC          logLik         logLik0              G2 
##     484.6403243    -227.3415008    -260.7864829      66.8899642

Interpretasi

  • R² McFadden > 0.2 menunjukkan model cukup baik dalam memprediksi kemungkinan seseorang mengadopsi panel surya berdasarkan faktor-faktor tersebut
  • R² Cox and Snell lebih konservatif dan tidak bisa mencapai nilai 1 penuh. Ini berarti informasi dari ketiga prediktor yang digunakan memang berkontribusi terhadap peningkatan akurasi model.
  • Nilai pseudo R-squared ini bukan interpretasi deterministik (tidak sama seperti R² di regresi linear), tapi tetap berguna sebagai ukuran komparatif untuk mengevaluasi dan membandingkan model klasifikasi logistik.

11.13 Kesimpulan Akhir

  • Model yang dibangun mampu memprediksi adopsi panel surya dengan cukup baik.
  • Pendekatan confirmatory dan eksploratory sama-sama memberi wawasan.
  • Stepwise bisa digunakan untuk menyederhanakan model.
  • Evaluasi dengan ROC, AUC, pseudo R², dan confusion matrix sangat penting.
  • Model sederhana namun informatif sebaiknya dipilih (parsimony).

12 Regresi Logistik Multinomial

Distribusi multinomial adalah perluasan dari distribusi binomial untuk lebih dari dua kategori. Jika dalam distribusi binomial kita hanya memiliki dua kemungkinan (misalnya sukses dan gagal), maka dalam distribusi multinomial kita bisa memiliki \(k\) kategori.

Misalnya, dalam suatu survei, responden dapat memilih salah satu dari tiga merek minuman ringan: A, B, atau C. Jika kita melakukan survei terhadap 100 orang, dan mencatat berapa orang yang memilih masing-masing merek, maka banyaknya orang yang memilih A, B, dan C akan mengikuti distribusi multinomial.

12.1 Notasi Matematis

Jika \(X_1, X_2, \dots, X_k\) menyatakan jumlah pengamatan dalam masing-masing dari \(k\) kategori, maka:

\[ P(X_1 = x_1, \dots, X_k = x_k) = \frac{n!}{x_1! x_2! \cdots x_k!} p_1^{x_1} p_2^{x_2} \cdots p_k^{x_k} \]

dengan syarat:

\[ \sum_{i=1}^{k} x_i = n, \quad \sum_{i=1}^{k} p_i = 1 \]

12.2 Studi Kasus : Pilihan Aktivitas Akhir Pekan

Distribusi multinomial digunakan ketika terdapat lebih dari dua kategori hasil yang saling eksklusif dalam suatu eksperimen dengan jumlah percobaan tetap. Berikut ini adalah sebuah studi kasus yang menggunakan distribusi ini.

Sebuah survei dilakukan terhadap 10 orang yang diminta memilih satu dari tiga jenis aktivitas akhir pekan favorit mereka: - Mendaki (M) - Menonton (N) - Berbelanja (B)

Hasil survei: - Mendaki: 2 orang - Menonton: 5 orang - Berbelanja: 3 orang

Probabilitas teoretik preferensi: - \(p_M = 0.2\) - \(p_N = 0.5\) - \(p_B = 0.3\)

Rumus Distribusi Multinomial \[ P(X_1 = x_1, \dots, X_k = x_k) = \frac{n!}{x_1! x_2! \cdots x_k!} p_1^{x_1} p_2^{x_2} \cdots p_k^{x_k} \]

dengan: - \(n = 10\) - \(x_1 = 2, x_2 = 5, x_3 = 3\) - \(p_1 = 0.2, p_2 = 0.5, p_3 = 0.3\)

Perhitungan Manual di R

n <- 10
x <- c(2, 5, 3)
p <- c(0.2, 0.5, 0.3)

# Hitung komponen-koefisien
faktorial_total <- factorial(n)
faktorial_x <- prod(factorial(x))
koefisien <- faktorial_total / faktorial_x

# Hitung peluang
peluang <- koefisien * prod(p^x)
peluang
## [1] 0.08505

Interpretasi : Probabilitas bahwa 2 orang memilih Mendaki, 5 Menonton, dan 3 Berbelanja (dengan proporsi preferensi masing-masing 0.2, 0.5, dan 0.3) adalah nilai dari peluang yang dihitung di atas. Distribusi multinomial digunakan untuk menghitung peluang dalam percobaan dengan beberapa kategori hasil. Rumus dasarnya merupakan generalisasi dari distribusi binomial untuk lebih dari dua kategori.

Distribusi ini dapat diperluas untuk kategori yang lebih banyak dan sampel yang lebih besar. Distribusi multinomial sangat umum digunakan dalam pengolahan data survei, polling, dan eksperimen yang hasilnya berupa kategori.

12.3 Multinomial Logistic Regression

Multinomial logistic regression digunakan untuk memodelkan hubungan antara satu variabel respon kategorik (dengan lebih dari dua kategori) dan satu atau lebih variabel prediktor.

Misalkan variabel respon \(Y\) memiliki \(K\) kategori, dan kita memilih kategori ke-\(K\) sebagai baseline, maka bentuk umum model logit untuk kategori \(j\) adalah: \[ \log\left(\frac{P(Y = j)}{P(Y = K)}\right) = \beta_{j0} + \beta_{j1}x_1 + \cdots + \beta_{jp}x_p, \quad \text{untuk } j = 1,2,\dots,K-1 \]

12.3.1 Baseline-category logit model

Baseline-category logit model adalah bentuk regresi logistik untuk variabel kategorik nominal dengan lebih dari dua kategori. Model ini menggunakan salah satu kategori sebagai baseline dan menghitung logit untuk membandingkan tiap kategori lain terhadap baseline tersebut:

\[ \log\left(\frac{\pi_j}{\pi_c}\right), \quad j = 1,\dots,c - 1 \]

dengan: - \(\pi_j\): probabilitas respon berada di kategori \(j\) - \(\pi_c\): probabilitas respon berada di kategori baseline \(c\)

Jumlah fungsi logit adalah sebanyak \(c - 1\).

Catatan: Default kategori baseline di R (fungsi multinom() dari paket nnet) adalah kategori terakhir, tapi bisa diubah dengan fungsi relevel().

Model Regresi Sederhana (1 Prediktor) Jika hanya terdapat satu prediktor \(x\), maka bentuk model logit: \[ \log\left(\frac{\pi_j}{\pi_c}\right) = \alpha_j + \beta_j x, \quad j = 1, \dots, c - 1 \]

Contoh Kasus (3 Kategori) Misalkan \(Y \in \{1, 2, 3\}\), dengan kategori ke-3 sebagai baseline. Maka: \[ \log\left(\frac{\pi_1}{\pi_3}\right) = \alpha_1 + \beta_1 x \] \[ \log\left(\frac{\pi_2}{\pi_3}\right) = \alpha_2 + \beta_2 x \]

Untuk membandingkan langsung antara kategori 1 dan 2: \[ \log\left(\frac{\pi_1}{\pi_2}\right) = \log\left(\frac{\pi_1/\pi_3}{\pi_2/\pi_3}\right) = (\alpha_1 - \alpha_2) + (\beta_1 - \beta_2)x \]

Karakteristik Model: - Digunakan untuk kategori nominal (> 2 kategori) - Menghasilkan \(c - 1\) fungsi logit terhadap baseline - Logit antar kategori selain baseline dapat dihitung dari selisih logit terhadap baseline

12.3.2 Estimasi Parameter

Estimasi parameter dilakukan menggunakan metode Maximum Likelihood Estimation (MLE), dengan algoritma iteratif seperti Newton-Raphson atau iteratively reweighted least squares (IRLS).

Log-likelihood umum untuk multinomial logit: \[ \ell(\beta) = \sum_{i=1}^n \sum_{j=1}^{K} y_{ij} \log(\pi_{ij}) \]

dengan: - \(\pi_{ij} = P(Y_i = j | x_i)\): probabilitas prediksi - \(y_{ij} = 1\) jika observasi ke-\(i\) berada pada kategori ke-\(j\), dan 0 jika tidak

12.4 Contoh Kasus

Sebuah lembaga transportasi ingin mengetahui preferensi warga kota terhadap moda transportasi utama yang digunakan: Bus, Kereta, atau Motor. Mereka melakukan survei terhadap 200 responden, dan mencatat beberapa informasi sebagai berikut:

  • Transport: Moda transportasi yang dipilih (Bus, Kereta, Motor)
  • Age: Usia responden
  • Area: Lokasi tempat tinggal (Urban, Suburban, Rural)
  • Income: Pendapatan bulanan (dalam juta rupiah)

Tujuan: Mengetahui bagaimana usia, lokasi tempat tinggal, dan pendapatan memengaruhi pilihan moda transportasi utama.

12.5 Simulasi Data

library(nnet)
library(dplyr)
set.seed(2025)
n <- 200

Area <- sample(c("Urban", "Suburban", "Rural"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 30, sd = 8))
Income <- round(rnorm(n, mean = 7, sd = 2), 1)

# Simulasi pilihan transportasi berdasarkan area
Transport <- sapply(Area, function(loc) {
  if (loc == "Urban") {
    sample(c("Bus", "Kereta", "Motor"), size = 1, prob = c(0.5, 0.4, 0.1))
  } else if (loc == "Suburban") {
    sample(c("Bus", "Kereta", "Motor"), size = 1, prob = c(0.3, 0.3, 0.4))
  } else {
    sample(c("Bus", "Kereta", "Motor"), size = 1, prob = c(0.2, 0.2, 0.6))
  }
})

transport_data <- data.frame(
  Transport = factor(Transport),
  Age,
  Area = factor(Area),
  Income
)

# Set baseline
transport_data$Transport <- relevel(transport_data$Transport, ref = "Bus")
head(transport_data)
##   Transport Age     Area Income
## 1       Bus  35    Urban    8.1
## 2    Kereta  31 Suburban    3.6
## 3    Kereta  25    Urban    3.2
## 4     Motor  40    Rural   11.8
## 5     Motor  30 Suburban    8.1
## 6     Motor  43    Rural    9.7

12.6 Estimasi Model

model_transport <- multinom(Transport ~ Age + Area + Income, data = transport_data)
## # weights:  18 (10 variable)
## initial  value 219.722458 
## iter  10 value 194.940253
## final  value 194.306803 
## converged
summary(model_transport)
## Call:
## multinom(formula = Transport ~ Age + Area + Income, data = transport_data)
## 
## Coefficients:
##        (Intercept)         Age AreaSuburban  AreaUrban     Income
## Kereta    1.602297 -0.05503436   -0.5493588 -0.4397256 0.03250565
## Motor     1.796763 -0.02001485   -1.3128581 -2.3475772 0.05535300
## 
## Std. Errors:
##        (Intercept)        Age AreaSuburban AreaUrban     Income
## Kereta   1.0995524 0.02596633    0.5386943 0.5276212 0.10030533
## Motor    0.9958387 0.02309125    0.4429547 0.4938124 0.09017309
## 
## Residual Deviance: 388.6136 
## AIC: 408.6136

12.7 Nilai P-Value dan Interpretasi

z <- summary(model_transport)$coefficients / summary(model_transport)$standard.errors
pval <- 2 * (1 - pnorm(abs(z)))
round(pval, 4)
##        (Intercept)    Age AreaSuburban AreaUrban Income
## Kereta      0.1451 0.0341       0.3078    0.4046 0.7459
## Motor       0.0712 0.3861       0.0030    0.0000 0.5393

Interpretasi: - Koefisien mengindikasikan pengaruh terhadap logit dibandingkan kategori referensi (Bus). - Misalnya, jika p-value untuk AreaRural signifikan terhadap kategori Motor, maka tinggal di daerah rural meningkatkan kemungkinan memilih motor dibanding bus.

12.8 Prediksi dan Validasi

transport_data$Predicted <- predict(model_transport)
table(Predicted = transport_data$Predicted, Actual = transport_data$Transport)
##          Actual
## Predicted Bus Kereta Motor
##    Bus     25     16    10
##    Kereta   4      5     5
##    Motor   31     28    76

Interpretasi: - Tabel prediksi vs aktual menunjukkan seberapa baik model mampu memprediksi kategori transportasi berdasarkan fitur input.

12.9 Kesimpulan

Model logistik multinomial berhasil digunakan untuk menganalisis pengaruh usia, pendapatan, dan lokasi tempat tinggal terhadap preferensi moda transportasi. - Area tempat tinggal terbukti sebagai variabel penting yang memengaruhi preferensi. - Model ini dapat digunakan untuk membantu perencanaan transportasi dan kebijakan berbasis wilayah.

12.9.1 Contoh Kasus 2 di R

Kita gunakan dataset mtcars yang telah disesuaikan untuk memodelkan jenis transmisi (am: 0 = automatic, 1 = manual, dan kita tambahkan kategori baru “semi-auto”) berdasarkan hp (horsepower) dan wt (weight).

library(nnet)
mtcars2 <- mtcars %>%
  mutate(am = as.character(am)) %>%
  mutate(am = ifelse(hp > 180, "semi-auto", ifelse(am == "0", "auto", "manual"))) %>%
  mutate(am = factor(am))

mtcars2$am <- relevel(mtcars2$am, ref = "auto")
model_case2 <- multinom(am ~ hp + wt, data = mtcars2)
## # weights:  12 (6 variable)
## initial  value 35.155593 
## iter  10 value 6.142895
## iter  20 value 5.405263
## iter  30 value 4.785759
## iter  40 value 4.696327
## iter  50 value 4.647685
## iter  60 value 4.627699
## iter  70 value 4.601272
## iter  80 value 4.582075
## iter  90 value 4.572830
## iter 100 value 4.559717
## final  value 4.559717 
## stopped after 100 iterations
summary(model_case2)
## Call:
## multinom(formula = am ~ hp + wt, data = mtcars2)
## 
## Coefficients:
##           (Intercept)         hp        wt
## manual       16.66255 0.03494577 -7.235147
## semi-auto   -47.15189 0.17383746  2.946640
## 
## Std. Errors:
##           (Intercept)        hp       wt
## manual       6.778766 0.0410300 3.100704
## semi-auto   36.878574 0.1419153 5.318602
## 
## Residual Deviance: 9.119435 
## AIC: 21.11943
z <- summary(model_case2)$coefficients / summary(model_case2)$standard.errors
p_values <- 2 * (1 - pnorm(abs(z)))
round(p_values, 4)
##           (Intercept)     hp     wt
## manual          0.014 0.3944 0.0196
## semi-auto       0.201 0.2206 0.5796
mtcars2$predicted <- predict(model_case2, newdata = mtcars2)
table(Predicted = mtcars2$predicted, Actual = mtcars2$am)
##            Actual
## Predicted   auto manual semi-auto
##   auto        13      1         0
##   manual       1     10         0
##   semi-auto    0      0         7
library(ggplot2)
ggplot(mtcars2, aes(x = hp, y = wt, color = predicted)) +
  geom_point(size = 3) +
  labs(title = "Multinomial Logistic Regression Predictions",
       x = "Horsepower", y = "Weight") +
  theme_minimal()

Kesimpulan:

  • Prediksi model terhadap jenis transmisi cukup akurat dalam kasus ini.
  • Multinomial logistic regression efektif digunakan untuk klasifikasi dengan >2 kategori.

13 Regresi Logistik Ordinal

Digunakan ketika variabel respon 𝑌 bersifat ordinal (memiliki urutan), misalnya tingkat kepuasan: Rendah, Sedang, Tinggi. Model ini berbeda dengan: - Regresi logistik biner yang hanya memiliki 2 kategori
- Regresi logistik multinomial yang memiliki kategori > 2 tetapi tidak berurutan

13.1 Konsep Cumulative Logit Model

Model regresi logistik ordinal digunakan ketika variabel respons memiliki urutan kategori. Salah satu pendekatan paling umum adalah Cumulative Logit Model, juga dikenal sebagai proportional odds model. Bentuk umum model ini adalah:

\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j + \beta x \]

Dengan: - \(\alpha_j\) adalah intercept untuk batas kumulatif ke-\(j\) - \(\beta\) adalah koefisien regresi yang sama untuk semua batas (asumsi proportional odds)

Jika ada \(c\) kategori ordinal, maka akan dibentuk sebanyak \(c - 1\) model logit kumulatif.

13.2 Interpretasi Koefisien

Koefisien \(\beta\) menginterpretasikan pengaruh dari variabel prediktor terhadap kemungkinan respons berada dalam kategori rendah atau sama (\(Y \leq j\)):

  • Jika \(\beta > 0\): Semakin besar nilai \(x\), semakin tinggi peluang berada dalam kategori lebih rendah.
  • Jika \(\beta < 0\): Semakin besar nilai \(x\), semakin tinggi peluang berada dalam kategori lebih tinggi.

13.2.1 Odds Ratio

Odds ratio dihitung sebagai:

\[ \text{OR} = e^{\beta} \]

  • \(\text{OR} > 1\): variabel prediktor meningkatkan peluang berada di kategori lebih rendah.
  • \(\text{OR} < 1\): variabel prediktor meningkatkan peluang berada di kategori lebih tinggi.

Model ini sangat cocok untuk data dengan respon ordinal, seperti tingkat kepuasan (rendah - sedang - tinggi), tingkat pendidikan, atau frekuensi kegiatan.

13.3 Contoh Data:Tingkat Kesulitan Mata Kuliah

Misalkan kita memiliki data fiktif tentang tingkat kesulitan mata kuliah yang dirasakan mahasiswa berdasarkan jumlah jam belajar mandiri per minggu. Skala kesulitan terdiri dari:

  1. Mudah
  2. Sedang
  3. Sulit
library(MASS)
set.seed(2025)
n <- 200
study_hours <- round(runif(n, 1, 10))
difficulty_level <- cut(5 - 0.4 * study_hours + rnorm(n),
                        breaks = c(-Inf, 3.5, 5.5, Inf),
                        labels = c("Mudah", "Sedang", "Sulit"),
                        ordered_result = TRUE)
df <- data.frame(difficulty_level, study_hours)
head(df)
##   difficulty_level study_hours
## 1            Mudah           8
## 2            Mudah           5
## 3           Sedang           6
## 4            Mudah           5
## 5            Mudah           8
## 6            Mudah           6

13.4 Estimasi Model Ordinal

model_ord <- polr(difficulty_level ~ study_hours, data = df, Hess = TRUE)
summary(model_ord)
## Call:
## polr(formula = difficulty_level ~ study_hours, data = df, Hess = TRUE)
## 
## Coefficients:
##               Value Std. Error t value
## study_hours -0.7282    0.09946  -7.322
## 
## Intercepts:
##              Value   Std. Error t value
## Mudah|Sedang -2.8454  0.4709    -6.0419
## Sedang|Sulit  1.1894  0.5844     2.0351
## 
## Residual Deviance: 187.2484 
## AIC: 193.2484

13.5 Nilai P-Value

(ctable <- coef(summary(model_ord)))
##                   Value Std. Error   t value
## study_hours  -0.7281874 0.09945644 -7.321671
## Mudah|Sedang -2.8453574 0.47094103 -6.041855
## Sedang|Sulit  1.1893617 0.58441951  2.035116
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = round(p, 4)))
##                   Value Std. Error   t value p value
## study_hours  -0.7281874 0.09945644 -7.321671  0.0000
## Mudah|Sedang -2.8453574 0.47094103 -6.041855  0.0000
## Sedang|Sulit  1.1893617 0.58441951  2.035116  0.0418

13.6 Prediksi Probabilitas

newdata <- data.frame(study_hours = 4:8)
predict(model_ord, newdata = newdata, type = "probs")
##       Mudah     Sedang        Sulit
## 1 0.5168416 0.46688969 0.0162686790
## 2 0.6890281 0.30305103 0.0079208966
## 3 0.8210925 0.17506771 0.0038398077
## 4 0.9048190 0.09332356 0.0018574840
## 5 0.9516689 0.04743350 0.0008976247

Interpretasi: - Probabilitas meningkat pada kategori “Mudah” seiring dengan peningkatan jumlah jam belajar. - Model logistik ordinal membantu mengukur pengaruh variabel numerik terhadap kategori bertingkat seperti persepsi kesulitan mata kuliah.

13.7 Goodness-of-Fit dan Proportional Odds

Model cumulative logit mengasumsikan efek prediktor sama untuk setiap cutoff. Jika tidak, per timbangkan model non-proportional odds seperti generalized ordinal model.

13.8 Alternatif Model Ordinal

Jika asumsi proportional odds tidak cocok dengan data, ada beberapa pendekatan lain:

  • Adjacent-category logit model: membandingkan peluang berada pada satu kategori dibanding yang berdekatan.
  • Continuation-ratio logit model: cocok untuk data bertahap atau urutan proses.

Model-model ini fleksibel dan dapat menangani situasi di mana efek prediktor berubah antar tingkat kategori.

13.9 Kesimpulan

  • Regresi logistik ordinal sangat berguna untuk data kategorik berurutan.
  • Model cumulative logit adalah pendekatan populer yang menghasilkan interpretasi log-odds kumulatif.
  • Implementasi umum di R menggunakan fungsi polr() dari paket MASS.
  • Validasi model dapat dilakukan melalui uji deviance atau likelihood ratio test untuk membandingkan kecocokan model.

13.10 Asumsi Paralelisme dalam Regresi Logistik Ordinal

Asumsi paralelisme (parallel lines assumption) menyatakan bahwa semua fungsi logit kumulatif memiliki koefisien regresi yang sama. Dalam bentuk matematis:

\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j + \beta x \]

  • Hanya intercept (\(\alpha_j\)) yang bervariasi untuk tiap cutoff kategori.
  • Koefisien \(\beta\) tetap konstan di semua model logit kumulatif.

13.10.1 Visualisasi

Dalam asumsi ini, garis logit dari setiap kategori akan memiliki kemiringan sejajar, meskipun posisinya berbeda.

13.10.2 Jika Asumsi Dilanggar

  • Efek variabel prediktor tidak konsisten di seluruh batas kategori.
  • Model cumulative logit menjadi tidak tepat.
  • Solusi: gunakan generalized ordinal logistic regression atau partial proportional odds model.

13.10.3 Pengujian Asumsi Paralelisme

Asumsi paralelisme dapat diuji melalui:

  • Likelihood Ratio Test antara model proportional dan non-proportional.
  • Brant Test menggunakan package brant.

13.10.4 Contoh Pengujian dengan Brant Test

Misalkan kita memiliki dataset tentang tingkat urgensi laporan kebakaran yang dikategorikan sebagai: “Rendah”, “Sedang”, “Tinggi”.

library(MASS)
library(brant)
## Warning: package 'brant' was built under R version 4.3.3
set.seed(2025)
n <- 150
respon_time <- round(runif(n, 1, 10))
fire_urgency <- cut(3 + 0.5 * respon_time + rnorm(n),
                    breaks = c(-Inf, 5.5, 7.5, Inf),
                    labels = c("Rendah", "Sedang", "Tinggi"),
                    ordered_result = TRUE)
dat <- data.frame(fire_urgency, respon_time)
model <- polr(fire_urgency ~ respon_time, data = dat, Hess = TRUE)
brant(model)
## -------------------------------------------- 
## Test for X2  df  probability 
## -------------------------------------------- 
## Omnibus      0.01    1   0.91
## respon_time  0.01    1   0.91
## -------------------------------------------- 
## 
## H0: Parallel Regression Assumption holds

Jika p-value dari hasil brant() kurang dari 0.05, maka asumsi tidak terpenuhi.

Kesimpulan: - Asumsi paralelisme penting untuk validitas model cumulative logit. - Menyederhanakan interpretasi dengan efek prediktor yang tetap. - Jika tidak terpenuhi, gunakan model alternatif yang lebih fleksibel.

14 Log Linear Model

Dalam dunia nyata, banyak data yang dikumpulkan berupa kategori, seperti jenis kelamin, tingkat pendidikan, status pekerjaan, atau hasil preferensi pelanggan. Untuk menganalisis jenis data ini, terdapat tiga pendekatan utama yang biasa digunakan, yaitu: tabel kontingensi, model log-linear, dan regresi logistik.

Tabel kontingensi digunakan untuk menyajikan jumlah observasi yang berada dalam kombinasi kategori tertentu dari dua atau lebih variabel. Ini sering menjadi langkah awal untuk mengeksplorasi hubungan antar variabel. Sebagai contoh, kita bisa melihat distribusi pasien yang mengalami atau tidak mengalami serangan jantung berdasarkan jenis obat yang dikonsumsi. Dari tabel ini, kita dapat menghitung ukuran asosiasi seperti odds ratio atau melakukan uji chi-square untuk menguji apakah ada hubungan antar variabel.

Namun, ketika kita ingin memahami hubungan yang lebih kompleks antara beberapa variabel sekaligus, model log-linear menjadi pilihan yang tepat. Model ini merupakan bagian dari keluarga Generalized Linear Model (GLM) dan biasanya digunakan untuk memodelkan frekuensi sel dalam tabel kontingensi. Asumsinya adalah bahwa data mengikuti distribusi Poisson. Yang membedakan model log-linear dengan regresi logistik adalah bahwa model log-linear memperlakukan semua variabel secara setara; tidak ada variabel yang secara eksplisit dianggap sebagai respon. Oleh karena itu, model ini lebih cocok digunakan ketika tujuan analisis adalah mengeksplorasi asosiasi dan independensi antar variabel.

Struktur dari model log-linear ditentukan berdasarkan efek utama dari masing-masing variabel dan interaksi antar variabel. Misalnya, dalam sebuah tabel tiga arah yang memuat jenis kelamin, status merokok, dan kejadian penyakit paru, kita bisa menguji apakah interaksi antara dua variabel sudah cukup menjelaskan data atau perlu menambahkan interaksi tiga arah. Perbandingan model yang lebih sederhana dan lebih kompleks dilakukan menggunakan likelihood ratio test.

Sebaliknya, jika kita memiliki satu variabel kategorik sebagai respon utama (seperti: apakah seseorang sakit atau tidak), dan beberapa variabel sebagai prediktor, maka pendekatan regresi logistik menjadi lebih sesuai. Model ini memodelkan peluang terjadinya suatu outcome dalam bentuk logit (log odds) dan sangat sering digunakan dalam studi observasional maupun eksperimental. Regresi logistik juga bisa diperluas untuk menangani lebih dari dua kategori, seperti pada regresi logistik multinomial dan ordinal.

Kesimpulannya, meskipun ketiganya digunakan pada data kategorik, pendekatan dan tujuannya sangat berbeda:

Pemilihan pendekatan yang tepat tergantung pada fokus analisis: apakah deskriptif, eksploratif, atau prediktif. Dalam praktiknya, ketiga metode ini sering digunakan secara komplementer untuk memperoleh pemahaman yang lebih menyeluruh terhadap struktur data kategorik.

Tabel Kontingensi

table_data <- matrix(c(40, 30, 25, 55), nrow=2,
                     dimnames = list(Status = c("Pekerja", "Pengangguran"),
                                     Pendidikan = c("Tinggi", "Rendah")))
table_data
##               Pendidikan
## Status         Tinggi Rendah
##   Pekerja          40     25
##   Pengangguran     30     55

Tabel kontingensi di atas menunjukkan distribusi frekuensi berdasarkan status pekerjaan dan tingkat pendidikan. Analisis tabel ini bersifat deskriptif dan tidak melibatkan model probabilistik.

Model Loglinear

Model loglinear memodelkan logaritma dari ekspektasi frekuensi sel dalam tabel kontingensi.

\[ \log(\mu_{ij}) = \lambda + \lambda^A_i + \lambda^B_j + \lambda^{AB}_{ij} \]

library(MASS)
loglm(~ Status * Pendidikan, data = table_data)
## Call:
## loglm(formula = ~Status * Pendidikan, data = table_data)
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1

Model Regresi Logistik Model regresi logistik biner:

\[ \log\left( \frac{p_1}{1 - p} \right) = \beta_0 + \beta_1 x \]

Contoh:

data_glm <- data.frame(
  Bekerja = c(1, 0, 1, 0),
  Pendidikan = factor(c("Tinggi", "Tinggi", "Rendah", "Rendah")),
  Frekuensi = c(40, 30, 25, 55)
)

model_logit <- glm(Bekerja ~ Pendidikan, weights = Frekuensi, family = binomial, data = data_glm)
summary(model_logit)
## 
## Call:
## glm(formula = Bekerja ~ Pendidikan, family = binomial, data = data_glm, 
##     weights = Frekuensi)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)       -0.7885     0.2412  -3.269  0.00108 **
## PendidikanTinggi   1.0761     0.3413   3.153  0.00162 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 205.27  on 3  degrees of freedom
## Residual deviance: 194.98  on 2  degrees of freedom
## AIC: 198.98
## 
## Number of Fisher Scoring iterations: 4

Perbandingan Ketiga Pendekatan

Aspek Tabel Kontingensi Model Loglinear Regresi Logistik
Tujuan Deskripsi frekuensi Deteksi asosiasi Prediksi probabilitas
Variabel Dependen Tidak ada Tidak ada (simetris) Ada (eksplisit)
Distribusi Tidak diasumsikan Poisson (frekuensi sel) Binomial (probabilitas)
Bentuk Model Tidak ada GLM: log(mu) ~ efek GLM: logit(p) ~ prediktor
Cocok untuk Eksplorasi awal Tabel > 2 variabel Studi prediktif

14.1 Model Log-Linear Dua Arah

Model log-linear digunakan untuk menganalisis frekuensi dari tabel kontingensi dengan dua atau lebih variabel kategorik. Untuk dua variabel kategorik, misalnya A dan B, model log-linear memiliki dua bentuk ekstrem yang penting untuk dikenali:

14.1.0.1 Model Saturated

Model saturated adalah model log-linier yang mencakup semua efek utama dan interaksi:

  • Cocok sempurna terhadap data
  • Tidak mengasumsikan independensi antar variabel
  • Tidak ada residual (deviasi) antara nilai observasi dan ekspektasi

Rumus umum untuk model saturated:

\[ \log(\mu_{ij}) = \lambda + \lambda^A_i + \lambda^B_j + \lambda^{AB}_{ij} \]

Model ini cocok jika ingin memodelkan semua variasi dalam data, termasuk interaksi antar variabel.

14.1.0.2 Model Independen

Model independen mengasumsikan tidak ada interaksi antara variabel (hanya efek utama):

\[ \log(\mu_{ij}) = \lambda + \lambda^A_i + \lambda^B_j \]

Model ini menguji hipotesis nol bahwa variabel A (Merokok) dan B (Status) saling independen, yaitu:

\[ H_0: \text{tidak ada asosiasi antara Merokok dan Status} \]

Dalam konteks data Anda:

  • Model saturated akan memodelkan efek penuh dari kombinasi Merokok dan Status.
  • Model independen hanya mempertimbangkan efek Merokok dan Status secara terpisah, tanpa mengakomodasi interaksinya.

Model saturated digunakan sebagai model acuan paling kompleks, sedangkan model independen digunakan untuk pengujian ketergantungan antar variabel melalui perbandingan deviance.

\[ \log(\mu_{ij}) = \mu + \lambda^T_i + \lambda^R_j \]

Model ini menguji hipotesis bahwa variabel X dan Y saling independen.

14.2 Estimasi Parameter

Sistem Persamaan Model Log-Linear:

\[ \begin{aligned} \log(\mu_{11}) &= \lambda + \lambda^A_1 + \lambda^B_1 + \lambda^{AB}_{11} \\ \log(\mu_{12}) &= \lambda + \lambda^A_1 + \lambda^B_2 + \lambda^{AB}_{12} \\ \log(\mu_{21}) &= \lambda + \lambda^A_2 + \lambda^B_1 + \lambda^{AB}_{21} \\ \log(\mu_{22}) &= \lambda + \lambda^A_2 + \lambda^B_2 + \lambda^{AB}_{22} \end{aligned} \]

Constraint Sum-to-Zero:

\[ \lambda^A_1 + \lambda^A_2 = 0 \\ \lambda^B_1 + \lambda^B_2 = 0 \\ \lambda^{AB}_{11} + \lambda^{AB}_{12} + \lambda^{AB}_{21} + \lambda^{AB}_{22} = 0 \]

Rumus Estimasi Parameter dengan Sum-to-Zero Constraint

\[ \lambda^A_1 = \frac{1}{2}[(\log\mu_{11} + \log\mu_{12}) - (\log\mu_{21} + \log\mu_{22})] \\ \lambda^B_1 = \frac{1}{2}[(\log\mu_{11} + \log\mu_{21}) - (\log\mu_{12} + \log\mu_{22})] \\ \lambda^{AB}_{12} = \frac{1}{4}[\log\mu_{12} - \log\mu_{11} - \log\mu_{22} + \log\mu_{21}] \]

14.3 Odds Ratio dan Interpretasi

Odds ratio digunakan untuk mengukur kekuatan hubungan antar dua kategori.

Odds ratio untuk tabel 2x2:

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

Interpretasi nilai OR: - OR = 1: Tidak ada asosiasi - OR > 1: Asosiasi positif - OR < 1: Asosiasi negatif

Log odds ratio (OR):

\[ \begin{aligned} \log(\text{OR}) &= \log\left(\frac{n_{11}n_{22}}{n_{12}n_{21}}\right) \end{aligned} \] Standard Error (SE):

\[ \begin{aligned} \text{SE} &= \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}} \end{aligned} \] 95% Confidence Interval untuk log(OR):

\[ \begin{aligned} \text{CI} &= \log(\text{OR}) \pm 1.96 \times \text{SE} \end{aligned} \] Back-transform to get CI for OR:

\[ \begin{aligned} \text{Lower bound} &= \exp\left[\log(\text{OR}) - z_{\alpha/2} \times \text{SE} \right] \\\\ \text{Upper bound} &= \exp\left[\log(\text{OR}) + z_{\alpha/2} \times \text{SE} \right] \end{aligned} \]

Estimasi odds ratio:

\[ \text{OR} = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} \quad \text{dengan CI: } \left( \exp[\log(\text{OR}) - 1.96 \cdot \text{SE}], \ \exp[\log(\text{OR}) + 1.96 \cdot \text{SE}] \right) \]

14.4 Model Lebih Sederhana dan Perbandingan Model

Perbandingan model independen dan saturated dilakukan menggunakan nilai deviance. Jika p-value < 0.05, model saturated lebih cocok menjelaskan data daripada model independen.

14.5 Studi Kasus: Tingkat Kebahagiaan dan Preferensi Aktivitas

Sebuah lembaga survei sosial ingin mengetahui apakah tingkat kebahagiaan seseorang berkaitan dengan preferensinya terhadap jenis aktivitas di akhir pekan: Aktivitas Sosial (nongkrong, kumpul keluarga) atau Aktivitas Personal (me time, nonton, baca buku). Data diperoleh dari 100 responden:

  • 35 orang dengan kebahagiaan tinggi memilih aktivitas sosial
  • 15 orang dengan kebahagiaan tinggi memilih aktivitas personal
  • 10 orang dengan kebahagiaan rendah memilih aktivitas sosial
  • 40 orang dengan kebahagiaan rendah memilih aktivitas personal

14.5.1 Data Studi Kasus

data <- matrix(c(35, 15, 10, 40), nrow = 2, byrow = TRUE,
               dimnames = list(Kebahagiaan = c("Tinggi", "Rendah"),
                               Aktivitas = c("Sosial", "Personal")))
data
##            Aktivitas
## Kebahagiaan Sosial Personal
##      Tinggi     35       15
##      Rendah     10       40

14.5.2 Tentukan Model: Saturated vs Independen

14.5.2.1 Bentuk Model Saturated (dengan interaksi):

\[ \log(\mu_{ij}) = \lambda + \lambda^A_i + \lambda^B_j + \lambda^{AB}_{ij} \]

14.5.2.2 Bentuk Model Independen (tanpa interaksi):

\[ \log(\mu_{ij}) = \lambda + \lambda^A_i + \lambda^B_j \]

14.5.3 Fitting Model dengan GLM

df <- as.data.frame(as.table(data))
colnames(df) <- c("Kebahagiaan", "Aktivitas", "Freq")

fit_indep <- glm(Freq ~ Kebahagiaan + Aktivitas, family = poisson, data = df)
fit_sat <- glm(Freq ~ Kebahagiaan * Aktivitas, family = poisson, data = df)

14.5.4 Uji Deviance (Perbandingan Model)

anova(fit_indep, fit_sat, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: Freq ~ Kebahagiaan + Aktivitas
## Model 2: Freq ~ Kebahagiaan * Aktivitas
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1         1     26.501                          
## 2         0      0.000  1   26.501 2.634e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretasi: - p-value (0.0000002634) < 0.05 → model saturated lebih baik - Maka, interaksi signifikan → gunakan model saturated untuk interpretasi parameter

14.5.5 Estimasi Parameter (Manual, Sum-to-zero)

14.5.5.1 Rata-rata log-frekuensi sel:

\[ \lambda = \frac{1}{4}[\log(35) + \log(15) + \log(10) + \log(40)] = 3.0637 \]

14.5.5.2 Efek utama Kebahagiaan:

\[ \lambda^A_1 = \frac{1}{2}[(\log(35) + \log(15)) - (\log(10) + \log(40))] = 0.3839 \\ \lambda^A_2 = -0.3839 \]

14.5.5.3 Efek utama Aktivitas:

\[ \lambda^B_1 = \frac{1}{2}[(\log(35) + \log(10)) - (\log(15) + \log(40))] = -0.2695 \\ \lambda^B_2 = 0.2695 \]

14.5.5.4 Efek interaksi:

\[ \lambda^{AB}_{11} = \frac{1}{4}[\log(35) - \log(15) - \log(10) + \log(40)] = 0.6891 \\ \lambda^{AB}_{12} = \lambda^{AB}_{21} = -0.6891, \quad \lambda^{AB}_{22} = 0.6891 \]

14.5.5.5 Interpretasi:

  • λ (intercept): Rata-rata log frekuensi sel = 3.0637
  • λ^A_1: Nilai positif untuk kebahagiaan tinggi menunjukkan bahwa orang bahagia secara keseluruhan muncul lebih banyak dibanding rata-rata, tetapi ini tidak berarti mereka menyukai satu jenis aktivitas tertentu — karena ada interaksi.
  • λ^B_1: Nilai negatif untuk aktivitas sosial berarti bahwa secara umum aktivitas ini terjadi lebih jarang dari rata-rata. Namun, karena interaksi signifikan, efek ini tidak bisa ditafsirkan secara terpisah.
  • λ^{AB}_{ij}:

Bahagia & Sosial (+0.6891): Orang bahagia cenderung memilih aktivitas sosial

Bahagia & Personal (–0.6891): Orang bahagia tidak terlalu memilih aktivitas personal

Tidak bahagia & Sosial (–0.6891): Orang tidak bahagia jarang memilih aktivitas sosial

Tidak bahagia & Personal (+0.6891): Orang tidak bahagia lebih sering memilih aktivitas personal

14.5.6 Uji Asosiasi : Odds Ratio dan Confidence Interval

14.5.6.1 Hitung OR:

\[ \text{OR} = \frac{35 \cdot 40}{15 \cdot 10} = 9.33 \] Maka:

\[ \text{OR} = \frac{35 \cdot 40}{15 \cdot 10} = \frac{1400}{150} = 9.33 \]

Interpretasi:
Orang yang bahagia memiliki kemungkinan 9,33 kali lebih besar untuk memilih aktivitas sosial dibandingkan orang yang tidak bahagia.

14.5.6.2 Log(OR) dan SE:

Log(OR) digunakan dalam model log-linear karena model bekerja di skala logaritma:

\[ \log(\text{OR}) = \log(9.33) = 2.2350 \]

Interpretasi:
Log(OR) setara dengan nilai efek interaksi dalam model log-linear saturated. Nilai positif menandakan bahwa kombinasi “Bahagia – Sosial” lebih sering terjadi dibandingkan jika tidak ada hubungan.

14.5.6.3 Hitung Standard Error (SE)

Standard Error digunakan untuk mengukur ketidakpastian dari estimasi log(OR). Rumusnya:

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

Substitusi:

\[ SE = \sqrt{\frac{1}{35} + \frac{1}{15} + \frac{1}{10} + \frac{1}{40}} = \sqrt{0.0286 + 0.0667 + 0.1000 + 0.0250} = \sqrt{0.2203} = 0.4693 \]

Interpretasi:
Semakin kecil SE, semakin presisi estimasi log(OR). Nilai 0.4693 masih cukup stabil untuk ukuran sampel 100.

14.5.6.4 CI log scale:

\[ 2.2350 \pm 1.96 \cdot 0.5123 = (1.232, 3.238) \]

14.5.6.5 Confidence Interval (log scale)

CI untuk log(OR) dihitung sebagai:

\[ \log(\text{OR}) \pm 1.96 \cdot SE = 2.2350 \pm 1.96 \cdot 0.4693 = 2.2350 \pm 0.9198 \]

Hasilnya:

\[ CI_{\log(\text{OR})} = (1.3152,\ 3.1548) \]

Interpretasi:
CI ini adalah rentang yang kemungkinan besar (95%) mengandung nilai log(OR) yang sebenarnya. Karena seluruh rentangnya positif, maka hubungan yang ditemukan sangat mungkin benar dan bukan kebetulan.

14.5.6.6 Back-transform ke Skala OR

Agar lebih mudah dipahami, kita transformasi kembali CI log(OR) ke skala OR:

\[ CI_{OR} = (\exp(1.3152),\ \exp(3.1548)) = (3.72,\ 23.42) \]

Interpretasi:
Dengan kepercayaan 95%, kita perkirakan bahwa orang bahagia memiliki peluang antara 3,7 hingga 23,4 kali lebih besar untuk memilih aktivitas sosial dibanding orang yang tidak bahagia. Karena CI tidak mencakup 1, hubungan ini signifikan secara statistik.

14.5.6.7 Interpretasi Uji Asosiasi:

  • OR = 9.33 menunjukkan bahwa terdapat asosiasi yang sangat kuat antara tingkat kebahagiaan dan aktivitas yang dipilih.
  • CI log(OR) dan CI OR membuktikan bahwa hubungan ini bukan terjadi secara acak.
  • Temuan ini konsisten dengan hasil model log-linear yang menunjukkan bahwa interaksi signifikan.
  • Maka, preferensi aktivitas akhir pekan seseorang sangat dipengaruhi oleh tingkat kebahagiaannya.

14.5.7 Kesimpulan

  • Model saturated lebih sesuai dibanding model independen, yang ditunjukkan oleh nilai p-value yang sangat kecil (0.0000002634). Ini berarti terdapat interaksi yang signifikan antara kebahagiaan dan aktivitas.

  • Estimasi parameter model menunjukkan bahwa kombinasi “Bahagia – Sosial” dan “Tidak Bahagia – Personal” muncul jauh lebih sering dibanding ekspektasi jika tidak ada hubungan. Sebaliknya, kombinasi “Bahagia – Personal” dan “Tidak Bahagia – Sosial” justru lebih jarang terjadi.

  • Hal ini mengindikasikan bahwa orang bahagia cenderung memilih aktivitas sosial, sementara orang yang kurang bahagia lebih sering memilih aktivitas personal.

  • Hasil ini diperkuat oleh nilai Odds Ratio sebesar 9.33, dengan confidence interval 95% (3.43 – 25.50), yang menunjukkan asosiasi yang sangat kuat dan signifikan secara statistik.

15 Referensi

  1. Agresti, A. (2013). Categorical Data Analysis (3rd ed.). John Wiley & Sons.

  2. Jaya, I. G. N. M. (2024). Analisis Data Kategori. FMIPA UNPAD.

  3. Hosmer, D. W., Lemeshow, S., & Sturdivant, R. X. (2013). Applied Logistic Regression (3rd ed.). John Wiley & Sons.

  4. Christensen, R. (1997). Log-Linear Models and Logistic Regression. Springer.

  5. Fox, J. (2015). Applied Regression Analysis and Generalized Linear Models (3rd ed.). SAGE Publications.

  6. Field, A., Miles, J., & Field, Z. (2012). Discovering Statistics Using R. SAGE Publications.

  7. Kleinbaum, D. G., & Klein, M. (2010). Logistic Regression: A Self-Learning Text (3rd ed.). Springer.

  8. Friendly, M. (2000). Visualizing Categorical Data. SAS Institute.

  9. Kuhn, M., & Johnson, K. (2013). Applied Predictive Modeling. Springer.

  10. Breiman, L. (2001). Random Forests. Machine Learning, 45(1), 5–32. https://doi.org/10.1023/A:1010933404324

  11. Wickham, H., & Grolemund, G. (2016). R for Data Science. O’Reilly Media.

  12. Fox, J., & Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). SAGE Publications.