Ebook Analisis Data Kategori

1. Pendahuluan

Analisis data kategorik adalah cabang statistik yang fokus pada pengolahan, pengambilan kesimpulan, dan interpretasi data yang terdiri dari kategori-kategori atau kelompok diskrit. Data kategorik muncul ketika variabel diukur dalam bentuk kategori seperti jenis kelamin (Laki-laki/Perempuan), status pekerjaan (Bekerja/Tidak), atau jawaban survei (Setuju/Tidak Setuju).

Berbeda dengan data numerik yang bersifat kuantitatif, data kategorik memberikan informasi kualitatif mengenai atribut, kelompok, atau klasifikasi objek atau individu. Oleh karena itu, metode statistika khusus dikembangkan untuk melakukan analisis data ini secara akurat dan efisien.

---

1.1 Tujuan Analisis Data Kategorik

Tujuan utama analisis data kategorik adalah:

• Memahami distribusi frekuensi dan proporsi dari kategori yang diamati.
  
• Menguji hubungan atau asosiasi antara dua atau lebih variabel kategori.
  
• Menilai derajat keberartian statistik hubungan antar variabel kategori.
  
• Membantu pengambilan keputusan berdasarkan karakteristik kelompok dalam data.
  
• Mengembangkan model prediktif atau inferensial dengan data kategori.

---

1.2 Definisi dan Ruang Lingkup

Data kategorik adalah data yang diklasifikasikan ke dalam kelompok atau kategori yang terpisah dan tidak memiliki nilai numerik intrinsik. Data ini terdiri dari:

• Nominal: Kategori tanpa urutan, contohnya warna favorit (merah, hijau, biru).
  
• Ordinal: Kategori dengan urutan, misalnya tingkat pendidikan (SD, SMP, SMA, Perguruan Tinggi).

Analisis data kategorik meliputi teknik statistik yang dirancang untuk menangani variabel jenis ini, seperti:

• Tabel kontingensi
  
• Uji Chi-Square
  
• Uji Fisher Exact
  
• Ukuran asosiasi (odds ratio, relative risk)
• Model regresi logistik dan sebagainya.

Ruang lingkup analisis ini mencakup semua bidang yang menggunakan data kategori untuk evaluasi dan pengambilan keputusan.

---

1.3 Perbedaan dengan Data Kualitatif

Sering kali, istilah data kategorik dan data kualitatif dipakai bergantian, namun ada perbedaan halus:

Aspek	Data Kategorik	Data Kualitatif
Sifat Data	Umumnya terstruktur dalam kategori yang jelas	Bisa berupa deskripsi, narasi, gambar, dll.
Fokus Analisis	Analisis statistik numerik dan komparatif	Analisis tematik, konten, subjektif
Contoh Data	Jenis kelamin, golongan darah, status nikah	Wawancara, observasi, catatan lapangan
Tujuan	Mendapatkan pola, hubungan statistik	Mendapatkan pemahaman mendalam, makna

Data kategorik adalah bagian dari data kuantitatif diskrit, sedangkan data kualitatif lebih bersifat eksploratif dan naratif.

---

1.4 Manfaat Analisis Data Kategorik dalam Berbagai Bidang

Analisis data kategorik memiliki banyak aplikasi penting, antara lain:

• Epidemiologi & Kesehatan Masyarakat: Menentukan faktor risiko penyakit berdasarkan atribut seperti umur, jenis kelamin, kebiasaan hidup.
  
• Sosial dan Psikologi: Memahami hubungan antar variabel sosial, misalnya tingkat pendidikan dan status pekerjaan.
  
• Pemasaran: Segmentasi pasar berdasarkan preferensi konsumen, perilaku membeli, dan demografi.
  
• Penelitian Klinis: Mengevaluasi efektivitas pengobatan dengan data outcome kategori (sukses/gagal).
  
• Pendidikan: Analisis hasil ujian atau survei yang dikategorikan seperti lulus/tidak lulus, tingkat kepuasan.
  
• Manajemen & Bisnis: Pengambilan keputusan berbasis kategori pelanggan, produk, atau performa.

Analisis yang tepat memungkinkan peneliti, praktisi, dan pengambil keputusan memahami pola, membuat prediksi, dan merancang intervensi yang efektif.

---

Kesimpulan

Analisis data kategorik merupakan komponen vital dalam statistika yang membantu mengubah data kategori menjadi informasi bermakna. Penguasaan konsep dan teknik analisis ini sangat penting untuk keberhasilan riset dan pengembangan di berbagai bidang ilmu dan praktik profesional.

2. Metode Dasar dalam Analisis Data Kategorik

2.1 Tabel Kontingensi

Tabel kontingensi adalah cara paling dasar untuk menampilkan frekuensi gabungan dari dua atau lebih variabel kategori.
Biasanya berbentuk tabel dua dimensi seperti tabel 2×2, 3×2, dan seterusnya, yang menggambarkan distribusi joint frequency antar kategori.

2.2 Ukuran Frekuensi dan Proporsi

• Frekuensi Absolut: Jumlah kasus di setiap kategori (contoh: 20 laki-laki, 30 perempuan).
  
• Proporsi dan Persentase: Menunjukkan bagian relatif dari total, memudahkan perbandingan antar kelompok.

2.3 Uji Chi-Square (Chi-Square Test)

Digunakan untuk menguji hipotesis korelasi atau asosiasi antara dua variabel kategori dalam tabel kontingensi.
Hipotesis nol menyatakan tidak ada asosiasi (independen).

• Kelebihan: Mudah digunakan, cocok untuk sampel besar.
  
• Keterbatasan: Tidak akurat jika frekuensi sel kecil (<5).

2.4 Uji Fisher’s Exact

Alternatif untuk Chi-Square saat ukuran sampel kecil atau frekuensi sel rendah, memberikan probabilitas exact tanpa asumsi asimptotik.
Sering digunakan dalam tabel 2×2 dengan data kecil.

2.5 Ukuran Asosiasi

Memperkuat analisis dengan mengukur seberapa kuat hubungan antara variabel kategori:

• Odds Ratio (OR)
  
• Relative Risk (RR)
  
• Risk Difference (RD)
  
• Koefisien Kontingensi, Phi, dan Cramér’s V (khusus untuk tabel lebih besar)

2.6 Regresi Logistik

Metode modeling untuk memprediksi variabel dependen kategori (biasanya biner) dari satu atau beberapa variabel independen (kategori maupun numerik).
Berguna untuk analisis multivariat dan kontrol variabel pengganggu.

---

2.7 Metode Visualisasi Data Kategorik

• Bar Plot (Diagram Batang): Visualisasi frekuensi atau proporsi kategori.
  
• Mosaic Plot: Visualisasi distribusi bersama (joint distribution) dalam tabel kontingensi.
  
• Pie Chart: Menunjukkan proporsi kategori, meskipun kurang direkomendasikan untuk perbandingan.
  
• Heatmap Tabel Kontingensi: Menunjukkan intensitas frekuensi di sel tabel.

---

2.8 Contoh Alur Analisis Data Kategorik

1. Eksplorasi Data: Hitung frekuensi dan proporsi kategori.
   
2. Tabulasi Bivariat: Buat tabel kontingensi antar variabel.
   
3. Uji Asosiasi: Lakukan uji Chi-Square atau Fisher’s Exact.
   
4. Ukuran Asosiasi: Hitung OR, RR, atau indeks kekuatan hubungan lainnya.
   
5. Modeling (Jika Perlu): Terapkan regresi logistik untuk analisis multivariat.
   
6. Visualisasi: Gunakan plot untuk membantu interpretasi dan komunikasi hasil.

3. Distribusi Probabilitas dalam Data Kategori

Distribusi probabilitas dalam data kategori menjelaskan peluang tiap kategori muncul dalam suatu variabel diskrit. Pemahaman terhadap berbagai distribusi probabilitas sangat penting dalam analisis data kategori, terutama untuk pemodelan dan inferensi statistik.

3.1 Distribusi Bernoulli

Distribusi Bernoulli adalah distribusi probabilitas diskrit paling sederhana, digunakan untuk variabel dengan dua kategori (biner), misal sukses/gagal, ya/tidak.

Definisi:

Untuk variabel acak \(X\) dengan nilai 1 (sukses) dengan peluang \(p\), dan nilai 0 (gagal) dengan peluang \(1-p\):

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

Karakteristik utama:

• Rata-rata: \(E(X) = p\)
  
• Varians: \(\text{Var}(X) = p(1-p)\)

Contoh:

Keluarnya angka kepala (success) saat melempar koin: \(p = 0.5\).

---

3.2 Distribusi Binomial

Distribusi Binomial memperluas Bernoulli ke jumlah percobaan \(n\) independen dengan peluang sukses sama \(p\), menghitung peluang \(k\) sukses dalam \(n\) percobaan.

Fungsi probabilitas:

\[ P(X = k) = \binom{n}{k} p^{k} (1-p)^{n-k}, \quad k = 0,1,\ldots,n \]

Karakteristik utama:

• Rata-rata: \(E(X) = np\)
  
• Varians: \(\text{Var}(X) = np(1-p)\)

Contoh:

Jumlah orang yang memilih produk A dari 10 pelanggan, dengan probabilitas memilih 0.3.

---

3.3 Distribusi Multinomial

Distribusi Multinomial adalah generalisasi distribusi Binomial untuk variabel kategori dengan lebih dari dua kemungkinan kategori.

Fungsi probabilitas:

Jika ada \(K\) kategori dengan peluang masing-masing \(p_1, p_2, \ldots, p_K\) dan \(n\) percobaan, peluang observasi hasil dengan frekuensi \((n_1, n_2, \ldots, n_K)\) adalah:

\[ P(n_1, n_2, \ldots, n_K) = \frac{n!}{n_1! n_2! \cdots n_K!} p_1^{n_1} p_2^{n_2} \cdots p_K^{n_K} \]

dengan \(\sum_{k=1}^K n_k = n\) dan \(\sum_{k=1}^K p_k = 1\).

Contoh:

Jumlah pelanggan memilih jenis transportasi dari tiga kategori: Mobil, Motor, Sepeda.

---

3.4 Distribusi Poisson

Distribusi Poisson digunakan untuk menghitung peluang jumlah kejadian dalam interval waktu atau ruang tertentu ketika kejadian tersebut jarang dan independen.

Fungsi probabilitas:

\[ P(X = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0,1,2, \ldots \]

dengan \(\lambda\) adalah rata-rata jumlah kejadian per interval.

Karakteristik utama:

• Rata-rata: \(E(X) = \lambda\)
  
• Varians: \(\text{Var}(X) = \lambda\)

Contoh:

Jumlah kecelakaan lalu lintas dalam sehari di sebuah persimpangan.

---

Contoh Implementasi Distribusi di R

# Distribusi Bernoulli: Probabilitas sukses = 0.4  
p <- 0.4  
dbinom(1, size = 1, prob = p)  # P(X=1)  

## [1] 0.4

dbinom(0, size = 1, prob = p)  # P(X=0)  

## [1] 0.6

# Distribusi Binomial: n=10 percobaan, p=0.3 sukses  
n <- 10  
p <- 0.3  
k <- 3  
dbinom(k, size = n, prob = p)  # P(X=3)  

## [1] 0.2668279

# Distribusi Multinomial: n=5, p kategori = c(0.2, 0.5, 0.3)  
n <- 5  
prob_kat <- c(0.2, 0.5, 0.3)  
n_kat <- c(1, 3, 1)  # jumlah observasi dalam tiap kategori  
dmultinom(n_kat, size = n, prob = prob_kat)  # P(X1=1,X2=3,X3=1)  

## [1] 0.15

# Distribusi Poisson: Lambda = 2, probabilitas 3 kejadian  
lambda <- 2  
k <- 3  
dpois(k, lambda = lambda)  # P(X=3)  

## [1] 0.180447

Berikut adalah versi yang sudah diperbaiki dan dimodifikasi agar lebih sistematis dan mudah dibaca, tetap lengkap, serta sesuai standar penulisan ilmiah:

---

4. Desain Sampling dalam Analisis Data Kategori

Dalam analisis data kategori, desain sampling memegang peran krusial dalam menentukan validitas dan reliabilitas hasil penelitian. Pemilihan desain sampling yang tepat sangat dipengaruhi oleh tujuan penelitian dan jenis data yang dikumpulkan. Secara umum, desain sampling dapat dikelompokkan ke dalam dua pendekatan utama, yaitu:

• Prospective Sampling
  
• Retrospective Sampling

Setiap pendekatan memiliki karakteristik dan metode sampling yang berbeda, yang sering kali digunakan baik dalam penelitian eksperimental maupun observasional seperti eksperimen, studi kohort, dan studi kasus-kontrol.

---

4.1 Prospective Sampling

Prospective sampling adalah metode pengambilan sampel di mana subjek penelitian diidentifikasi dan diikuti dalam suatu periode waktu tertentu untuk memantau perkembangan variabel yang diteliti. Pendekatan ini memungkinkan peneliti mengontrol variabel bebas sebelum pengukuran hasil dilakukan, sehingga banyak digunakan dalam studi kausal dan eksperimental. Beberapa jenis desain sampling dalam pendekatan ini adalah:

4.1.1 Eksperimen

Dalam studi eksperimen, subjek secara acak dialokasikan ke dalam kelompok perlakuan dan kontrol untuk mengendalikan variabel pengganggu. Teknik sampling yang umum digunakan meliputi:

• Simple Random Sampling (SRS): Setiap individu dalam populasi memiliki probabilitas yang sama untuk dipilih.
• Stratified Random Sampling: Populasi dibagi menjadi strata berdasarkan karakteristik tertentu, lalu diambil sampel secara acak dari setiap strata.
• Cluster Sampling: Populasi dibagi menjadi kelompok-kelompok (cluster), kemudian beberapa cluster dipilih secara acak untuk dianalisis.

4.1.2 Studi Kohort

Studi kohort adalah penelitian observasional yang mengikuti kelompok individu dengan karakteristik tertentu dari waktu ke waktu untuk mengamati kejadian yang dipelajari. Teknik sampling dalam studi kohort antara lain:

• Census Sampling: Seluruh anggota populasi tertentu diikutsertakan dalam penelitian.
• Systematic Sampling: Subjek dipilih berdasarkan interval tertentu dari daftar populasi.
• Matched Sampling: Setiap individu dalam kelompok kohort dipasangkan dengan individu serupa dalam kelompok lain berdasarkan variabel tertentu.

---

4.2 Retrospective Sampling

Retrospective sampling mengumpulkan data dari peristiwa yang telah terjadi sebelumnya. Metode ini umum digunakan pada studi observasional untuk mengkaji hubungan antara faktor risiko dan hasil yang sudah terjadi.

4.2.1 Studi Kasus-Kontrol

Dalam studi kasus-kontrol, sekelompok individu dengan kondisi tertentu (kasus) dibandingkan dengan kelompok tanpa kondisi tersebut (kontrol). Teknik sampling umum meliputi:

• Purposive Sampling: Pemilihan sampel berdasarkan karakteristik yang relevan dengan tujuan penelitian.
• Snowball Sampling: Responden awal merekrut subjek lain yang memiliki karakteristik serupa.
• Incidence Density Sampling: Kasus dan kontrol dipilih dari populasi yang sama dengan memperhitungkan periode waktu munculnya kasus.

4.2.2 Studi Kohort Retrospektif

Dalam studi kohort retrospektif, data historis digunakan untuk mengelompokkan individu berdasarkan paparan, kemudian dianalisis hasil yang terjadi. Teknik sampling yang sering digunakan adalah:

• Convenience Sampling: Subjek dipilih berdasarkan ketersediaan data.
• Quota Sampling: Sampel dipilih untuk mencerminkan proporsi tertentu dalam populasi.
• Case-Based Sampling: Sampel diambil berdasarkan karakteristik kasus yang telah terjadi.

4.3 Tabel Perbandingan Desain Sampling

Jenis Studi	Pendekatan	Metode Sampling	Keuntungan	Kelemahan
Eksperimen	Prospektif	SRS, Stratified, Cluster	Kontrol tinggi terhadap variabel, dapat analisis sebab-akibat	Biaya tinggi, etika dan validitas perlu diperhatikan
Studi Kohort	Prospektif	Census, Systematic, Matched	Dapat mengamati perkembangan kejadian dalam jangka panjang	Membutuhkan waktu lama, risiko kehilangan partisipan
Studi Kasus-Kontrol	Retrospektif	Purposive, Snowball, Incidence Density	Mudah dan cepat, efisien untuk penyakit langka	Sulit mengontrol variabel pengganggu, bias recall rentan
Studi Kohort Retrospektif	Retrospektif	Convenience, Quota, Case-Based	Memanfaatkan data historis, lebih murah dari studi kohort	Ketergantungan pada kualitas data historis, risiko missing data

---

Desain sampling dalam analisis data kategori sangat bergantung pada pendekatan yang digunakan, apakah prospective atau retrospective. Pemilihan metode sampling yang tepat dalam eksperimen, studi kohort, dan studi kasus-kontrol dapat meningkatkan validitas penelitian serta memastikan hasil yang dapat digeneralisasikan. Pemahaman mendalam mengenai karakteristik tiap metode sampling penting dalam merancang penelitian yang robust dan berkualitas.

---

5 Tabel Kontingensi

Tabel kontingensi adalah tabel yang digunakan untuk menyajikan distribusi frekuensi dari dua atau lebih variabel kategori. Tabel ini membantu dalam memahami hubungan antara variabel kategori dengan menampilkan jumlah kejadian dalam setiap kombinasi kategori.

Terdapat beberapa jenis Tabel Kontingensi

1. Tabel 2×2 Tabel ini memiliki dua variabel dengan dua kategori masing-masing. Ini adalah bentuk paling sederhana dari tabel kontingensi. Contoh:Hubungan antara merokok (Ya/Tidak) dan kejadian penyakit paru-paru (Ya/Tidak).

2. Tabel RxC (General) Tabel ini memiliki R baris dan C kolom, yang memungkinkan lebih banyak kategori. Contoh:Hubungan antara tingkat pendidikan (SD, SMP, SMA, S1, S2) dan status pekerjaan (Pengangguran, Pekerja, Wiraswasta).

3. Tabel Kontingensi Berganda Digunakan untuk lebih dari dua variabel kategori. Contoh:Hubungan antara jenis kelamin, status pekerjaan, dan tingkat pendidikan.

5.1 Pengertian Tabel Kontingensi 2x2

Tabel kontingensi adalah tabel yang digunakan untuk menampilkan frekuensi kejadian dua variabel kategori sekaligus. Tabel kontingensi 2x2 menggambarkan hubungan antara dua variabel kategori yang masing-masing memiliki dua kategori.

5.2 Struktur Tabel Kontingensi 2x2

	Kategori 1 (Variabel B)	Kategori 2 (Variabel B)	Jumlah Baris
Kategori 1 (Variabel A)	\(a\)	\(b\)	\(a+b\)
Kategori 2 (Variabel A)	\(c\)	\(d\)	\(c+d\)
Jumlah Kolom	\(a+c\)	\(b+d\)	\(n=a+b+c+d\)

---

Simbol dan Notasi

• \(a\), \(b\), \(c\), dan \(d\) adalah frekuensi atau jumlah kasus pada setiap sel.
  
• \(n\) adalah total seluruh observasi.
  
• Baris dan kolom menunjukkan kategori dari dua variabel berbeda.

---

Fungsi dan Kegunaan

• Untuk mengetahui distribusi frekuensi dari dua variabel kategori.
  
• Sebagai dasar pengujian statistik seperti uji Chi-Square untuk independensi.
  
• Menilai hubungan atau asosiasi antara dua variabel kategori.

---

5.3 Distribusi Peluang dalam Tabel Kontingensi 2 × 2

---

5.3.1 Peluang Bersama

Peluang bersama adalah probabilitas bahwa kedua variabel kategori terjadi secara bersamaan dalam suatu sel tabel kontingensi. Secara matematis, peluang bersama antara kejadian \(A_i\) (misalnya kategori baris) dan \(B_j\) (kategori kolom) dihitung sebagai:

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

dimana:
- \(n_{ij}\) adalah frekuensi pengamatan pada sel di baris \(i\) dan kolom \(j\).
- \(n\) adalah total seluruh pengamatan dalam tabel.

---

5.3.2 Peluang Marginal

Peluang marginal adalah probabilitas kejadian suatu variabel tanpa mempertimbangkan variabel lain. Peluang marginal dapat dihitung untuk baris maupun kolom.

• Peluang marginal baris:

\[ P(A_i) = \frac{n_{i \cdot}}{n} \]

dimana \(n_{i \cdot}\) adalah jumlah total pada baris ke-\(i\).

• Peluang marginal kolom:

\[ P(B_j) = \frac{n_{\cdot j}}{n} \]

dimana \(n_{\cdot j}\) adalah jumlah total pada kolom ke-\(j\).

---

5.3.3 Peluang Bersyarat

Peluang bersyarat adalah probabilitas suatu kejadian terjadi dengan syarat kejadian lain telah terjadi. Dalam konteks tabel kontingensi, peluang bersyarat variabel kolom \(B_j\) diberikan baris \(A_i\) adalah:

\[ P(B_j \mid A_i) = \frac{P(A_i, B_j)}{P(A_i)} = \frac{n_{ij}}{n_{i \cdot}} \]

---

Contoh Perhitungan Manual

Tabel kontingensi di bawah ini menunjukkan data hubungan antara Kebiasaan Berjalan Kaki (Sering, Jarang) dan Hipertensi (Ya, Tidak).

Berjalan Kaki / Hipertensi	Ya	Tidak	Total
Sering	20	80	100
Jarang	50	50	100
Total	70	130	200

---

Hitung Peluang Bersama \(P(A_i, B_j)\):

Sel	Frekuensi \(n_{ij}\)	Peluang Bersama \(P(A_i, B_j) = \frac{n_{ij}}{n}\)
(Sering, Ya)	20	20 / 200 = 0.10
(Sering, Tidak)	80	80 / 200 = 0.40
(Jarang, Ya)	50	50 / 200 = 0.25
(Jarang, Tidak)	50	50 / 200 = 0.25

---

Hitung Peluang Marginal Baris \(P(A_i)\):

Kategori Berjalan Kaki	Jumlah Baris \(n_{i\cdot}\)	Peluang \(P(A _i) = \frac{n_{i \cdot}}{n}\)
Sering	100	\(\frac{100}{200} = 0.5\)
Jarang	100	\(\frac{100}{200} = 0.5\)

---

Hitung Peluang Marginal Kolom \(P(B_j)\):

Kategori Hipertensi	Jumlah Kolom \(n_{\cdot j}\)	Peluang \(P( B _j) = \frac{n_{\cdot j}}{n}\)
Ya	70	\(\frac{70}{200} = 0.35\)
Tidak	130	\(\frac{130}{200} = 0.65\)

---

Hitung Peluang Bersyarat \(P(B_j \mid A_i)\):

Kategori Berjalan Kaki	Hipertensi Ya	Hipertensi Tidak
Sering	\(\ f r ac{20}{100} = 0.20\)	\(\ f r ac{80}{100} = 0.80\)
Jarang	\(\ f r ac{50}{100} = 0.50\)	\(\ f r ac{50}{100} = 0.50\)

---

Interpretasi

• Peluang bersyarat menunjukkan bahwa peluang seseorang memiliki hipertensi jika berjalan kaki sering adalah 20%, sedangkan jika berjalan kaki jarang, peluangnya meningkat menjadi 50%.
  
• Ini mengindikasikan bahwa kebiasaan berjalan kaki mungkin berpengaruh terhadap risiko hipertensi.
  
• Peluang marginal memberikan gambaran proporsi masing-masing kategori dalam populasi keseluruhan.

---

Implementasi Perhitungan di R

# Membuat tabel kontingensi 2x2  
tabel <- matrix(c(20, 80, 50, 50), nrow = 2, byrow = TRUE)  
rownames(tabel) <- c("Sering", "Jarang")  
colnames(tabel) <- c("Ya", "Tidak")  

print("Tabel Kontingensi:")  

## [1] "Tabel Kontingensi:"

print(tabel)  

##        Ya Tidak
## Sering 20    80
## Jarang 50    50

# Total observasi  
n <- sum(tabel)  

# Peluang Bersama  
peluang_bersama <- tabel / n  

cat("\nPeluang Bersama (Joint Probability):\n")  

## 
## Peluang Bersama (Joint Probability):

print(round(peluang_bersama, 3))  

##          Ya Tidak
## Sering 0.10  0.40
## Jarang 0.25  0.25

# Peluang Marginal Baris  
peluang_marginal_baris <- rowSums(tabel) / n  

cat("\nPeluang Marginal Baris:\n")  

## 
## Peluang Marginal Baris:

print(round(peluang_marginal_baris, 3))  

## Sering Jarang 
##    0.5    0.5

# Peluang Marginal Kolom  
peluang_marginal_kolom <- colSums(tabel) / n  

cat("\nPeluang Marginal Kolom:\n")  

## 
## Peluang Marginal Kolom:

print(round(peluang_marginal_kolom, 3))  

##    Ya Tidak 
##  0.35  0.65

# Peluang Bersyarat (Baris)  
peluang_bersyarat <- prop.table(tabel, margin = 1)  

cat("\nPeluang Bersyarat (Conditional Probability P(Hipertensi | Berjalan Kaki)):\n")  

## 
## Peluang Bersyarat (Conditional Probability P(Hipertensi | Berjalan Kaki)):

print(round(peluang_bersyarat, 3))

##         Ya Tidak
## Sering 0.2   0.8
## Jarang 0.5   0.5

5.4 Ukuran Asosiasi

Beda Peluang (Risk Difference)

Beda peluang mengukur perbedaan probabilitas kejadian di antara dua kelompok:

\[ RD = P(\text{Hipertensi | Jarang Berjalan}) - P(\text{Hipertensi | Sering Berjalan}) \]

\[ RD = 0.50 - 0.20 = 0.30 \]

Interpretasi: Individu yang jarang berjalan kaki memiliki kemungkinan hipertensi lebih tinggi sebesar 30% dibandingkan individu yang sering berjalan kaki.

Risiko Relatif (Relative Risk)

Risiko relatif mengukur rasio kemungkinan kejadian antara dua kelompok:

\[ RR = \frac{P(\text{Hipertensi | Jarang Berjalan})}{P(\text{Hipertensi | Sering Berjalan})} \]

\[ RR = \frac{0.50}{0.20} = 2.5 \]

Interpretasi: Risiko hipertensi pada individu yang jarang berjalan kaki 2.5 kali lebih tinggi dibandingkan dengan individu yang sering berjalan kaki.

Odd Ratio (OR)

Odd Ratio (OR) membandingkan peluang kejadian relatif antara dua kelompok:

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

Di mana: - \(a = 50\) (Jarang Berjalan dan Hipertensi) - \(b = 50\) (Jarang Berjalan dan Tidak Hipertensi) - \(c = 20\) (Sering Berjalan dan Hipertensi) - \(d = 80\) (Sering Berjalan dan Tidak Hipertensi)

\[ OR = \frac{\frac{50}{50}}{\frac{20}{80}} = \frac{1}{0.25} = 4 \]

Interpretasi: Individu yang jarang berjalan kaki memiliki kemungkinan 4 kali lebih besar untuk mengalami hipertensi dibandingkan individu yang sering berjalan kaki.

Kesimpulan

Dari hasil perhitungan, dapat disimpulkan bahwa:

1. Peluang seseorang mengalami hipertensi lebih tinggi jika mereka jarang berjalan kaki (50%) dibandingkan dengan mereka yang sering berjalan kaki (20%).

2. Beda peluang menunjukkan bahwa individu yang jarang berjalan memiliki risiko hipertensi lebih tinggi sebesar 30%.

3. Risiko relatif menunjukkan bahwa individu yang jarang berjalan memiliki risiko hipertensi 2.5 kali lebih besar dibandingkan yang sering berjalan.

4. Odd Ratio menunjukkan bahwa individu yang jarang berjalan kaki memiliki kemungkinan 4 kali lebih besar mengalami hipertensi dibandingkan yang sering berjalan.

5. Kebiasaan berjalan kaki tampaknya memiliki hubungan dengan kejadian hipertensi, di mana individu yang lebih sering berjalan kaki memiliki risiko hipertensi yang lebih rendah.

Dengan memahami ukuran asosiasi ini, kita dapat menginterpretasikan hubungan antara kebiasaan berjalan kaki dan risiko hipertensi secara lebih mendalam.

6 Inferensi Tabel Kontingensi Dua Arah

Inferensi dalam statistik mengacu pada proses pengambilan kesimpulan mengenai populasi berdasarkan sampel data. Dalam konteks tabel kontingensi dua arah, inferensi digunakan untuk menganalisis hubungan antara dua variabel kategorikal yang disusun dalam tabel kontingensi.

Tabel kontingensi adalah tabel yang menyajikan distribusi frekuensi dari dua variabel kategorikal dalam bentuk matriks. Tujuan utama dari inferensi pada tabel ini adalah untuk memahami hubungan antara variabel-variabel tersebut melalui estimasi dan pengujian hipotesis.

6.1 Estimasi

Estimasi bertujuan untuk memperkirakan parameter populasi berdasarkan data sampel. Estimasi dibagi menjadi dua jenis:

6.1.1 Estimasi Titik

Estimasi titik digunakan untuk menentukan satu nilai spesifik sebagai perkiraan terbaik dari parameter populasi. Rumus untuk estimasi titik proporsi adalah:

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

Keterangan: - \(\hat{p}\): estimasi titik proporsi, - \(x\): jumlah individu dalam kategori tertentu, - \(n\): total jumlah individu dalam sampel.

6.1.2 Estimasi Interval

Estimasi interval memberikan rentang nilai yang diyakini mengandung parameter populasi dengan tingkat kepercayaan tertentu. Rumusnya:

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

Keterangan: - \(Z_{\alpha/2}\): nilai dari distribusi normal standar untuk tingkat kepercayaan tertentu (misal, 1.96 untuk 95% confidence level), - \(\hat{p}\): estimasi titik proporsi, - \(n\): ukuran sampel.

6.3 Uji Hipotesis

6.3.1 Uji Proporsi

Uji proporsi dua sampel digunakan untuk menguji apakah ada perbedaan signifikan antara dua proporsi populasi.

a. Rumus Uji Proporsi Dua Sampel

Estimasi proporsi dalam masing-masing kelompok diberikan oleh: \[ \hat{p_1} = \frac{n_{11}}{n_1.},\ \hat{p_2} = \frac{n_{21}}{n_2.} \]

Estimasi proporsi gabungan (pooling proportion): \[ \hat{p} = \frac{n_{11}+n_{21}}{n_1.+n_2.} \]

Statistik uji untuk uji proporsi dua sampel:

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

Perhitungan

# Data Observasi
data <- matrix(c(20, 80, 50, 50), nrow = 2, byrow = TRUE)
dimnames(data) <- list(Berjalan_Kaki = c("Sering", "Jarang"), Hipertensi = c("Ya", "Tidak"))
print(data)

##              Hipertensi
## Berjalan_Kaki Ya Tidak
##        Sering 20    80
##        Jarang 50    50

p1 <- data[1,1] / sum(data[1,])  # Proporsi hipertensi pada kelompok sering berjalan kaki
p2 <- data[2,1] / sum(data[2,])  # Proporsi hipertensi pada kelompok jarang berjalan kaki
p_hat <- (data[1,1] + data[2,1]) / sum(data)  # Proporsi gabungan

SE_pooled <- sqrt(p_hat * (1 - p_hat) * (1/sum(data[1,]) + 1/sum(data[2,])))
Z_stat <- (p1 - p2) / SE_pooled

cat("Proporsi sering berjalan kaki dengan hipertensi:", p1, "\n")

## Proporsi sering berjalan kaki dengan hipertensi: 0.2

cat("Proporsi jarang berjalan kaki dengan hipertensi:", p2, "\n")

## Proporsi jarang berjalan kaki dengan hipertensi: 0.5

cat("Z Statistik:", Z_stat, "\n")

## Z Statistik: -4.447496

alpha <- 0.05
p_value <- 2 * (1 - pnorm(abs(Z_stat)))  # Uji dua sisi
Z_critical <- qnorm(1 - alpha/2)

cat("Z kritis:", Z_critical, "\n")

## Z kritis: 1.959964

cat("P-value:", p_value, "\n")

## P-value: 8.687712e-06

if(abs(Z_stat) > Z_critical) {
  cat("Keputusan: Tolak H0, terdapat perbedaan proporsi hipertensi yang signifikan antara kelompok.\n")
} else {
  cat("Keputusan: Gagal menolak H0, tidak ada perbedaan signifikan dalam proporsi hipertensi antara kelompok.\n")
}

## Keputusan: Tolak H0, terdapat perbedaan proporsi hipertensi yang signifikan antara kelompok.

6.3.2 Uji Asosiasi

Uji asosiasi digunakan untuk mengukur kekuatan dan arah hubungan antara dua variabel kategorik. Dengan tiga ukuran utama yang diujikan, yakni Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR)

Hipotesis

• H₀ : Kedua variabel tidak saling berhubungan (independen).

• H₁ : Kedua variabel saling berhubungan (tidak independen)

Statistik Uji: Z

Untuk statistik Risk Difference (RD): Rumus: \[RD = p_1 - p_2, \text{ dengan } p_1 = \frac{n_{11}}{n_{1.}}, \quad p_2 = \frac{n_{21}}{n_{2.}}\]

Standard Error:

\[SE_{RD} = \sqrt{\frac{p_1(1-p_1)}{n_{1.}} + \frac{p_2(1-p_2)}{n_{2.}}}\]

Z-statistik:

\[Z_{RD} = \frac{RD}{SE_{RD}}\]

Didapatkan

• Grup 1 → \[p_1 = \frac{50}{80} = 0.625\] • Grup 2 → \[p_2 = \frac{30}{80} = 0.375\]

Risk Difference:

\[RD = p_1 - p_2 = 0.625 - 0.375 = 0.25\]

Standard Error RD:

\[SE_{RD} = \sqrt{\frac{0.625(1-0.625)}{80} + \frac{0.375(1-0.375)}{80}}\] \[= \sqrt{\frac{0.234375}{80} + \frac{0.234375}{80}} = \sqrt{\frac{0.46875}{80}} \approx \sqrt{0.00586} \approx 0.0766\]

Z-statistik:

\[Z_{RD} = \frac{RD}{SE_{RD}} = \frac{0.25}{0.0766} \approx 3.26\]

Untuk statistik Relative Risk(RR):

• Rumus: \[RR = \frac{p_1}{p_2}\] • Standard Error (dari log RR): \[SE_{ln \ RR} = \sqrt{(\frac{1}{n_{11}} - \frac{1}{n_{1.}})+(\frac{1}{n_{21}} - \frac{1}{n_{2.}})}\] • Z-statistik: \[Z_{ln \ RR} = \frac{ln(RR)}{SE_{ln \ RR}}\] Didapatakan \[RR = \frac{p_1}{p_2} = \frac{0.625}{0.375} = \boxed{1.\overline{6}} \approx 1.67\] Standard Error log(RR): \[SE_{ln\, RR} = \sqrt{\left(\frac{1}{50} - \frac{1}{80}\right) + \left(\frac{1}{30} - \frac{1}{80}\right)}\] \[= \sqrt{(0.02 - 0.0125) + (0.0333 - 0.0125)} = \sqrt{0.0283} \approx \boxed{0.168}\] Z-statistik: \[Z_{ln\, RR} = \frac{ln(1.67)}{0.168} \approx \frac{0.511}{0.168} \approx \boxed{3.04}\] Untuk statistik Odds Ratio • Rumus: \[OR = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}}\] • Standard Error (dari log OR): \[SE_{ln\,OR} = \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}}\] • Z-statistik: \[Z_{ln\,OR} = \frac{ln(OR)}{SE_{ln\,OR}}\] Didapatkan \[OR = \frac{50 \times 50}{30 \times 30} = \frac{2500}{900} = 2.78\] Standard Error log(OR): \[SE_{ln \ OR} = \sqrt{\frac{1}{50} + \frac{1}{30} + \frac{1}{30} + \frac{1}{50}} = \sqrt{0.02 + 0.0333 + 0.0333 + 0.02} = \sqrt{0.1066} \approx 0.3265\] Z-statistik: \[Z_{ln \ OR} = \frac{ln(2.78)}{0.3265} \approx \frac{1.022}{0.3265} \approx 3.13\] Perhitungan dengan software R

# Input data
n11 <- 50  # Grup 1, Kejadian
n12 <- 30  # Grup 1, Tidak Kejadian
n21 <- 30  # Grup 2, Kejadian
n22 <- 50  # Grup 2, Tidak Kejadian

# Marginal total
n1. <- n11 + n12  # Total Grup 1
n2. <- n21 + n22  # Total Grup 2

# Risk Difference (RD)
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)
rr <- p1 / p2
se_ln_rr <- sqrt((1 / n11) - (1 / n1.) + (1 / n21) - (1 / n2.))
z_rr <- log(rr) / se_ln_rr

# Odds Ratio (OR)
or <- (n11 * n22) / (n12 * n21)
se_ln_or <- sqrt(1 / n11 + 1 / n12 + 1 / n21 + 1 / n22)
z_or <- log(or) / se_ln_or

# Output hasil
list(
  Risk_Difference = rd,
  SE_RD = se_rd,
  Z_RD = z_rd,

  Relative_Risk = rr,
  SE_Log_RR = se_ln_rr,
  Z_Log_RR = z_rr,

  Odds_Ratio = or,
  SE_Log_OR = se_ln_or,
  Z_Log_OR = z_or
)

## $Risk_Difference
## [1] 0.25
## 
## $SE_RD
## [1] 0.07654655
## 
## $Z_RD
## [1] 3.265986
## 
## $Relative_Risk
## [1] 1.666667
## 
## $SE_Log_RR
## [1] 0.1683251
## 
## $Z_Log_RR
## [1] 3.034756
## 
## $Odds_Ratio
## [1] 2.777778
## 
## $SE_Log_OR
## [1] 0.3265986
## 
## $Z_Log_OR
## [1] 3.128155

Titik Kritis: Dengan α = 0.05, maka Ztabel= ±1,96

Kriteria Uji: Tolak H₀ apabila, Z hitung > |Ztabel|

Keputusan :

asosiasi_data <- data.frame(
  Ukuran = c("Risk Difference (RD)", "Relative Risk (RR)", "Odds Ratio (OR)"),
  Nilai = c(0.25, 1.67, 2.78),
  Z_Statistik = c(3.26, 3.04, 3.13),
  Signifikan = c("Ya", "Ya", "Ya")
)

print(asosiasi_data, row.names = FALSE)

##                Ukuran Nilai Z_Statistik Signifikan
##  Risk Difference (RD)  0.25        3.26         Ya
##    Relative Risk (RR)  1.67        3.04         Ya
##       Odds Ratio (OR)  2.78        3.13         Ya

Kesimpulan: 1. Terdapat perbedaan risiko absolut sebesar 25% antara kelompok sering berjalan dan jarang berjalan terhadap hipertensi. 2.Kelompok sering berjalan memiliki risiko 1.67 kali lebih besar dibanding kelompok 2. 3. Odds kejadian pada kelompok 1 adalah 2.78 kali dibanding kelompok 2.

6.3.3 Uji Independensi

Tujuan: Mengetahui apakah dua variabel kategori saling bebas (independen) atau tidak (berhubungan).

Hipotesis

-   **H₀** : Kedua variabel **tidak saling berhubungan**
    (independen).

-   **H₁** : Kedua variabel **saling berhubungan** (tidak
    independen)

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

Dengan \[E_{ij} = \frac{(total\:baris) \times (total\:kolom)}{n}\]

Didapatkan: \[E_{(Sering,Ya)} = \frac{100 \times 70}{200} = 35\] \[E_{(Sering,Tidak)} = \frac{100 \times 130}{200} = 65\] \[E_{(Jarang,Ya)} = \frac{100 \times 70}{200} = 35\] \[E_{(Jarang,Tidak)} = \frac{100 \times 130}{200} = 65\]

Nilai Chi-Square: \[\chi^2 = \frac{(20-35)^2}{35} + \frac{(80-65)^2}{65} + \frac{(50-35)^2}{35} + \frac{(50-65)^2}{65}\] \[ \chi^2 = \frac{225}{35} + \frac{225}{65} + \frac{225}{35} + \frac{225}{65} \] \[\chi^2 \approx 6.43 + 3.46 + 6.43 + 3.46 = 19.78\]

Perhitungan dengan software R

# Buat matriks data: baris = kebiasaan berjalan kaki (Sering, Jarang), 
# kolom = kondisi hipertensi (Ya, Tidak)

data_matrix <- matrix(c(20, 80, 50, 50), nrow = 2, byrow = TRUE)
colnames(data_matrix) <- c("Ya", "Tidak")
rownames(data_matrix) <- c("Sering", "Jarang")

# Uji Chi-Square
chi_test <- chisq.test(data_matrix, correct = FALSE)

# Hasil
list(
  Chi_Square = chi_test$statistic,
  P_Value = chi_test$p.value,
  Decision = ifelse(chi_test$p.value < 0.05, "Tolak H0 (Ada hubungan)", "Gagal Tolak H0 (Tidak ada hubungan)")
)

## $Chi_Square
## X-squared 
##  19.78022 
## 
## $P_Value
## [1] 8.687712e-06
## 
## $Decision
## [1] "Tolak H0 (Ada hubungan)"

Wilayah Kritis:

\[df = (2-1)(2-1) = 1\] \[X^2 tabel(α = 0.05, df = 1) = 3.841 \]

Kriteria Uji: Tolak H₀ apabila χ²Hitung > χ²tabel. Terima untuk selainnya.

Keputusan : Karena χ²hitung = 19.78 > 3.841, maka tolak H₀.

Kesimpulan : Ada hubungan yang signifikan antara kebiasaan berjalan kaki dan hipertensi.

7. Tabel Kontingensi Tiga Arah

Tabel kontingensi 3 arah digunakan untuk menganalisis hubungan antara tiga variabel kategori secara bersamaan. Tabel ini memperluas tabel kontingensi dua arah dengan menambahkan satu dimensi kategori tambahan.

Contoh Kasus
Hipertensi atau tekanan darah tinggi merupakan salah satu penyakit kronis yang banyak ditemukan pada orang dewasa. Gaya hidup, khususnya aktivitas fisik seperti kebiasaan berjalan kaki, diduga memiliki pengaruh terhadap kondisi hipertensi. Selain itu, faktor usia juga dikenal sebagai salah satu determinan penting dalam kejadian hipertensi.

Dalam analisis ini, kita ingin mengetahui apakah terdapat hubungan antara kebiasaan berjalan dan hipertensi, serta bagaimana pengaruh faktor umur terhadap hubungan tersebut. Oleh karena itu, digunakan pendekatan tabel kontingensi tiga arah untuk mengamati hubungan antara tiga variabel kategori:

• Kebiasaan berjalan (rutin / tidak rutin),

• Status hipertensi (hipertensi / tidak hipertensi), dan

• Kelompok umur (dewasa muda / dewasa tua).

library(knitr)
library(gmodels)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

# Membuat data array
data_array <- array(
  data = c(5, 35, 12, 28,   # Dewasa Muda: Rutin, Tidak Rutin
           15, 25, 30, 10), # Dewasa Tua: Rutin, Tidak Rutin
  dim = c(2, 2, 2),
  dimnames = list(
    "Kebiasaan" = c("Rutin", "Tidak Rutin"),
    "Hipertensi" = c("Hipertensi", "Tidak Hipertensi"),
    "Umur" = c("Dewasa Muda", "Dewasa Tua")
  )
)



marginal_total <- addmargins(apply(data_array, c(1,2), sum))
kable(marginal_total, "html", caption = "Tabel Marginal: Kebiasaan × Hipertensi (Gabungan Semua Umur)") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "bordered"), full_width = F)

Tabel Marginal: Kebiasaan × Hipertensi (Gabungan Semua Umur)

	Hipertensi	Tidak Hipertensi	Sum
Rutin	20	42	62
Tidak Rutin	60	38	98
Sum	80	80	160

parsial_muda <- addmargins(data_array[,, "Dewasa Muda"])
kable(parsial_muda, "html", caption = "Tabel Parsial: Dewasa Muda (18–40 Tahun)") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "bordered", "condensed"), full_width = F, position = "center")

Tabel Parsial: Dewasa Muda (18–40 Tahun)

	Hipertensi	Tidak Hipertensi	Sum
Rutin	5	12	17
Tidak Rutin	35	28	63
Sum	40	40	80

parsial_tua <- addmargins(data_array[,, "Dewasa Tua"])
kable(parsial_tua, "html", caption = "Tabel Parsial: Dewasa Tua (>40 Tahun)") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "bordered", "condensed"), full_width = F, position = "center")

Tabel Parsial: Dewasa Tua (>40 Tahun)

	Hipertensi	Tidak Hipertensi	Sum
Rutin	15	30	45
Tidak Rutin	25	10	35
Sum	40	40	80

7.1 Perhitungan Asosiasi

Ukuran asosiasi digunakan untuk mengukur kuat dan arah hubungan antara dua variabel kategori. Dalam konteks tabel kontingensi 2x2, ukuran asosiasi paling umum adalah:

• Risk Ratio (RR) atau Relative Risk

• Risk Difference (RD)

• Odds Ratio (OR)

Untuk tabel kontingensi 3 arah (misalnya A × B × C), maka kita bisa menghitung ukuran asosiasi parsial yakni antara A dan B dengan mengendalikan (mengstratifikasi) variabel C dan juga ukuran asosiasi marginal.

hitung_asosiasi <- function(mat) {
  a <- mat[1,1] # Rutin & Hipertensi
  b <- mat[1,2] # Rutin & Tidak Hipertensi
  c <- mat[2,1] # Tidak Rutin & Hipertensi
  d <- mat[2,2] # Tidak Rutin & Tidak Hipertensi
  
  rr <- (a / (a + b)) / (c / (c + d))
  rd <- (a / (a + b)) - (c / (c + d))
  or <- (a * d) / (b * c)
  
  data.frame(
    `Risk Ratio (RR)` = round(rr, 2),
    `Risk Difference (RD)` = round(rd, 2),
    `Odds Ratio (OR)` = round(or, 2)
  )
}

Ukuran Asosiasi Parsial

asosiasi_muda <- hitung_asosiasi(data_array[,, "Dewasa Muda"])
kable(asosiasi_muda, caption = "Ukuran Asosiasi - Dewasa Muda") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "bordered"))

Ukuran Asosiasi - Dewasa Muda

Risk.Ratio..RR.	Risk.Difference..RD.	Odds.Ratio..OR.
0.53	-0.26	0.33

asosiasi_tua <- hitung_asosiasi(data_array[,, "Dewasa Tua"])
kable(asosiasi_tua, caption = "Ukuran Asosiasi - Dewasa Tua") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "bordered"))

Ukuran Asosiasi - Dewasa Tua

Risk.Ratio..RR.	Risk.Difference..RD.	Odds.Ratio..OR.
0.47	-0.38	0.2

Interpretasi:

Dewasa Muda (18-40 Tahun)

• Risk Ratio (RR) = 0.53
  Peluang hipertensi pada kelompok dengan kebiasaan rutin adalah 0.53 kali peluang hipertensi pada kelompok yang tidak rutin. Ini menunjukkan bahwa kebiasaan rutin berhubungan dengan penurunan risiko hipertensi sebesar 47% pada kelompok dewasa muda.

• Risk Difference (RD) = -0.26
  Selisih risiko hipertensi antara kelompok dengan kebiasaan rutin dan tidak rutin adalah -26%. Artinya, terdapat penurunan absolut kejadian hipertensi sebesar 26% pada kelompok dengan kebiasaan rutin dibandingkan yang tidak.

• Odds Ratio (OR) = 0.33
  Peluang hipertensi pada kelompok dengan kebiasaan rutin adalah seper tiga dari peluang hipertensi pada kelompok yang tidak rutin, menandakan hubungan protektif yang cukup kuat di antara keduanya.

Dewasa Tua (>40 Tahun)

• Risk Ratio (RR) = 0.47
  Pada kelompok dewasa tua, kebiasaan rutin dikaitkan dengan penurunan risiko hipertensi hingga 53% dibandingkan kelompok yang tidak rutin.

• Risk Difference (RD) = -0.38
  Penurunan absolut risiko hipertensi pada kelompok kebiasaan rutin adalah sebesar 38%.

• Odds Ratio (OR) = 0.20
  Odds hipertensi pada kelompok dengan kebiasaan rutin hanya 20% dari odds pada kelompok tidak rutin, menunjukkan hubungan protektif yang lebih kuat daripada pada kelompok dewasa muda.

Kesimpulan Umum

Pada kedua kelompok umur, kebiasaan rutin berperan sebagai faktor protektif terhadap hipertensi. Efek protektif ini tampak sedikit lebih kuat pada kelompok dewasa tua dibandingkan dewasa muda, baik secara relatif (RR dan OR) maupun absolut (RD). Hal ini menunjukkan pentingnya menjaga kebiasaan sehat khususnya pada kelompok yang lebih tua sebagai upaya pencegahan hipertensi.

Ukuran Asosiasi Marginal

marginal_matrix <- apply(data_array, c(1,2), sum)
asosiasi_marginal <- hitung_asosiasi(marginal_matrix)

kable(asosiasi_marginal, caption = "Ukuran Asosiasi - Marginal (Gabungan Semua Umur)") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "bordered"))

Ukuran Asosiasi - Marginal (Gabungan Semua Umur)

Risk.Ratio..RR.	Risk.Difference..RD.	Odds.Ratio..OR.
0.53	-0.29	0.3

Interpretasi Ukuran Asosiasi Marginal (Gabungan Semua Umur)

• Risk Ratio (RR) = 0.53
  Peluang hipertensi pada kelompok dengan kebiasaan rutin adalah 0.53 kali peluang hipertensi pada kelompok yang tidak rutin. Ini menunjukkan kebiasaan rutin terkait dengan penurunan risiko hipertensi sebesar 47% secara keseluruhan, tanpa membedakan kelompok umur.

• Risk Difference (RD) = -0.29
  Perbedaan risiko absolut hipertensi antara kelompok kebiasaan rutin dan tidak rutin adalah sebesar 29%, yang berarti penurunan kejadian hipertensi absolut sebesar 29% pada mereka yang rutin dibandingkan yang tidak.

• Odds Ratio (OR) = 0.30
  Peluang hipertensi pada kelompok dengan kebiasaan rutin sekitar 30% dari peluang pada kelompok yang tidak rutin, menunjukkan hubungan protektif yang kuat secara umum di seluruh kelompok umur.

Kesimpulan

Secara keseluruhan, kebiasaan rutin memiliki efek protektif yang signifikan terhadap hipertensi pada populasi gabungan tanpa memandang kelompok umur. Ini memperkuat pentingnya kebiasaan rutin dalam pencegahan hipertensi di berbagai kelompok umur.

7.2 Pengujian Conditional Independence

Dalam konteks tabel kontingensi 3 arah (misalnya variabel A, B, dan C), conditional independence menguji apakah dua variabel A dan B independen satu sama lain setelah dikontrol oleh variabel ketiga C.

Kita ingin menguji apakah Kebiasaan Berjalan dan Hipertensi independen, setelah dikendalikan oleh umur.

Hipotesis:

• H₀: Tidak terdapat hubungan antara kebiasaan berjalan dan hipertensi, setelah dikontrol oleh umur tertentu.

• H₁: Terdapat hubungan antara kebiasaan berjalan dan hipertensi, setelah dikontrol oleh umur tertentu.

Taraf Signifikansi: α = 0.05

Statistik Uji: Chi-Square

Perhitugan dengan Software R

# Uji Chi-Square untuk masing-masing strata usia

chisq.test(data_array[,, "Dewasa Muda"]) # untuk kontrol umur = dewasa muda

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data_array[, , "Dewasa Muda"]
## X-squared = 2.6891, df = 1, p-value = 0.101

chisq.test(data_array[,, "Dewasa Tua"]) # untuk kontrol umur =dewasa tua

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data_array[, , "Dewasa Tua"]
## X-squared = 9.9556, df = 1, p-value = 0.001604

Kriteria Uji : Tolak H₀ apabila p-value < 0.05. Terima selainnya.

Keputusan :

1. P-value pada uji conditional independence pada strata umur dewasa muda adalah sebesar 0.101 > 0.05. sehingga terima H₀.
2. P-value pada uji conditional independence pada strata umur dewasa tua adalah sebesar 0.001604 > 0.05. sehingga tolak H₀.

Kesimpulan: Terdapat hubungan antara kebiasaan berjalan dan hipertensi atau hubungan keduanya tidaklah kondisional independen., setelah dikontrol oleh umur tertentu.

7.3 Inferensia Tabel Kontingensi Tiga Arah

Dalam analisis data kategorik dua variabel dengan satu variabel pengganggu (strata), penting untuk memahami bagaimana hubungan antar variabel utama dapat berbeda pada tiap strata. Oleh karena itu, langkah pertama yang dilakukan adalah memisahkan data ke dalam tabel kontingensi 2×2 berdasarkan strata, seperti kelompok usia atau jenis kelamin. Setelah itu, dilakukan Uji Homogenitas Odds Ratio menggunakan Statistik Breslow-Day untuk menguji apakah odds ratio (OR) antar strata bersifat homogen (sama) atau heterogen (berbeda).

Jika odds ratio antar strata terbukti homogen, maka analisis dapat dilanjutkan dengan Uji Mantel-Haenszel, yang digunakan untuk menguji hubungan antara dua variabel utama (misalnya, kebiasaan olahraga dan hipertensi) sambil mengendalikan efek variabel pengganggu (strata). Uji ini memberikan estimasi gabungan odds ratio yang lebih akurat dan sah karena telah mempertimbangkan struktur stratifikasi dalam data. Namun, apabila hasil uji Breslow-Day menunjukkan bahwa odds ratio antar strata berbeda secara signifikan (tidak homogen), maka penggunaan uji Mantel-Haenszel tidak disarankan karena penggabungan data antar strata bisa menyesatkan dan menghasilkan kesimpulan yang tidak valid.

Uji Homogenitas Odds Ratio (Menggunakan Statistik Breslow-Day)

Uji ini digunakan untuk mengevaluasi apakah hubungan antara dua variabel (misalnya: kebiasaan dan hipertensi) konsisten di seluruh strata dari variabel ketiga (misalnya: kelompok umur). Uji ini penting untuk memastikan asumsi homogenitas odds ratio sebelum menggunakan uji Cochran-Mantel-Haenszel (CMH).

Langkah-langkah Pengerjaan:

1. Estimasi Rasio Odds Gabungan (Mantel-Haenszel OR)

   Langkah pertama adalah menghitung estimasi odds ratio gabungan (\(\widehat{OR}_{MH}\)) dari seluruh strata dengan rumus: \[ \widehat{OR}_{MH} = \frac{\sum_{j=1}^{K} \frac{a_j d_j}{n_j}}{\sum_{j=1}^{K} \frac{b_j c_j}{n_j}} \] Keterangan:

• \(a_j\): jumlah kasus terpapar di strata ke-\(j\)

• \(b_j\): jumlah kontrol terpapar di strata ke-\(j\)

• \(c_j\): jumlah kasus tidak terpapar di strata ke-\(j\)

• \(d_j\): jumlah kontrol tidak terpapar di strata ke-\(j\)

• \(n_j = a_j + b_j + c_j + d_j\): total pengamatan di strata ke-\(j\)

• \(K\): jumlah strata

2. Menentukan Nilai Ekspektasi \(\tilde{a}_j\)

Setelah \(\widehat{OR}_{MH}\) diperoleh, kita menghitung nilai ekspektasi untuk sel \(a\) dalam setiap strata (\(\tilde{a}_j\)) berdasarkan asumsi odds ratio gabungan tersebut benar. Caranya dengan menyelesaikan persamaan kuadrat berikut: \[-m_{1j}n_{1j}\widehat{OR}_{MH} + (n_{2j} – m_{1j} + \widehat{OR}_{MH}(n_{1j} + m_{1j}))\tilde{a}_j + (1 – \widehat{OR}_{MH})\tilde{a}_j^2 = \] Keterangan:

• \(m_{1j}\): total subjek terpapar pada strata ke-\(j\) \((a_j + b_j)\)

• \(m_{2j}\): total subjek tidak terpapar pada strata ke-\(j\) \((c_j + d_j)\)

•\(n_{1j}\): total kasus pada strata ke-\(j\) \((a_j + c_j)\)

•\(n_{2j}\): total non-kasus pada strata ke-\(j\) \((b_j + d_j)\)

Dari duasolusi kuadrat, hanya solusi positif yang berada dalam batas realistis data, yaitu 0 < \(\tilde{a}_j\) ≤ min (\(n_{1j}\), \(m_{1j}\))

Nilai \(\tilde{a}_j\) ini mewakili jumlah kasus terpapar yang diharapkan jika odds ratio homogen benar.

3. Menghitung Varians Kondisional Setelah mendapatkan nilai ekspektasi \(\tilde{a}_j\), hitung juga varians dari nilai tersebut: \[Var(a_j|\widehat{OR}_{MH}) = \left(\frac{1}{\tilde{a}_j} + \frac{1}{\tilde{b}_j} + \frac{1}{\tilde{c}_j} + \frac{1}{\tilde{d}_j}\right)^{-1}\] dengan:

• \(\tilde{b}j = m{1j} - \tilde{a}_j\)

• \(\tilde{c}j = n{1j} - \tilde{a}_j\)

• \(\tilde{d}j = m{2j} - \tilde{c}_j\)

  Varians ini merefleksikan seberapa besar penyimpangan yang bisa diterima dari nilai ekspektasi \(\tilde{a}_j\).

4. Menghitung Statistik Uji Breslow-Day

   Gunakan rumus berikut untuk menghitung nilai statistik uji utama: \[X_{HBD}^{2} = \sum_{j=1}^{K} \frac{(a_{j} - \tilde{a}_{j})^{2}}{Var(a_{j}|\widehat{OR}_{MH})}\] Semakin besar nilai ini, semakin besar ketidaksesuaian antara data yang diamati dan ekspektasi jika odds ratio homogen.

5. Koreksi Tarone Karena statistik Breslow-Day bisa sedikit bias pada ukuran sampel kecil, maka koreksi Tarone bisa diterapkan: \[X_{HBDT}^2 = X_{HBD}^2 - \frac{(\sum_{j=1}^{K}(a_j - \tilde{a}_j))^2}{\sum_{j=1}^{K}Var(a_j|\widehat{OR}_{MH})}\]

   Koreksi ini mengurangi potensi kesalahan Type I, terutama saat penyimpangan kecil secara sistematik muncul di seluruh strata.

6. Distribusi Statistik Baik \(X^2_{HBD}\) maupun \(X^2_{HBDT}\) mengikuti distribusi chi-square dengan derajat bebas sebanyak \((K - 1)\).

---

Penerapan Dalam Contoh Kasus

Data

Tabel Parsial Dewasa Muda (18-40 Tahun)

Kebiasaan	Hipertensi	Tidak Hipertensi	Jumlah
Rutin	5	12	17
Tidak Rutin	35	28	63
Jumlah	40	40	80

Tabel Parsial Dewasa Tua (> 40 Tahun)

Kebiasaan	Hipertensi	Tidak Hipertensi	Jumlah
Rutin	15	30	45
Tidak Rutin	25	10	35
Jumlah	40	40	80

Hipotesis

• \(H_0: OR_1 = OR_2 = \dots = OR_K\)
  (Odds ratio pada semua strata adalah sama / homogen)

• \(H_1: \exists j,k \text{ sehingga } OR_j \neq OR_k\)
  (Ada minimal dua strata yang memiliki odds ratio berbeda / tidak homogen)

Taraf signifikansi: \(\alpha = 0.05\) (5%)

Statistik Uji 1. Perhitungan Estimasi Rasio Odds Gabungan (Mantel-Haenszel OR)

Strata 1 (Dewasa Muda):

• \(a_1 = 5\), \(b_1 = 12\), \(c_1 = 35\), \(d_1 = 28\)
• \(n_1 = 80\)

Hitung:

\[ \frac{a_1 d_1}{n_1} = \frac{5 \times 28}{80} = \frac{140}{80} = 1.75 \]

\[ \frac{b_1 c_1}{n_1} = \frac{12 \times 35}{80} = \frac{420}{80} = 5.25 \]

Strata 2 (Dewasa Tua):

• \(a_2 = 15\), \(b_2 = 30\), \(c_2 = 25\), \(d_2 = 10\)
• \(n_2 = 80\)

Hitung:

\[ \frac{a_2 d_2}{n_2} = \frac{15 \times 10}{80} = \frac{150}{80} = 1.875 \]

\[ \frac{b_2 c_2}{n_2} = \frac{30 \times 25}{80} = \frac{750}{80} = 9.375 \]

Rasio Odds Gabungan (Mantel-Haenszel OR):

\[ \widehat{OR}_{MH} = \frac{1.75 + 1.875}{5.25 + 9.375} = \frac{3.625}{14.625} \approx 0.248 \]

2. Menghitung Nilai Ekspektasi \(\tilde{a}_j\)

Gunakan persamaan kuadrat:

\[ -(m_{1j} n_{1j}) \widehat{OR}_{MH} + (n_{2j} - m_{1j} + \widehat{OR}_{MH}(n_{1j} + m_{1j})) \tilde{a}_j + (1 - \widehat{OR}_{MH}) \tilde{a}_j^2 = 0 \]

Untuk Strata 1:

• \(m_{1j} = 17\), \(m_{2j} = 63\)
• \(n_{1j} = 40\), \(n_{2j} = 40\)
• \(\widehat{OR}_{MH} = 0.248\)

Substitusi:

\[ -(17)(40)(0.248) + (40 - 17 + 0.248(40 + 17)) \tilde{a}_1 + (1 - 0.248) \tilde{a}_1^2 = 0 \]

Hitung koefisien:

• \(A = 1 - 0.248 = 0.752\)
• \(B = 23 + 0.248 \times 57 = 23 + 14.136 = 37.136\)
• \(C = -17 \times 40 \times 0.248 = -168.32\)

Persamaan:

\[ 0.752 \tilde{a}_1^2 + 37.136 \tilde{a}_1 - 168.32 = 0 \]

Gunakan rumus kuadrat:

\[ \tilde{a}_1 = \frac{-B \pm \sqrt{B^2 - 4AC}}{2A} \]

Hitung diskriminan:

\[ \sqrt{37.136^2 + 4 \times 0.752 \times 168.32} \approx \sqrt{1379.1 + 506.3} = \sqrt{1885.4} \approx 43.43 \]

Solusi positif:

\[ \tilde{a}_1 = \frac{-37.136 + 43.43}{2 \times 0.752} = \frac{6.294}{1.504} \approx 4.19 \]

Untuk Strata 2:

• \(m_{1j} = 45\), \(m_{2j} = 35\)
• \(n_{1j} = 40\), \(n_{2j} = 40\)

Persamaan:

\[ -(45)(40)(0.248) + (40 - 45 + 0.248 (40 + 45)) \tilde{a}_2 + (1 - 0.248) \tilde{a}_2^2 = 0 \]

Hitung:

• \(C = -448.8\)
• \(B = -5 + 0.248 \times 85 = -5 + 21.08 = 16.08\)
• \(A = 0.752\)

Persamaan kuadrat:

\[ 0.752 \tilde{a}_2^2 + 16.08 \tilde{a}_2 - 448.8 = 0 \]

Diskriminan:

\[ \sqrt{16.08^2 + 4 \times 0.752 \times 448.8} = \sqrt{258.6 + 1350.2} = \sqrt{1608.8} \approx 40.1 \]

Solusi positif:

\[ \tilde{a}_2 = \frac{-16.08 + 40.1}{2 \times 0.752} = \frac{24.02}{1.504} \approx 15.97 \]

3. Hitung Varians \(Var(a_j | \widehat{OR}_{MH})\)

Rumus:

\[ Var(a_j | \widehat{OR}_{MH}) = \left( \frac{1}{\tilde{a}_j} + \frac{1}{\tilde{b}_j} + \frac{1}{\tilde{c}_j} + \frac{1}{\tilde{d}_j} \right)^{-1} \]

Untuk Strata 1:

• \(\tilde{a}_1 = 4.19\)
• \(\tilde{b}_1 = 17 - 4.19 = 12.81\)
• \(\tilde{c}_1 = 40 - 4.19 = 35.81\)
• \(\tilde{d}_1 = 63 - 35.81 = 27.19\)

Hitung varians:

\[ Var_1 = \left( \frac{1}{4.19} + \frac{1}{12.81} + \frac{1}{35.81} + \frac{1}{27.19} \right)^{-1} \approx (0.239 + 0.078 + 0.028 + 0.037)^{-1} = 0.382^{-1} \approx 2.62 \]

Untuk Strata 2:

• \(\tilde{a}_2 = 15.97\)
• \(\tilde{b}_2 = 45 - 15.97 = 29.03\)
• \(\tilde{c}_2 = 40 - 15.97 = 24.03\)
• \(\tilde{d}_2 = 35 - 24.03 = 10.97\)

Hitung varians:

\[ Var_2 = \left( \frac{1}{15.97} + \frac{1}{29.03} + \frac{1}{24.03} + \frac{1}{10.97} \right)^{-1} \approx (0.063 + 0.034 + 0.042 + 0.091)^{-1} = 0.23^{-1} \approx 4.35 \]

4. Hitung Statistik Uji Breslow-Day

Statistik:

\[ X_{BD}^2 = \sum \frac{(a_j - \tilde{a}_j)^2}{Var(a_j | \widehat{OR}_{MH})} \]

Hitung:

\[ X_{BD}^2 = \frac{(5 - 4.19)^2}{2.62} + \frac{(15 - 15.97)^2}{4.35} = \frac{0.6561}{2.62} + \frac{0.9409}{4.35} \approx 0.25 + 0.216 = 0.466 \]

Titik Kritis dan Keputusan

• Chi-square dengan derajat kebebasan \(df = K - 1 = 2 - 1 = 1\): 3.841
• Aturan uji: Tolak \(H_0\) jika \(X^2_{BD} > 3.841\) atau p-value < 0.05

Karena \(X^2_{BD} = 0.466 < 3.841\), maka gagal menolak \(H_0\).

Interpretasi Hasil

Tidak terdapat bukti yang cukup pada tingkat signifikansi 5% untuk menolak hipotesis bahwa odds ratio antar strata umur homogen. Dengan kata lain, odds ratio dapat dianggap sama di kedua strata umur.

Perhitungan dengan Software R

library(DescTools)  

## Registered S3 method overwritten by 'DescTools':
##   method         from 
##   reorder.factor gdata

# Misal data_array adalah array data 3 dimensi untuk strata
BreslowDayTest(data_array)

## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  data_array
## X-squared = 0.44566, df = 1, p-value = 0.5044

7.4 Pengujian Statistik untuk Independensi Bersyarat

Uji Cochran-Mantel

Tujuan Uji CMH

• Menguji apakah terdapat hubungan antara dua variabel kategori setelah mengendalikan variabel perancu.
  
• Mengatasi bias yang disebabkan oleh variabel pengganggu (confounding).
  
• Menggabungkan hasil tabel-tabel kontingensi berlapis menjadi satu statistik uji.

---

Konsep Dasar

• Data terdiri atas beberapa tabel kontingensi 2×2 yang mewakili strata berbeda berdasarkan variabel pengganggu.
  
• Hipotesis null menyatakan bahwa hubungan antara dua variabel tidak ada (independen) di setiap strata.
  
• Hipotesis alternatif menyatakan ada hubungan pada sedikitnya satu strata.
  
• Uji ini menggabungkan statistik dari masing-masing strata untuk menghasilkan satu hasil uji keseluruhan.

---

Hipotesis

• \(H_0\): Tidak ada asosiasi antara dua variabel kategori dalam setiap strata, atau odds ratio strata sama dengan 1.

• \(H_1\): Minimal ada satu strata dengan asosiasi antara dua variabel kategori, odds ratio strata tidak sama dengan 1.

---

Statistik Uji CMH

Statistik CMH dirumuskan sebagai

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

dengan

• \(n_{11k}\): frekuensi sel baris 1 kolom 1 pada strata \(k\),
  
• \(\mu_{11k} = \frac{n_{1.k} \times n_{.1k}}{n_{..k}}\): ekspektasi nilai \(n_{11k}\) pada strata \(k\),
  
• \(\operatorname{var}(n_{11k}) = \frac{n_{1.k} \times n_{2.k} \times n_{.1k} \times n_{.2k}}{n_{..k}^2 (n_{..k} - 1)}\): varians nilai \(n_{11k}\),
  
• \(n_{1.k}\), \(n_{2.k}\), \(n_{.1k}\), \(n_{.2k}\), dan \(n_{..k}\): marginal dan total pada strata \(k\).

---

Distribusi dan Keputusan

• Statistik CMH mengikuti distribusi Chi-square (\(\chi^2\)) dengan derajat kebebasan (df) = 1.
  
• Tolak \(H_0\) jika statistik CMH lebih besar dari nilai kritis \(\chi^2_{1,\alpha}\) pada tingkat signifikansi \(\alpha\), atau jika p-value < \(\alpha\).

---

Kelebihan Uji CMH

• Mengontrol sekaligus memperhitungkan efek variabel perancu tanpa memerlukan model yang kompleks.
  
• Lebih valid dibandingkan uji independensi sederhana bila ada variabel pengganggu.
  
• Dapat diaplikasikan pada berbagai bidang, misal epidemiologi, ilmu sosial, dan eksperimen.

---

Langkah-langkah Perhitungan Manual

---

1. Tentukan Nilai \(n_{11k}\), Marginal Baris dan Kolom, dan Total

Strata	\(n _ {11k}\)	\(n _ {1.k}\)	\(n _ {.1k}\)	\(n _ {2.k}\)	\(n _ {.2k}\)	\(n _ {..k}\)
Dewasa Muda	5	17	40	63	40	80
Dewasa Tua	15	45	40	35	40	80

---

2. Hitung Ekspektasi \(\mu_{11k}\) untuk tiap strata

\[ \mu_{11k} = \frac{n_{1.k} \times n_{.1k}}{n_{..k}} \]

\[ \begin{cases} \mu_{11,1} = \frac{17 \times 40}{80} = 8.5\\[6pt] \mu_{11,2} = \frac{45 \times 40}{80} = 22.5 \end{cases} \]

---

3. Hitung Varians \(\operatorname{var}(n_{11k})\)

\[ \operatorname{var}(n_{11k}) = \frac{n_{1.k} \times n_{2.k} \times n_{.1k} \times n_{.2k}}{n_{..k}^2 (n_{..k} - 1)} \]

\[ \begin{cases} \operatorname{var}(n_{11,1}) = \frac{17 \times 63 \times 40 \times 40}{80^2 \times 79} \approx 3.40 \\[6pt] \operatorname{var}(n_{11,2}) = \frac{45 \times 35 \times 40 \times 40}{80^2 \times 79} \approx 4.99 \end{cases} \]

---

4. Hitung nilai \((n_{11k} - \mu_{11k})\) dan kuadratkan jumlahnya

\[ \begin{cases} 5 - 8.5 = -3.5 \\[6pt] 15 - 22.5 = -7.5 \end{cases} \]

Jumlah:

\[ -3.5 + (-7.5) = -11 \]

Kuadratkan:

\[ (-11)^2 = 121 \]

---

5. Hitung jumlah varians

\[ 3.40 + 4.99 = 8.39 \]

---

6. Hitung Statistik CMH

\[ CMH = \frac{121}{8.39} \approx 14.42 \]

---

7. Keputusan Uji

• Nilai kritis \(\chi^2_{1,0.05} = 3.841\)
  
• Karena \(14.42 > 3.841\), tolak hipotesis nol
  
• Artinya, ada hubungan yang signifikan antara kebiasaan rutin dan hipertensi setelah mengontrol umur.

---

Kesimpulan

Hasil uji CMH menunjukkan ada hubungan signifikan antara kebiasaan rutin dan hipertensi meskipun setelah mengontrol strata umur. Oleh karena itu, kebiasaan rutin berpengaruh terhadap kejadian hipertensi di kelompok umur yang berbeda.

---

Perhitungan dengan Software R

library(vcdExtra)  

## Loading required package: vcd

## Loading required package: grid

## Loading required package: gnm

## 
## Attaching package: 'vcdExtra'

## The following object is masked from 'package:dplyr':
## 
##     summarise

# Jalankan uji CMH  
 cmh_base <- mantelhaen.test(data_array, correct = FALSE)
 cmh_base

## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  data_array
## Mantel-Haenszel X-squared = 14.45, df = 1, p-value = 0.0001439
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.1184442 0.5186932
## sample estimates:
## common odds ratio 
##         0.2478632

8. Generalized Linear Model (GLM)

8.1 Pendahuluan

Generalized Linear Model (GLM) adalah perluasan dari model regresi linear klasik yang memungkinkan penanganan berbagai jenis data dan hubungan yang lebih kompleks antara variabel. Model regresi linear tradisional mengasumsikan dua hal penting: variabel respons (dependent variable) terdistribusi normal, dan hubungan antara variabel prediktor dan respons bersifat linear secara langsung. Namun, dalam banyak kasus nyata, asumsi-asumsi ini tidak terpenuhi.

Sebagai contoh:

Variabel respons bisa berupa jumlah kejadian (data count), yang biasanya tidak berdistribusi normal, melainkan berdistribusi Poisson. Variabel respons bisa berupa kategori (kualitatif), misalnya keberhasilan atau kegagalan, yang mengikuti distribusi binomial. Hubungan antara variabel bebas dan variabel respons terkadang tidak linear secara langsung, melainkan perlu transformasi tertentu agar model dapat merepresentasikannya dengan baik. Untuk mengatasi keterbatasan tersebut, GLM menyediakan framework yang fleksibel dan kuat untuk membangun model statistik yang sesuai dengan berbagai jenis data dan pola hubungan. GLM menyatukan konsep statistik dasar seperti distribusi probabilitas dari keluarga eksponensial, fungsi link, dan linear predictor dalam satu kerangka yang koheren.

Melalui GLM, kita dapat:

Memodelkan variabel respons yang beragam seperti data biner, count, proporsi, dan data kontinu yang tidak normal. Menggunakan fungsi link untuk menghubungkan mean dari variabel respons dengan kombinasi linier dari prediktor, sehingga memungkinkan pemodelan hubungan nonlinier secara langsung dalam domain parameter. Melakukan estimasi parameter dan inferensi statistik dengan pendekatan maximum likelihood dan metode lain yang relevan. Penggunaan GLM sangat luas, mulai dari bidang epidemiologi, biostatistik, ilmu sosial, ekonomi, sampai ilmu komputer. Model ini memungkinkan peneliti dan praktisi melakukan analisis yang tepat terhadap data yang kompleks dan heterogen.

8.2 Jenis-jenis Generalized Linear Models (GLM)

Generalized Linear Model (GLM) adalah framework umum yang mencakup berbagai model regresi yang sesuai dengan tipe data dan distribusi berbeda. Pemilihan jenis GLM bergantung pada sifat variabel respons dan hubungan yang ingin dimodelkan.

Berikut adalah jenis-jenis GLM yang paling umum digunakan:

1. Regresi Linear Klasik

• Variabel Respons: Kontinu, biasanya diasumsikan berdistribusi Normal.
  
• Fungsi Link: Identitas (identity), artinya hubungan langsung linear:
  \[ \mu = \eta = X\beta \]
  
• Penggunaan: Model regresi linear biasa.
  
• Contoh Kasus: Prediksi tinggi badan, berat badan, suhu.

2. Regresi Logistik (Logistic Regression)

• Variabel Respons: Biner (dua kategori), mengikuti distribusi Binomial.
  
• Fungsi Link: Logit, menghubungkan probabilitas \(p\) dengan prediktor linear:
  \[ \text{logit}(p) = \log \left(\frac{p}{1-p}\right) = X\beta \]
  
• Penggunaan: Mengestimasi probabilitas kejadian suatu peristiwa.
  
• Contoh Kasus: Diagnosa penyakit (sakit/sehat), keberhasilan percobaan (ya/tidak).

3. Regresi Binomial (Binomial Regression)

• Variabel Respons: Proporsi atau jumlah sukses dalam sejumlah percobaan (binomial).
  
• Fungsi Link: Logit, Probit, atau Complementary log-log (cloglog), sesuai model.
  
• Penggunaan: Situasi dengan data proporsi, misal proporsi pasien sembuh.
  
• Contoh Kasus: Persentase keberhasilan vaksin dari total peserta.

4. Regresi Poisson (Poisson Regression)

• Variabel Respons: Data hitung (count data), mengikuti distribusi Poisson.
  
• Fungsi Link: Logaritma, yaitu:
  \[ \log(\mu) = X\beta \]
  
• Penggunaan: Memodelkan jumlah kejadian dalam interval waktu/ruang.
  
• Contoh Kasus: Jumlah kecelakaan lalu lintas per hari, jumlah panggilan telepon per jam.

5. Regresi Multinomial (Multinomial Logistic Regression)

• Variabel Respons: Kategori lebih dari dua kelas (multikategori tanpa urutan).
  
• Fungsi Link: Logit terhadap referensi kategori.
  
• Penggunaan: Mengestimasi probabilitas keanggotaan dalam beberapa kategori.
  
• Contoh Kasus: Preferensi produk (produk A, B, C), pilihan transportasi.

6. Regresi Ordinal (Ordinal Logistic Regression)

• Variabel Respons: Kategori dengan urutan (ordinal).
  
• Fungsi Link: Logit atau probit kumulatif, mengkaji peluang kumulatif kategori.
  
• Penggunaan: Memodelkan data skor/rating yang berurutan.
  
• Contoh Kasus: Tingkat kepuasan pelanggan (rendah, sedang, tinggi).

7. Regresi Gamma dan Inverse Gaussian

• Variabel Respons: Data kontinu positif, dengan distribusi miring (skewed).
  
• Fungsi Link: Log atau inverse.
  
• Penggunaan: Memodelkan variabel positif seperti waktu bertahan, biaya.
  
• Contoh Kasus: Biaya perawatan rumah sakit, waktu tunggu pelayanan.

8.3 Fungsi Link Umum pada Beberapa GLM

Fungsi link dalam konteks generalized linear model (GLM), termasuk regresi logistik, adalah fungsi yang menghubungkan nilai rata-rata variabel respon dengan kombinasi linear dari prediktor.
Dalam GLM, kita tidak langsung memodelkan nilai rata-rata \(\mu = E(Y)\) sebagai fungsi linear dari variabel prediktor \(X\). Melainkan kita memodelkan transformasi \(\mu\) melalui fungsi link agar hasilnya berbentuk garis linear, yaitu:
\[ g(\mu) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots \]

Model	Distribusi	Fungsi Link	Bentuk Fungsi Link \(g(\mu)\)
Regresi Linier	Normal	Identitas	\(g(\mu) = \mu\)
Regresi Logistik	Binomial	Logit	\(g(p ) = \log \frac{p}{1-p}\)
Regresi Poisson	Poisson	Log	\(g(\mu) = \log \mu\)
Regresi Probit	Binomial	Probit	\(g(p) = \Phi^{-1}(p)\) (invers CDF normal)

8.4 Distribusi Eksponensial

Distribusi exponential family adalah kelas distribusi probabilitas yang memiliki bentuk matematis tertentu. Distribusi ini sangat penting dalam Generalized Linear Models (GLM) karena memungkinkan pemodelan berbagai tipe data dengan pendekatan yang seragam dan kompak.

8.4.1 Peran Exponential Family dalam GLM

Generalized Linear Model dibangun dengan asumsi:

• Variabel respons mengikuti distribusi dari exponential family.
  
• Fungsi link yang menghubungkan rata-rata respons ke prediktor linear berdasarkan parameter natural \(\theta\).
  
• Dengan asumsi ini, langkah estimasi parameter menjadi lebih sistematis dan efisien melalui metode maximum likelihood yang didukung sifat matematis distribusi ini.

Definisi Distribusi Exponential Family

Sebuah variabel acak \(Y\) dikatakan mengikuti distribusi dari exponential family jika fungsi densitas (atau fungsi massa probabilitas) dapat dituliskan dalam bentuk:

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

di mana:
- \(\theta\) adalah parameter natural (natural parameter), berhubungan dengan rata-rata distribusi.
- \(\phi\) adalah parameter dispersi (dispersion parameter).
- \(a(\cdot)\), \(b(\cdot)\), dan \(c(\cdot)\) adalah fungsi yang bergantung pada distribusi tertentu.

Komponen Utama

• Parameter natural \(\theta\): Menghubungkan distribusi dengan mean (ekspektasi) variabel acak.
  
• Fungsi \(b(\theta)\): Fungsi cumulant generating, menentukan mean dan varians distribusi.
  
• Parameter dispersi \(\phi\): Mengontrol varians relatif dari distribusi.

Koneksi ke Mean dan Varians

Mean atau rata-rata \(\mu = E(Y)\) dan varians \(Var(Y)\) dapat dihitung melalui fungsi \(b(\theta)\):

\[ \mu = b'(\theta) \quad \text{dan} \quad Var(Y) = a(\phi) b''(\theta) \]

Dimana \(b'(\theta)\) dan \(b''(\theta)\) adalah turunan pertama dan kedua fungsi \(b(\theta)\).

8.5 Model Regresi Logistik

Regresi logistik adalah metode penting dalam statistik dan machine learning yang digunakan untuk memodelkan hubungan antara satu atau lebih variabel independen dengan variabel respons biner (dua kelas). Berbeda dengan regresi linear yang memprediksi nilai kontinu, regresi logistik memprediksi probabilitas suatu kejadian dengan rentang nilai antara 0 dan 1.

Pengertian

Regresi logistik mengkombinasikan nilai input secara linear dengan koefisien (bobot) untuk menghasilkan prediksi awal. Karena targetnya adalah nilai biner (0 atau 1), fungsi aktivasi sigmoid digunakan untuk mengubah prediksi tersebut menjadi probabilitas yang berada di dalam rentang 0 hingga 1.

Fungsi sigmoid didefinisikan sebagai:

\[ \sigma(z) = \frac{1}{1 + e^{-z}} \]

dengan

\[ z = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p \]

dimana \(x_1, \ldots, x_p\) adalah variabel prediktor dan \(\beta_0, \ldots, \beta_p\) adalah koefisien regresi.

Kegunaan Regresi Logistik

Regresi logistik digunakan untuk klasifikasi biner, di mana variabel target hanya memiliki dua kemungkinan nilai, seperti:

• Sukses/Gagal
  
• Ya/Tidak
  
• Positif/Negatif

Contoh Penerapan Regresi Logistik

• Medis: Memprediksi probabilitas pasien mengalami serangan jantung berdasarkan variabel seperti berat badan, usia, tekanan darah.
  
• Penerimaan Mahasiswa: Memperkirakan kemungkinan diterima di program studi berdasarkan nilai GRE, GMAT, TOEFL.
  
• Spam Email: Mengklasifikasikan email sebagai spam atau bukan dengan menganalisis fitur isi email.

Keunggulan Regresi Logistik

1. Cocok untuk Data Biner: Efektif untuk hasil dengan dua kelas serta data yang dapat dipisahkan secara linear.
   
2. Interpretasi Mudah: Koefisien model berhubungan dengan log-odds kejadian target.
   
3. Komputasi Efisien: Cepat dan ringan dibandingkan metode machine learning lain.

Persamaan Regresi Logistik

Model memetakan input menjadi probabilitas menggunakan fungsi sigmoid:

\[ p = P(Y=1|X) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)}} \]

atau

Fungsi sigmoid dirumuskan sebagai:

\[ f(x) = \frac{1}{1 + e^{-x}} \]

Di mana:

• \(x\) adalah input dari fungsi, berupa kombinasi linear variabel prediktor dan koefisien regresi (\(x = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p\)).
  
• \(f(x)\) menghasilkan nilai antara 0 dan 1, yang merupakan probabilitas.

Persamaan log-odds (logit):

\[ \log \frac{p}{1-p} = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p \]

• Probabilitas \(p\) berada dalam rentang 0 sampai 1.
  
• Ambang batas 0.5 biasanya digunakan untuk klasifikasi:
  • Jika \(p \geq 0.5\), prediksi kelas 1 (positif)
    
  • Jika \(p < 0.5\), prediksi kelas 0 (negatif)

Kesimpulan

Regresi logistik adalah teknik analisis statistik yang efektif dan populer untuk masalah klasifikasi biner. Model ini memberikan probabilitas yang mudah diinterpretasikan dan efisien secara komputasi.

8.5.1 Estimasi Parameter pada Regresi Logistik

Estimasi parameter \(\beta_0, \beta_1, \ldots, \beta_p\) dalam model regresi logistik berbeda dengan regresi linear biasa karena fungsi link yang digunakan adalah fungsi logit (log odds), yang bersifat non-linear.

Metode Maximum Likelihood Estimation (MLE)

Parameter regresi logistik diestimasi menggunakan metode Maximum Likelihood Estimation (MLE).
MLE mencari nilai parameter \(\beta\) yang memaksimalkan fungsi likelihood, yaitu probabilitas mengamati data yang sebenarnya dengan model yang diberikan.

Fungsi Likelihood

Misalkan dataset terdiri dari \(n\) observasi \((x_i, y_i)\) dengan \(y_i \in \{0,1\}\).

Probabilitas prediksi untuk observasi ke-\(i\):

\[ p_i = P(Y=1|x_i) = \frac{e^{\eta_i}}{1 + e^{\eta_i}} \quad \text{dengan} \quad \eta_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip} \]

Likelihood data \(L(\beta)\) adalah:

\[ L(\beta) = \prod_{i=1}^n p_i^{y_i} (1 - p_i)^{1 - y_i} \]

Fungsi Log-Likelihood

Untuk mempermudah optimasi, digunakan fungsi log-likelihood:

\[ \ell(\beta) = \log L(\beta) = \sum_{i=1}^n \Big( y_i \log p_i + (1 - y_i) \log (1 - p_i) \Big) \]

MLE bertujuan untuk menemukan \(\hat{\beta}\) yang memaksimalkan \(\ell(\beta)\).

Cara Estimasi

• Fungsi log-likelihood tidak dapat diselesaikan secara analitik sehingga digunakan metode numerik (iteratif) seperti Newton-Raphson atau Fisher Scoring.
  
• Proses ini ada di balik layar fungsi glm() di R dengan family binomial.

Interpretasi

• Parameter \(\hat{\beta}_j\) mengukur perubahan log odds akibat perubahan satu unit pada variabel \(x_j\), dengan variabel lain tetap.
  
• Nilai parameter yang signifikan menunjukan variabel tersebut berpengaruh terhadap probabilitas kejadian.

Contoh sederhana estimasi parameter dengan R

# Data simulasi
set.seed(123)
n <- 100
x <- rnorm(n)
# Probabilitas logit dengan beta0= -1, beta1=2
p <- 1 / (1 + exp(-(-1 + 2 * x)))
y <- rbinom(n, 1, p)

# Estimasi model regresi logistik
model <- glm(y ~ x, family = binomial)

# Melihat hasil estimasi
summary(model)

## 
## Call:
## glm(formula = y ~ x, family = binomial)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.9961     0.2887  -3.451 0.000559 ***
## x             2.0262     0.4205   4.819 1.44e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 131.79  on 99  degrees of freedom
## Residual deviance:  90.54  on 98  degrees of freedom
## AIC: 94.54
## 
## Number of Fisher Scoring iterations: 5

---

8.5.2 Contoh Kasus

Prediksi Jenis Transmisi Mobil (Automatic vs Manual)

Kasus ini menggunakan dataset mtcars, sebuah dataset klasik dalam analisis regresi yang berisi data teknis dan performa berbagai mobil.

Salah satu variabel penting di dataset ini adalah am, yang merupakan variabel kategori (faktor) yang menunjukkan jenis transmisi mobil dengan nilai:

• 0 = Automatic (transmisi otomatis)
• 1 = Manual (transmisi manual)

Tujuan utama analisis adalah untuk memodelkan dan memprediksi kemungkinan sebuah mobil menggunakan transmisi manual berdasarkan variabel prediktor teknis, khususnya dua variabel:

• mpg (miles per gallon): Mengindikasikan tingkat efisiensi bahan bakar mobil, semakin besar berarti mobil lebih irit.
• hp (horsepower): Menunjukkan daya atau tenaga mesin mobil.

Analisis menggunakan regresi logistik karena variabel target (am) bersifat kategorikal biner (manual atau otomatis).

Melalui model ini, kita ingin mengetahui seberapa besar pengaruh mpg dan hp terhadap kemungkinan sebuah mobil memakai transmisi manual, serta mengevaluasi seberapa akurat model tersebut memperkirakan jenis transmisi berdasarkan karakteristik mobil.

# Load data  
data(mtcars)  
df <- mtcars  
df$am <- factor(df$am, labels = c("Automatic", "Manual"))  

# Tampilkan ringkasan data  
summary(df[, c("am", "mpg", "hp")])  

##          am          mpg              hp       
##  Automatic:19   Min.   :10.40   Min.   : 52.0  
##  Manual   :13   1st Qu.:15.43   1st Qu.: 96.5  
##                 Median :19.20   Median :123.0  
##                 Mean   :20.09   Mean   :146.7  
##                 3rd Qu.:22.80   3rd Qu.:180.0  
##                 Max.   :33.90   Max.   :335.0

# Fit model regresi logistik  
model_real <- glm(am ~ mpg + hp, data = df, family = binomial)  
summary(model_real)  

## 
## Call:
## glm(formula = am ~ mpg + hp, family = binomial, data = df)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -33.60517   15.07672  -2.229   0.0258 *
## mpg           1.25961    0.56747   2.220   0.0264 *
## hp            0.05504    0.02692   2.045   0.0409 *
## ---
## 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: 19.233  on 29  degrees of freedom
## AIC: 25.233
## 
## Number of Fisher Scoring iterations: 7

# Odds ratio  
exp_coef <- exp(coef(model_real))  
print(exp_coef)  

##  (Intercept)          mpg           hp 
## 2.543664e-15 3.524063e+00 1.056588e+00

# Prediksi probabilitas dan kelas  
df$prob <- predict(model_real, type = "response")  
df$pred_class <- ifelse(df$prob > 0.5, "Manual", "Automatic")  

# Visualisasi prediksi dan probabilitas  
if (!requireNamespace("ggplot2", quietly = TRUE))  
  install.packages("ggplot2")  
library(ggplot2)  

mpg_seq <- seq(min(df$mpg), max(df$mpg), length.out = 100)  
hp_seq <- seq(min(df$hp), max(df$hp), length.out = 100)  
grid <- expand.grid(mpg = mpg_seq, hp = hp_seq)  

grid$prob <- predict(model_real, newdata = grid, type = "response")  

ggplot() +  
  geom_tile(data = grid, aes(x = mpg, y = hp, fill = prob), alpha = 0.3) +  
  geom_point(data = df, aes(x = mpg, y = hp, color = pred_class, shape = am), size = 3) +  
  scale_fill_gradient(low = "blue", high = "red", name = "Prob Manual") +  
  labs(title = "Prediksi dan Probabilitas Transmisi Manual pada Mobil",  
       x = "Miles per Galon (mpg)",  
       y = "Tenaga Kuda (hp)",  
       color = "Prediksi",  
       shape = "Fakta") +  
  theme_minimal() +  
  theme(legend.position = "right")  

# Evaluasi dengan ROC dan AUC  
if (!requireNamespace("pROC", quietly = TRUE))  
  install.packages("pROC")  
library(pROC)  

## Type 'citation("pROC")' for a citation.

## 
## Attaching package: 'pROC'

## The following object is masked from 'package:gmodels':
## 
##     ci

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

roc_obj <- roc(response = df$am, predictor = df$prob, levels = rev(levels(df$am)))  

## Setting direction: controls > cases

auc_value <- auc(roc_obj)  
cat("AUC:", round(auc_value, 3), "\n")  

## AUC: 0.931

plot(roc_obj, main = "ROC Curve untuk Model Prediksi Transmisi Mobil")  
abline(a = 0, b = 1, col = "red", lty = 2)  
legend("bottomright", legend = paste("AUC =", round(auc_value, 3)), col = "black", lwd = 2)  

# Confusion matrix  
print(table(Predicted = df$pred_class, Actual = df$am))

##            Actual
## Predicted   Automatic Manual
##   Automatic        16      3
##   Manual            3     10

---

Interpretasi Hasil Model Regresi Logistik

Model regresi logistik yang dipasang untuk memprediksi jenis transmisi mobil (manual atau automatic) berdasarkan variabel mpg (miles per gallon) dan hp (horsepower) memberikan hasil sebagai berikut.

Koefisien dan Signifikansi

Nilai estimasi koefisien pada model adalah:

• Intercept: -33.61 (p = 0.0258),
  
• mpg: 1.26 (p = 0.0264),
  
• hp: 0.055 (p = 0.0409).

Ketiga koefisien tersebut memiliki nilai p kurang dari 0.05, sehingga semua variabel signifikan secara statistik pada tingkat kepercayaan 95%. Artinya, terdapat bukti yang cukup untuk menyatakan bahwa mpg dan hp berpengaruh signifikan terhadap peluang mobil menggunakan transmisi manual.

Interpretasi Koefisien

• Koefisien mpg sebesar 1.26 berarti setiap kenaikan 1 mil per galon meningkatkan log-odds peluang menggunakan transmisi manual sebesar 1.26. Dalam istilah odds ratio, peluangnya meningkat sekitar 3.5 kali lipat. Hal ini menunjukkan mobil yang lebih irit bahan bakar lebih mungkin menggunakan transmisi manual.

• Koefisien hp sebesar 0.055 menunjukkan setiap kenaikan 1 tenaga kuda menaikkan log-odds transmisi manual sebesar 0.055. Dalam angka odds ratio, ini berarti peluang transmisi manual meningkat sekitar 5.6% untuk setiap tambahan 1 tenaga kuda.

Perlu diperhatikan bahwa secara intuitif hp yang lebih besar sering diasosiasikan dengan transmisi otomatis, namun hasil model ini menunjukkan pengaruh positif kecil. Hal ini bisa disebabkan oleh karakteristik data atau korelasi dengan variabel lain di dalam model.

Kinerja Model

• Nilai deviance residual turun dari 43.23 (model tanpa prediktor) menjadi 19.23 (dengan mpg dan hp), menandakan model lebih baik dibanding model tanpa prediktor.

• AIC (Akaike Information Criterion) sebesar 25.23 merupakan indikasi keseimbangan antara kompleksitas model dan goodness-of-fit.

• Model konvergen setelah 7 iterasi perhitungan.

Evaluasi Prediksi

Model menghasilkan AUC (Area Under Curve) sebesar 0.931. Nilai AUC ini di atas 0.9 menandakan kinerja model sangat baik dalam membedakan mobil transmisi manual dan otomatis, dengan tingkat sensitivitas dan spesifisitas yang tinggi.

Confusion Matrix

Prediksi model terhadap data asli menghasilkan matriks kebingungan sebagai berikut:

Prediksi  Aktual	Automatic	Manual
Automatic	16	3
Manual	3	10

• Dari 32 mobil, 26 terklasifikasi dengan benar, dan 6 salah klasifikasi.
  
• Akurasi keseluruhan sekitar 81.25%.
  
• Terdapat 3 false negative (manual dikira automatic), dan 3 false positive (automatic dikira manual).

Kesimpulan

• Model regresi logistik ini berhasil menangkap hubungan signifikan antara mpg dan hp dengan tipe transmisi mobil, dengan mpg memiliki pengaruh kuat.
  
• Kinerja model sangat baik dilihat dari AUC > 0.9 dan tingkat akurasi yang memadai.
  
• Interpretasi ini dapat menjadi dasar untuk pemahaman faktor-faktor yang mempengaruhi pilihan jenis transmisi dan juga untuk prediksi transmisi mobil berdasarkan variabel teknis tersebut.

---

8.6 Model Regresi Poisson

Regresi Poisson adalah metode statistik yang digunakan untuk memodelkan variabel respons yang berupa data diskrit berupa hitungan (count data), yaitu bilangan bulat non-negatif seperti jumlah kejadian, frekuensi, atau banyaknya suatu peristiwa dalam periode waktu atau area tertentu. Contohnya termasuk jumlah kecelakaan lalu lintas per bulan, banyaknya panggilan telepon yang diterima suatu call center per hari, atau jumlah pasien yang datang ke rumah sakit.

Model regresi ini merupakan bagian dari Generalized Linear Model (GLM), di mana data respons diasumsikan mengikuti distribusi Poisson. Berbeda dengan regresi linear biasa, model ini tidak memprediksi nilai terus menerus melainkan rata-rata jumlah kejadian yang bersifat diskrit, dengan menghubungkannya melalui fungsi link ke variabel prediktor.

Distribusi Poisson

Distribusi Poisson merupakan distribusi probabilitas yang paling umum digunakan untuk data hitungan, dengan fungsi probabilitas:

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

di mana:

• \(Y\) adalah variabel random jumlah kejadian,
• \(\lambda > 0\) adalah parameter rata-rata (mean) dan varians dari distribusi Poisson,
• \(y!\) adalah faktorial dari \(y\).

Distribusi ini memiliki ciri bahwa rata-rata dan variansnya sama, yaitu \(E(Y) = Var(Y) = \lambda\).

Regresi Poisson dalam Kerangka Exponential Family

Distribusi Poisson termasuk dalam keluarga exponential family, yang memungkinkan penerapan Generalized Linear Model (GLM). Bentuk umum fungsi kepadatan dari exponential family dapat ditulis sebagai:

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

Untuk distribusi Poisson:

• Parameter natural (canonical parameter) \(\theta = \log \lambda\),
• Fungsi skala \(a(\phi) = 1\) (parameter dispersi \(\phi = 1\)),
• Fungsi \(b(\theta) = e^\theta = \lambda\),
• Fungsi \(c(y, \phi) = - \log y!\).

Sehingga, distribusi Poisson dapat diekspresikan sebagai:

\[ f(y; \theta) = \exp\{ y \theta - e^\theta - \log y! \} \]

Fungsi Link dan Model Regresi

Pada regresi Poisson, hubungan antara variabel respons \(Y\) dan variabel prediktor \(X = (X_1, X_2, ..., X_p)\) dimodelkan melalui fungsi link kanonik, yaitu fungsi logaritma:

\[ g(\mu) = \log(\mu) \]

dengan

\[ \mu_i = E(Y_i) = \lambda_i, \quad i = 1, 2, ..., n \]

Sehingga modelnya menjadi:

\[ \log(\mu_i) = \eta_i = \mathbf{x}_i^\top \boldsymbol{\beta} \]

atau ekuivalen:

\[ \mu_i = \exp(\mathbf{x}_i^\top \boldsymbol{\beta}) \]

di mana:

• \(\mathbf{x}_i\) adalah vektor nilai prediktor pada observasi ke-\(i\),
• \(\boldsymbol{\beta} = (\beta_0, \beta_1, ..., \beta_p)^\top\) adalah vektor parameter regresi yang akan diestimasi,
• \(\eta_i\) adalah nilai prediktor linear (linear predictor),
• \(\mu_i\) adalah nilai ekspektasi variabel respons untuk observasi ke-\(i\).

Fungsi link log memastikan bahwa prediksi \(\mu_i\) selalu positif, sesuai kebutuhan data hitungan.

8.6.1 Estimasi Parameter

Parameter \(\boldsymbol{\beta}\) dalam model regresi Poisson diestimasi dengan metode Maximum Likelihood Estimation (MLE). Likelihood fungsi berdasarkan distribusi Poisson adalah:

\[ L(\boldsymbol{\beta}) = \prod_{i=1}^n \frac{e^{-\mu_i} \mu_i^{y_i}}{y_i!} \]

Log-likelihood fungsi (yang akan dimaksimalkan) adalah:

\[ l(\boldsymbol{\beta}) = \sum_{i=1}^n \big[y_i \log(\mu_i) - \mu_i - \log(y_i!) \big] \]

Dengan substitusi \(\mu_i = \exp(\mathbf{x}_i^\top \boldsymbol{\beta})\), diperoleh:

\[ l(\boldsymbol{\beta}) = \sum_{i=1}^n \big[ y_i (\mathbf{x}_i^\top \boldsymbol{\beta}) - \exp(\mathbf{x}_i^\top \boldsymbol{\beta}) - \log(y_i!) \big] \]

Karena log-likelihood tersebut tidak memiliki bentuk tertutup untuk solusi analitis, parameter \(\boldsymbol{\beta}\) diestimasi menggunakan metode numerik, seperti algoritma iteratif Newton-Raphson atau Fisher Scoring.

Asumsi Model Regresi Poisson

1. Variabel respons berupa bilangan bulat non-negatif (count data).
   
2. Kejadian bersifat independen, artinya setiap observasi adalah independen satu sama lain.
   
3. Rata-rata dan variansi data respons diasumsikan sama (ekspektasi dan varian sama \(\lambda\)) — ciri khas distribusi Poisson.
   
4. Fungsi link log menghubungkan rata-rata respons dan prediktor melalui hubungan linear pada skala logaritmik.

Implementasi dan Interpretasi

• Setelah estimasi parameter dilakukan, nilai \(\hat{\boldsymbol{\beta}}\) dapat diinterpretasikan sebagai efek logaritmik dari variabel prediktor terhadap nilai ekspektasi respons.
  
• Misalnya, koefisien regresi \(\beta_j\) dapat ditafsirkan bahwa kenaikan satu unit pada \(X_j\) akan mengalikan ekspektasi respons dengan faktor \(e^{\beta_j}\).
  
• Model ini banyak digunakan untuk analisis data frekuensi dalam bidang medis, ekologi, pemasaran, dan lain-lain.

Model regresi Poisson adalah alat penting untuk menganalisis data cacah. Ia memberikan hubungan log linear antara prediktor dan rata-rata kejadian. Namun, perlu diperhatikan kemungkinan overdispersion dalam penerapannya

---

8.6.2 Contoh Kasus 1:

Prediksi Jumlah Kunjungan Toko Berdasarkan Anggaran Iklan

Misalkan kita ingin memodelkan jumlah kunjungan ke sebuah toko (variabel count) berdasarkan anggaran iklan yang dikeluarkan (dalam ribuan dolar).

set.seed(123)
n <- 150
ad_budget <- runif(n, min = 0, max = 10)       # Anggaran iklan (0-10 ribu dolar)
lambda <- exp(0.5 + 0.4 * ad_budget)           # Mean jumlah kunjungan sesuai anggaran
visits <- rpois(n, lambda)                      # Simulasi banyak kunjungan toko
data <- data.frame(visits, ad_budget)

## Estimasi Regresi Poisson

poisson_model <- glm(visits ~ ad_budget, data = data, family = poisson)
summary(poisson_model)

## 
## Call:
## glm(formula = visits ~ ad_budget, family = poisson, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 0.553174   0.064969   8.514   <2e-16 ***
## ad_budget   0.394689   0.008135  48.518   <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: 3360.95  on 149  degrees of freedom
## Residual deviance:  148.75  on 148  degrees of freedom
## AIC: 799.61
## 
## Number of Fisher Scoring iterations: 4

# Plot Hasil Prediksi

plot(data$ad_budget, data$visits, 
     pch = 16, col = "darkgreen", 
     main = "Data Kunjungan Toko dan Prediksi Model Poisson",
     xlab = "Anggaran Iklan (ribuan dolar)",
     ylab = "Jumlah Kunjungan")

newdata <- data.frame(ad_budget = seq(min(data$ad_budget), max(data$ad_budget), length.out = 100))
pred <- predict(poisson_model, newdata = newdata, type = "response")
lines(newdata$ad_budget, pred, col = "blue", lwd = 2)

#Diagnostik dan Overdispersion


dispersion <- sum(residuals(poisson_model, type = "pearson")^2) / poisson_model$df.residual
print(dispersion)

## [1] 0.9692365

Berikut interpretasi output model regresi Poisson untuk kasus jumlah kunjungan toko berdasarkan anggaran iklan:

Koefisien dan Signifikansi

Variabel	Estimate	Std. Error	z value	Pr(>
( Intercept)	0.553174	0.064969	8.514	< 2e-16 (sangat signifikan)
ad_budget	0.394689	0.008135	48.518	< 2e-16 (sangat signifikan)

• Kedua koefisien signifikan secara statistik (p-value sangat kecil).
• Intercept 0.553174 menunjukkan log rata-rata kunjungan saat anggaran iklan = 0.
• Ad_budget 0.394689 berarti kenaikan 1 ribu dolar anggaran iklan meningkatkan log rata-rata kunjungan sebanyak 0.3947.

Interpretasi Koefisien dalam Skala Asli

• Dalam skala logaritma, koefisien ad_budget \(\beta = 0.3947\).
• Dalam skala rata-rata kunjungan:

\[ e^{0.3947} \approx 1.484 \]

Artinya, setiap penambahan 1 ribu dolar anggaran iklan mengalikan ekspektasi jumlah kunjungan sebanyak sekitar 1.48 kali (48% peningkatan).

Goodness of Fit

• Null deviance: 3360.95 (model dengan intercept saja)

• Residual deviance: 148.75 (model dengan prediktor ada)

• Penurunan deviance sangat besar menunjukkan model dengan ad_budget jauh lebih baik.

• AIC (Akaike Information Criterion) sebesar 799.61 — ini ukuran kualitas model, semakin kecil semakin baik (untuk perbandingan antar model).

• Model konvergen setelah 4 iterasi Fisher scoring.

Diagnostik Overdispersion

• Nilai dispersion = 0.9692 (mendekati 1)
• Artinya, asumsi varians sama dengan mean (karakteristik Poisson) terpenuhi dengan baik. Tidak ada indikasi overdispersion pada data ini.

Jika dispersion > 1, berarti ada indikasi overdispersion (varians data lebih besar dari rata-rata). Jika overdispersion terjadi, alternatif model seperti Negative Binomial Regression perlu dipertimbangkan.

Kesimpulan

• Anggaran iklan berpengaruh positif dan signifikan terhadap jumlah kunjungan toko.
  
• Setiap kenaikan anggaran sebesar 1 ribu dolar meningkatkan rata-rata kunjungan sekitar 48%.
  
• Model regresi Poisson yang dibuat cukup tepat dan tidak menunjukkan gejala overdispersion.
  
• Dengan model ini, prediksi rata-rata kunjungan toko dapat dilakukan secara efektif berdasarkan nilai anggaran iklan.

8.6.3 Contoh Kasus 2

Dataset quine berasal dari paket MASS dan berisi data tentang jumlah hari siswa absen dari sekolah. Variabel utama adalah jumlah hari absen (Days) yang merupakan data count (bilangan bulat non-negatif). Variabel prediktornya meliputi:

• Eth: etnis siswa (A = Aboriginal, N = Non-Aboriginal)
• Sex: jenis kelamin (M/F)
• Age: kategori usia (0, 1, 2, 3)
• Lrn: status pembelajaran (AL = Average Learner, SL = Slow Learner, ML = More Slow Learner)

Tujuannya adalah memodelkan hubungan antara karakteristik siswa dan jumlah hari absen dengan menggunakan model regresi Poisson.

Model ini cocok karena Days adalah data hitungan. Fungsi link logaritma menghubungkan rata-rata hari absen dengan variabel prediktor melalui model linier pada skala log.

Hasil model memungkinkan untuk menginterpretasi seberapa faktor-faktor seperti etnis, umur, jenis kelamin, dan status pembelajaran memengaruhi jumlah hari absen siswa. Misalnya, koefisien log-odds untuk kelompok tertentu menunjukkan besarnya peningkatan atau penurunan ekspektasi jumlah absen relatif terhadap kelompok referensi.

Model ini juga dapat dipakai untuk prediksi rata-rata hari absen berdasarkan profil siswa dan membantu dalam pengambilan keputusan pendidikan.

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

data("quine")
# Pemodelan Regresi Poisson

model_quine <- glm(Days ~ Eth + Sex + Age + Lrn, family = poisson, data = quine)
summary(model_quine)

## 
## Call:
## glm(formula = Days ~ Eth + Sex + Age + Lrn, family = poisson, 
##     data = quine)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.71538    0.06468  41.980  < 2e-16 ***
## EthN        -0.53360    0.04188 -12.740  < 2e-16 ***
## SexM         0.16160    0.04253   3.799 0.000145 ***
## AgeF1       -0.33390    0.07009  -4.764 1.90e-06 ***
## AgeF2        0.25783    0.06242   4.131 3.62e-05 ***
## AgeF3        0.42769    0.06769   6.319 2.64e-10 ***
## LrnSL        0.34894    0.05204   6.705 2.02e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2073.5  on 145  degrees of freedom
## Residual deviance: 1696.7  on 139  degrees of freedom
## AIC: 2299.2
## 
## Number of Fisher Scoring iterations: 5

#Visualisasi Prediksi vs Data

library(ggplot2)
quine$predicted <- predict(model_quine, type = "response")

ggplot(quine, aes(x = Age, y = Days, color = Eth)) +
  geom_jitter(width = 0.2, alpha = 0.6) +
  geom_point(aes(y = predicted), shape = 18, size = 3) +
  facet_wrap(~ Lrn) +
  labs(title = "Jumlah Hari Absen berdasarkan Usia, Etnis dan Status Pembelajaran",
       x = "Usia (kategori)",
       y = "Hari Absen",
       color = "Etnis") +
  theme_minimal()

#Diagnostik Overdispersion

dispersion <- sum(residuals(model_quine, type = "pearson")^2) / model_quine$df.residual
dispersion

## [1] 13.16684

#Jika nilai dispersion >> 1, berarti mungkin ada overdispersion dan alternatif model seperti negative binomial perlu dipertimbangkan.

Berikut interpretasi output model regresi Poisson untuk kasus jumlah hari absen dari dataset quine:

Koefisien dan Signifikansi

Variabel	Estimate	Std. Error	z value	Pr(>	z
(Intercept)	2.71538	0.06468	41.980	< 2e-16	*** (sangat signifikan)
EthN	-0.53360	0.04188	-12.740	< 2e-16	***
SexM	0.16160	0.04253	3.799	0.000145	***
AgeF1	-0.33390	0.07009	-4.764	1.90e-06	***
AgeF2	0.25783	0.06242	4.131	3.62e-05	***
AgeF3	0.42769	0.06769	6.319	2.64e-10	***
LrnSL	0.34894	0.05204	6.705	2.02e-11	***

• Semua variabel signifikan secara statistik (p-value < 0.001).

---

Interpretasi Koefisien (dalam skala log dan eksponensial)

• Intercept (2.71538): log rata-rata hari absen untuk kelompok referensi (Etnis A, perempuan, kategori AgeF0, dan pembelajar AL).
  

• EthN (-0.53360): siswa etnis Non-Aboriginal memiliki rata-rata hari absen yang secara log berkurang 0.5336 dibanding etnis Aboriginal (referensi). Dalam skala eksponensial:

  \[ e^{-0.5336} \approx 0.586, \]

  artinya rata-rata hari absen siswa etnis Non-Aboriginal sekitar 58.6% dari etnis Aboriginal, atau menurun sekitar 41.4%.

• SexM (0.16160): siswa laki-laki memiliki rata-rata hari absen 17.5% lebih tinggi dibanding perempuan (karena \(e^{0.1616} \approx 1.175\)).
  

• AgeF1 (-0.33390), AgeF2 (0.25783), AgeF3 (0.42769): perubahan rata-rata hari absen relatif terhadap kategori usia referensi.

  • AgeF1: lebih rendah (sekitar 71.6% dari referensi),
    
  • AgeF2 dan AgeF3: lebih tinggi (sekitar 129.4% dan 153.3%).
    

• LrnSL (0.34894): siswa dengan status pembelajaran SL rata-rata absen meningkat sekitar 41.7% dibanding AL (\(e^{0.34894} \approx 1.417\)).

Goodness of Fit

• Null deviance: 2073.5 (model intercept saja).
  
• Residual deviance: 1696.7, penurunan deviance cukup signifikan menunjukkan model lebih baik dengan variabel prediktor.
• AIC: 2299.2 — kriteria untuk seleksi model (semakin kecil semakin baik).

Diagnostik Overdispersion

• Nilai dispersion: 13.17 (jauh lebih besar dari 1)
  
• Ini menandakan overdispersion, artinya varians jumlah absen jauh lebih besar dari rata-rata yang diasumsikan oleh model Poisson.
  
• Overdispersion ini menunjukkan model Poisson standar kurang tepat, dan model alternatif seperti Negative Binomial Regression perlu dipertimbangkan untuk menangani variansi berlebih.

Kesimpulan Singkat

• Variabel etnis, jenis kelamin, usia, dan status pembelajaran berpengaruh signifikan terhadap lama ketidakhadiran siswa.
  
• Terdapat perbedaan signifikan dalam rata-rata hari absen antar kategori.
  
• Namun, asumsi varians sama dengan mean pada model Poisson tidak terpenuhi (terjadi overdispersion), sehingga disarankan menggunakan model Negative Binomial atau metode lain untuk hasil yang lebih akurat.

Berikut materi lengkap dan terperinci mengenai Inferensi dalam Generalized Linear Model (GLM) yang telah dimodifikasi dan dikembangkan agar lebih komprehensif namun tetap jelas.

---

Kesimpulan Penting dari Inferensi GLM

• Ekspektasi estimasi digunakan untuk mengukur ketidaks biasan estimator, sehingga hasil yang diperoleh dapat dipercaya sebagai representasi dari parameter sebenarnya.
  
• Varians estimasi mengindikasikan presisi dan digunakan dalam menyusun uji hipotesis serta interval kepercayaan.
  
• Distribusi asimtotik normal dari estimator \(\hat{\beta}\) menjadi landasan metode inferensi seperti uji Wald dan perhitungan p-value.
  
• Dalam GLM, varians observasi sangat bergantung pada bentuk distribusi eksponensial yang digunakan, sehingga memerlukan fungsi varians khusus untuk setiap model.
  
• Parameter dispersi \(\phi\) memengaruhi varians dan kadang-kadang perlu diestimasi atau diperhitungkan jika tidak bernilai 1, terutama pada distribusi yang mengizinkan overdispersi.

Dengan pemahaman teori inferensi ini, penerapan GLM menjadi lebih valid dan handal dalam pengambilan keputusan berdasarkan data statistik.

Berikut materi lanjutan mengenai Ekspektasi, Varians, dan Metode Penaksiran Parameter dalam Generalized Linear Model (GLM), dengan penyesuaian penyajian agar tetap lengkap dan jelas tanpa menyederhanakan konsep pentingnya.

---

8.7 Mencari Ekspektasi dan Varians dalam GLM

8.7.1 Ekspektasi dalam GLM

Ekspektasi adalah nilai rata-rata dari variabel acak \(Y\), yang secara matematis dapat diekspresikan melalui fungsi momen probabilitas:

\[ \mathbb{E}[Y] = \int y \, f(y; \theta) \, dy = \mu \]

di mana \(f(y; \theta)\) adalah fungsi densitas (atau fungsi massa) probabilitas dengan parameter \(\theta\), dan \(\mu\) adalah nilai harapan atau rata-rata dari \(Y\).

Untuk kelas distribusi eksponensial (exponential family) yang menjadi dasar GLM, fungsi logaritma dari fungsi densitasnya memiliki bentuk umum:

\[ \log f(y; \theta) = a(y) + b(\theta) y + c(\theta) \]

atau bisa juga dirumuskan sebagai:

\[ \log f(y; \theta) = y \theta - b(\theta) + c(y) \]

Dimana fungsi \(b(\theta)\) dan \(c(\theta)\) merupakan fungsi tertentu yang bergantung pada parameter dan data.

8.7.2 Fungsi Skor dan Ekspektasi Turunan Pertama

Misalkan \(\ell(\theta) = \log f(y; \theta)\) adalah fungsi log-likelihood. Turunan log-likelihood terhadap \(\theta\) sering disebut fungsi skor, yang ditulis sebagai:

\[ U(\theta) = \frac{\partial \ell}{\partial \theta} = y - b'(\theta) \]

di mana \(b'(\theta)\) adalah turunan pertama dari \(b(\theta)\) terhadap \(\theta\).

Karena fungsi skor adalah deret statistik yang digunakan untuk pengujian dan estimasi, nilai harapan dari fungsi skor ini memberikan:

\[ \mathbb{E}[U(\theta)] = \mathbb{E}[y - b'(\theta)] = \mu - b'(\theta) = 0 \]

sehingga:
  
  \[
    \mu = b'(\theta)

]

Artinya, nilai harapan \(Y\) berhubungan langsung dengan turunan pertama fungsi \(b(\theta)\).

8.7.3 Varians dalam GLM

Turunan kedua dari log-likelihood \(\ell(\theta)\) adalah:

\[ \frac{\partial^2 \ell}{\partial \theta^2} = -b''(\theta) \]

yang memberikan informasi tentang kelengkungan fungsi log-likelihood.

Varians \(Y\) kemudian dapat dituliskan sebagai:

\[ \text{Var}(Y) = b''(\theta) \cdot \phi \cdot V(\mu) \]

atau lebih sederhana:

\[ \text{Var}(Y) = \phi \cdot V(\mu) \]

di mana:

• \(b''(\theta)\) adalah turunan kedua dari \(b(\theta)\),
• \(\phi\) adalah parameter dispersi (dispersion parameter),
• \(V(\mu)\) adalah fungsi varians yang spesifik untuk distribusi yang dipakai,

Fungsi \(V(\mu)\) merepresentasikan bentuk varians sebagai fungsi dari nilai rata-rata \(\mu\), yang merupakan ciri khas dari distribusi eksponensial.

---

8.8 Metode Penaksiran Parameter dalam GLM

Penaksiran parameter \(\beta\) pada GLM umumnya dilakukan dengan metode Maximum Likelihood Estimation (MLE). Prinsip dasar MLE adalah mencari nilai parameter yang memaksimalkan fungsi likelihood — yaitu fungsi yang menunjukkan seberapa besar kemungkinan data yang diamati terjadi berdasarkan model dan parameter.

Langkah-langkah dasar dalam MLE:

• Turunan Pertama (Score):

  Hitung turunan fungsi log-likelihood \(\frac{\partial \ell}{\partial \beta} = 0\) untuk mencari titik maksimum.

• Turunan Kedua:

  Pastikan turunan kedua negatif (\(\frac{\partial^2 \ell}{\partial \beta^2} < 0\)) agar titik tersebut adalah maksimum lokal.

Namun, mudahnya, dalam kasus GLM, bentuk solusi analitik sering tidak tersedia, sehingga diperlukan metode numerik.

Beberapa metode optimasi yang populer digunakan untuk mencari nilai \(\hat{\beta}\) yang memaksimalkan likelihood:

Newton-Raphson

• Menggunakan score vector (gradien fungsi log-likelihood).
• Menggunakan Hessian matrix (matriks turunan kedua).
• Iterasi dilakukan dengan formula:

\[ \beta^{(t+1)} = \beta^{(t)} - H^{-1}(\beta^{(t)}) \cdot U(\beta^{(t)}) \]

dimana:

• \(U(\beta^{(t)}) = \frac{\partial \ell(\beta^{(t)})}{\partial \beta}\) adalah score vector,
• \(H(\beta^{(t)}) = \frac{\partial^2 \ell(\beta^{(t)})}{\partial \beta \partial \beta^\top}\) adalah Hessian matrix.

Fisher Scoring

• Modifikasi dari metode Newton-Raphson,
• Mengganti Hessian matrix dengan matriks informasi Fisher, yang merupakan nilai harapan Hessian.
• Formula iterasi menjadi:

\[ \beta^{(t+1)} = \beta^{(t)} + I^{-1}(\beta^{(t)}) \cdot U(\beta^{(t)}) \]

dengan \(I(\beta^{(t)}) = \mathbb{E}[ -H(\beta^{(t)})]\).

Iteratively Reweighted Least Squares (IRLS)

• Metode ini adalah bentuk modifikasi Fisher Scoring,
• Pendekatan ini menyamakan estimasi GLM dengan masalah weighted least squares yang diperbaharui iteratif,
• Setiap iterasi merevisi bobot berdasarkan informasi varians,
• IRLS menawarkan stabilitas numerik dan kecepatan konvergensi dalam implementasi praktis GLM.

---

8.8.1 Implementasi Newton-Raphson dalam GLM

Langkah-langkah dasar Newton-Raphson dalam konteks GLM adalah sebagai berikut:

• Statistik skor pada parameter ke-j:

\[ U_j(\beta) = \frac{\partial \log L(\beta)}{\partial \beta_j} \]

• Matriks Hessian (turunan kedua):

\[ H_{jk}(\beta) = \frac{\partial^2 \log L(\beta)}{\partial \beta_j \partial \beta_k} \]

• Gunakan perluasan Taylor sekitar \(\beta^*\):

\[ U(\beta^*) \approx U(\beta) + H(\beta)(\beta^* - \beta) \]

• Perkirakan nilai parameter dengan iterasi:

\[ \hat{\beta} \approx \beta^{(t)} - H^{-1}(\beta^{(t)}) U(\beta^{(t)}) \]

Di mana iterasi dilakukan hingga konvergensi (perubahan parameter antar iterasi cukup kecil).

---

8.9 Diagnostik Model dalam Generalized Linear Model (GLM)

Diagnostik model bertujuan untuk mengevaluasi apakah model yang dibangun sudah tepat dan sesuai dengan data. Evaluasi ini bisa dilakukan dengan beberapa pendekatan, antara lain:

• Uji Formal: Pengujian statistik tertentu untuk menilai kesesuaian model.
• Visualisasi: Grafik perbandingan nilai prediksi model dengan nilai aktual dari data.

8.9.1 Statistik Devians

Devians adalah ukuran kesesuaian model yang digunakan untuk membandingkan model yang sedang diuji dengan model yang paling kompleks (saturated model) yang memiliki jumlah parameter sama dengan jumlah observasi.

• Devians mengukur perbedaan antara model yang diuji dengan saturated model.
• Jika nilai deviannya besar, ini menunjukkan model kurang cocok.
• Rumus deviannya adalah:

\[ D = 2 \sum_i \left[ y_i \log\left(\frac{y_i}{\hat{\mu}_i}\right) - (y_i - \hat{\mu}_i) \right] \]

di mana \(y_i\) adalah nilai observasi dan \(\hat{\mu}_i\) nilai prediksi model.

• Devians yang kecil mengindikasikan model yang lebih baik sesuai data.

8.9.2 Statistik Chi-Kuadrat Pearson

Statistik ini digunakan untuk menguji apakah model lebih baik dibandingkan dengan tidak menggunakan model (null model).

• Rumus statistiknya:

\[ X^2 = \sum_i \frac{(y_i - \hat{\mu}_i)^2}{\hat{\mu}_i} \]

• Jika nilai statistik signifikan (berdasarkan distribusi chi-kuadrat), maka model yang digunakan lebih baik dibandingkan tanpa model sama sekali.

Catatan Penting

• Untuk data yang terkelompok (clustered data), nilai devians dan chi-kuadrat Pearson mengikuti distribusi chi-kuadrat dengan derajat kebebasan tertentu.
• Untuk data tidak terkelompok, distribusi tersebut mungkin tidak berlaku secara ketat.
• Devians diminimalkan oleh metoda MLE, sehingga cocok digunakan untuk evaluasi kualitas model GLM.

8.9.3 Analisis Residual

• Residual adalah selisih antara nilai observasi dengan nilai prediksi.
• Analisis residual dapat mendeteksi penyimpangan sistematis dari model.
• Plot residual membantu menilai asumsi model seperti linearitas, homoskedastisitas, dan distribusi error.

Berikut lanjutan materi tentang Diagnostik Model GLM dan Detail Metode Estimasi serta Inferensi Regresi Logistik dengan penyesuaian agar lebih jelas dan rinci tanpa disingkat:

8.10 Detail Metode Estimasi dan Inferensi pada Regresi Logistik

Regresi logistik dipakai untuk memodelkan probabilitas kejadian dari variabel respons biner (0 atau 1) berdasarkan satu atau lebih variabel prediktor.

Estimasi parameter dilakukan dengan Maximum Likelihood Estimation (MLE) karena model regresi logistik tidak linear pada parameternya.

8.10.1 Fungsi Model Logistik

Fungsi probabilitas dari keluaran model logistik adalah:

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

yang menggambarkan probabilitas respons \(Y=1\) berdasarkan variabel prediktor \(x\).

8.10.2 Fungsi Log-Likelihood

Untuk \(n\) observasi, fungsi log-likelihood adalah:

\[ \ell(\beta) = \sum_{i=1}^n \left[ y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \right] \]

dengan:

\[ \pi_i = P(Y_i=1 | X_i) = \frac{1}{1 + \exp(-X_i^\top \beta)} \]

8.10.3 Estimasi Parameter dengan Metode Newton-Raphson

Metode Newton-Raphson digunakan untuk mencari parameter \(\beta\) yang memaksimalkan fungsi log-likelihood tersebut dengan langkah iteratif sebagai berikut:

1. Turunan Pertama (Score Function):

\[ U(\beta) = \frac{\partial \ell(\beta)}{\partial \beta} = X^\top (y - \pi) \]

di mana \(y\) adalah vektor observasi dan \(\pi\) adalah vektor probabilitas model.

2. Turunan Kedua (Hessian Matrix):

\[ H(\beta) = \frac{\partial^2 \ell(\beta)}{\partial \beta \partial \beta^\top} = -X^\top W X \]

dengan \(W = \text{diag}(\pi_i (1 - \pi_i))\) adalah matriks bobot diagonal berdasarkan probabilitas prediksi.

3. Iterasi Newton-Raphson:

Nilai parameter diperbarui menggunakan:

\[ \beta^{(t+1)} = \beta^{(t)} + (X^\top W X)^{-1} X^\top (y - \pi^{(t)}) \]

dengan \(\pi^{(t)} = \frac{1}{1 + \exp(-X \beta^{(t)})}\).

Iterasi dilakukan sampai konvergen, yaitu nilai \(\beta\) tidak berubah signifikan antara iterasi.

Kesimpulan

• Newton-Raphson merupakan metode iteratif efisien untuk estimasi parameter dalam regresi logistik.
• Fungsi score dan Hessian memainkan peran utama dalam arah dan ukuran langkah estimasi.
• Diagnostik dan validasi model GLM penting untuk memastikan kualitas dan kesesuaian model.

---

Contoh Implementasi Newton-Raphson Regresi Logistik (R)

# Data simulasi
set.seed(2025)
n <- 100
x <- rnorm(n)
beta_true <- c(-0.5, 1.2)
X <- cbind(1, x)  # desain matriks dengan intercept

# Fungsi logit dan probabilitas
logit <- function(beta, X) {
  eta <- X %*% beta
  1 / (1 + exp(-eta))
}

# Simulasi response binomial
prob <- logit(beta_true, X)
y <- rbinom(n, size = 1, prob = prob)

# Newton-Raphson Manual
beta <- c(0,0)       # inisialisasi parameter
tol <- 1e-6          # toleransi konvergensi
max_iter <- 100
for (i in 1:max_iter) {
  pi <- logit(beta, X)
  W <- diag(as.vector(pi * (1 - pi)))
  score <- t(X) %*% (y - pi)
  Hessian <- -t(X) %*% W %*% X
  step <- solve(Hessian, score)
  
  beta_new <- beta - step
  
  if (max(abs(beta_new - beta)) < tol) {
    cat("Konvergen setelah", i, "iterasi\n")
    break
  }
  beta <- beta_new
}

## Konvergen setelah 5 iterasi

cat("Estimasi parameter:\n")

## Estimasi parameter:

print(beta)

##         [,1]
##   -0.3841078
## x  1.0170038

# Bandingkan dengan hasil glm standar
model_glm <- glm(y ~ x, family = binomial)
summary(model_glm)

## 
## Call:
## glm(formula = y ~ x, family = binomial)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3841     0.2261  -1.699   0.0893 .  
## x             1.0170     0.2605   3.904 9.48e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 136.06  on 99  degrees of freedom
## Residual deviance: 116.36  on 98  degrees of freedom
## AIC: 120.36
## 
## Number of Fisher Scoring iterations: 3

8.10.4 Analisis Residual dan Plot Diagnostik

Setelah model diestimasi, langkah selanjutnya untuk evaluasi:

# Prediksi probabilitas dan residual deviance
yhat <- predict(model_glm, type="response")
res_dev <- residuals(model_glm, type = "deviance")
res_pearson <- residuals(model_glm, type = "pearson")

# Plot residual deviance vs fitted values
plot(yhat, res_dev, main="Residual Deviance vs Predicted Probability",
     xlab="Predicted Probability", ylab="Deviance Residuals")
abline(h=0, col="red", lty=2)

# Plot residual Pearson vs fitted values
plot(yhat, res_pearson, main="Pearson Residual vs Predicted Probability",
     xlab="Predicted Probability", ylab="Pearson Residuals")
abline(h=0, col="blue", lty=2)

# Plot ROC curve (gunakan package pROC jika ingin lebih lanjut)
# library(pROC)
# roc_obj <- roc(y, yhat)
# plot(roc_obj, main="ROC Curve")

---

Interpretasi Berikut interpretasi hasil output regresi logistik tersebut secara lengkap:

1. Proses Estimasi dan Konvergensi

• Konvergen setelah 5 iterasi
  Metode Newton-Raphson berhasil mencapai solusi estimasi parameter dalam 5 iterasi, menandakan proses optimasi log-likelihood berjalan efisien dan stabil.

2. Estimasi Parameter

Parameter	Estimate	Std. Error	z value	Pr(>
Intercept	-0.3841	0.2261	-1.699	0.0893 (.)
x	1.0170	0.2605	3.904	9.48e-05 (***)

• Intercept (-0.3841)
  Nilai ini adalah log odds ketika \(x=0\). Artinya, peluang baseline \(P(Y=1|x=0) = \exp(-0.3841) / (1 + \exp(-0.3841)) \approx 0.405\).

• Koefisien \(x\) (1.0170)
  Setiap kenaikan satu satuan pada \(x\), log-odds keberhasilan (\(Y=1\)) meningkat sebesar 1.017. Dalam istilah odds ratio:

  \[ \text{Odds Ratio} = e^{1.017} \approx 2.77 \]

  Ini berarti peluang keberhasilan menjadi sekitar 2.77 kali lebih besar setiap kenaikan 1 unit pada \(x\).

3. Signifikansi Statistik

• Koefisien \(x\) sangat signifikan dengan nilai p-value sangat kecil (\(\approx 0.0001\)), menandakan variabel \(x\) memiliki pengaruh kuat terhadap probabilitas respons.
• Intercept menunjukkan tanda signifikansi marginal (p-value ~0.089), yang bisa dianggap borderline tergantung konteks penelitian.

4. Diagnostik Model

• Null Deviance: 136.06 (df=99)
  Devians model tanpa variabel prediktor (hanya intercept).

• Residual Deviance: 116.36 (df=98)
  Devians setelah memasukkan variabel \(x\). Penurunan devians ini mengindikasikan model dengan \(x\) lebih baik daripada tanpa \(x\).

• AIC (Akaike Information Criterion): 120.36
  Ukuran kualitas model dengan penalti kompleksitas model; semakin kecil AIC, model semakin baik.

• Fisher Scoring Iterations: 3
  Menunjukkan estimasi parameter pada fungsi glm menggunakan metode Fisher scoring (mirip Newton-Raphson) membutuhkan 3 iterasi untuk konvergen.

5. Intisari

• Variabel \(x\) adalah prediktor signifikan yang meningkatkan peluang terjadinya outcome positif.
• Model dengan \(x\) secara statistik lebih baik dibandingkan model hanya intercept.
• Proses estimasi stabil dan efisien sesuai jumlah iterasi konvergensi.

Penjelasan:

• Newton-Raphson iteratif menghitung skor dan Hessian, lalu memperbarui parameter sampai konvergen.
• Estimasi parameter manual dan glm seharusnya menghasilkan nilai yang sangat mirip.
• Residual deviance dan Pearson residual dipakai untuk mengidentifikasi pola penyimpangan dan mengecek asumsi model.
• Plot residual membantu mendeteksi heteroskedastisitas atau ketidakcocokan model.
• Plot ROC (Receiver Operating Characteristic) bisa digunakan untuk mengevaluasi kemampuan prediksi model.

Interpretasi Plot

1. Plot Residual Deviance vs Predicted Probability

• Residual deviance menunjukkan selisih antara nilai observasi dan model yang diprediksi dalam bentuk deviance residual.
• Dari plot terlihat residual tersebar di sekitar garis nol (garis merah putus-putus) tapi terdapat pola jelas: residual positif dan negatif terpisah pada asal probabilitas tertentu.
• Pola ini menunjukkan adanya penyimpangan sistematis, yang menandakan kemungkinan model regresi logistik saat ini belum sepenuhnya menangkap hubungan antara \(x\) dan respons.
• Namun, tidak ada residual ekstrem yang sangat besar, sehingga model secara umum masih cukup layak.

2. Plot Pearson Residual vs Predicted Probability

• Residual Pearson mengukur penyimpangan terstandarisasi antara observasi dan prediksi.
• Plot ini memperlihatkan pola mirip dengan plot deviance residual, yaitu residual tersebar secara simetris di atas dan bawah garis nol (garis biru putus-putus).
• Pola sistematik ini mengindikasikan ada kemungkinan aspek model (seperti hubungan non-linear, variabel yang hilang, atau kesalahan distribusi) masih perlu diperbaiki.
• Tidak adanya titik residual yang sangat besar juga menunjukkan model tidak terlalu terdistorsi oleh outlier.

Kesimpulan untuk keduanya:

• Model sudah cukup baik dalam memprediksi probabilitas, tapi ada indikasi penyimpangan sistematis yang bisa jadi akibat model yang kurang lengkap atau prediktor belum merepresentasikan pola data sepenuhnya.
• Untuk perbaikan, dapat dilakukan misalnya:
  • Menambahkan variabel prediktor atau transformasi variabel (misal polinomial, interaksi),
  • Memeriksa outlier atau leverage points lebih lanjut,
  • Mempertimbangkan model non-linear atau alternatif lain.

Berikut penjelasan lengkap, terstruktur, dan detail tentang Inferensi Parameter pada Regresi Logistik dengan pendekatan uji Wald, likelihood ratio test, serta evaluasi model menggunakan AIC dan BIC. Materi ini adalah kelanjutan yang komprehensif tanpa ringkasan berlebihan.

8.11 Inferensi Parameter pada Regresi Logistik

Inferensi parameter bertujuan untuk menguji keabsahan dan signifikansi koefisien regresi \(\beta_j\) dalam model regresi logistik. Hal ini penting untuk mengetahui apakah variabel prediktor yang dimasukkan berpengaruh signifikan terhadap probabilitas kejadian.

8.11.1 Uji Wald

Tujuan
Menguji hipotesis nol:
\[ H_0: \beta_j = 0 \quad \text{(parameter tidak berpengaruh)} \] melawan alternatif
\[ H_1: \beta_j \neq 0 \]

Teori Uji Wald

• Dari teori estimasi Maximum Likelihood Estimation (MLE), estimator \(\hat{\beta}_j\) mendekati distribusi normal:
  \[ \hat{\beta}_j \sim N(\beta_j, \operatorname{Var}(\hat{\beta}_j)) \]

• Jika hipotesis nol benar (\(\beta_j = 0\)), maka statistik berikut:
  \[ Z = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \sim \mathcal{N}(0,1) \]

• Statistik Wald yang digunakan adalah kuadrat dari \(Z\):
  \[ W = Z^2 = \left(\frac{\hat{\beta}_j}{SE(\hat{\beta}_j)}\right)^2 \sim \chi^2_1 \]

Ini mengikuti distribusi Chi-Square dengan 1 derajat kebebasan.

Contoh Simulasi dan Uji Wald (R)

set.seed(123)
n <- 100
x <- rnorm(n)
log_odds <- -0.5 + 1.2 * x
p <- 1 / (1 + exp(-log_odds))
y <- rbinom(n, 1, p)
data <- data.frame(x, y)

model <- glm(y ~ x, data = data, family = binomial)
summary(model)

## 
## Call:
## glm(formula = y ~ x, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3097     0.2296  -1.349    0.177    
## x             1.2663     0.3080   4.111 3.94e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 137.99  on 99  degrees of freedom
## Residual deviance: 114.76  on 98  degrees of freedom
## AIC: 118.76
## 
## Number of Fisher Scoring iterations: 4

# Ambil koefisien dan standar error
beta_hat <- coef(model)["x"]
se_beta <- summary(model)$coefficients["x", "Std. Error"]

# Hitung statistik Z
Z <- beta_hat / se_beta
print(Z)

##        x 
## 4.110965

# Hitung statistik Wald
Wald_stat <- Z^2
print(Wald_stat)

##        x 
## 16.90003

8.11.2 Uji Likelihood Ratio (Likelihood Ratio Test / LRT)

Tujuan Menggunakan perbandingan antara model penuh dengan model nol (tanpa variabel prediktor) untuk menguji pengaruh variabel prediktor secara keseluruhan.

• Model nol: hanya intercept.
• Model penuh: mengandung variabel prediktor \(x\).

Langkah Pengujian dengan fungsi anova di R

model_null <- glm(y ~ 1, data = data, family = binomial)
anova(model_null, model, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: y ~ 1
## Model 2: y ~ x
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        99     137.99                          
## 2        98     114.76  1   23.229 1.438e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretasi

• Penurunan deviance yang signifikan dari model nol ke model penuh menunjukkan model dengan variabel prediktor lebih baik.
• Nilai p-value kecil (< 0.05) berarti variabel prediktor memberikan kontribusi signifikan ke model.

8.11.3 Evaluasi Kebaikan Model: AIC dan BIC

Untuk menilai kualitas dan kompleksitas model, digunakan:

• Akaike Information Criterion (AIC)
  Semakin kecil nilai AIC, semakin baik model.

• Bayesian Information Criterion (BIC)
  Penalti kompleksitas model lebih berat dibanding AIC, untuk memilih model yang sederhana tapi efektif.

Contoh Kode Evaluasi

AIC_value <- AIC(model)
BIC_value <- BIC(model)

cat("AIC model:", AIC_value, "\n")

## AIC model: 118.7598

cat("BIC model:", BIC_value, "\n")

## BIC model: 123.9701

# Hitung p-value
p_value <- 1 - pchisq(Wald_stat, df = 1)
print(p_value)

##            x 
## 3.940095e-05

8.12 Detail Metode Estimasi dan Inferensi Regresi Poisson

Model regresi Poisson digunakan khusus untuk memodelkan data hitungan (count data), yaitu jumlah kejadian yang mengikuti distribusi Poisson. Model ini cocok ketika variabel respons adalah data diskrit yang berupa jumlah kejadian dalam interval waktu, ruang, atau grup tertentu.

Model Regresi Poisson

Distribusi probabilitas variabel respons \(Y_i\) yang mengikuti distribusi Poisson dinyatakan dengan:

\[ P(Y_i = y_i) = \frac{e^{-\lambda_i} \lambda_i^{y_i}}{y_i!} \]

di mana \(\lambda_i\) adalah parameter rata-rata (mean) dan varians dari distribusi Poisson, yang menggambarkan intensitas kejadian.

Model regresi Poisson memodelkan secara logaritmik hubungan antara rata-rata kejadian dengan variabel prediktor:

\[ \log(\lambda_i) = \mathbf{x}_i^\top \boldsymbol{\beta} \quad \Rightarrow \quad \lambda_i = \exp(\mathbf{x}_i^\top \boldsymbol{\beta}) \]

8.12.1 Estimasi Parameter

Maximum Likelihood Estimation (MLE)

Fungsi log-likelihood untuk \(n\) observasi adalah:

\[ \ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[y_i \log(\lambda_i) - \lambda_i - \log(y_i!)\right] \]

dengan \(\lambda_i = \exp(\mathbf{x}_i^\top \boldsymbol{\beta})\).

Estimasi parameter \(\boldsymbol{\beta}\) diperoleh dengan memaksimalkan fungsi log-likelihood tersebut.

Metode Iteratif: Iteratively Reweighted Least Squares (IRLS)

Karena log-likelihood tidak dapat diselesaikan langsung, MLE diperkirakan secara iteratif menggunakan metode IRLS (Iteratively Reweighted Least Squares). Metode ini sangat mirip dengan Newton-Raphson, namun dikhususkan untuk model GLM dengan fungsi link log dan distribusi Poisson.

Langkah-Langkah IRLS Manual

1. Definisikan Model
   \[ \eta_i = \log(\lambda_i) = \mathbf{x}_i^\top \boldsymbol{\beta} \] Sehingga \(\lambda_i = \exp(\eta_i)\).

2. Formulasi Fungsi Log-Likelihood
   Maksimalkan: \[ \ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[y_i \log(\lambda_i) - \lambda_i - \log(y_i!)\right] \]

3. Formulasi Iteratif untuk Update Parameter
   Pada iterasi ke-\(t+1\), parameter diperbarui dengan:

   \[ \boldsymbol{\beta}^{(t+1)} = \left(\mathbf{X}^\top \mathbf{W}^{(t)} \mathbf{X}\right)^{-1} \mathbf{X}^\top \mathbf{W}^{(t)} \mathbf{z}^{(t)} \]

   dengan:

   • \(\mathbf{W} = \operatorname{diag}(\lambda_i)\), matriks bobot diagonal,
     
   • \(\mathbf{z} = \boldsymbol{\eta} + \frac{\mathbf{y} - \boldsymbol{\lambda}}{\boldsymbol{\lambda}}\) merupakan variabel yang di-“working response”,
     
   • \(\boldsymbol{\eta} = \mathbf{X} \boldsymbol{\beta}\).

4. Ulangi hingga konvergen (perubahan parameter kecil antara iterasi).

Contoh Simulasi Data dan Estimasi IRLS Manual (R)

set.seed(123)
n <- 100
x <- rnorm(n)
X <- cbind(Intercept=1, x)  # Desain matriks dengan intercept
beta_true <- c(0.5, 0.8)

eta <- X %*% beta_true
lambda <- exp(eta)
y <- rpois(n, lambda)

# Inisialisasi variabel untuk iterasi IRLS
beta <- c(0, 0)
tol <- 1e-6
max_iter <- 100

for (i in 1:max_iter) {
  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)) < tol) {
    cat("Konvergen pada iterasi ke-", i, "\n")
    break
  }
  
  beta <- beta_new
}

## Konvergen pada iterasi ke- 8

cat("Estimasi parameter:\n")

## Estimasi parameter:

print(beta)

##                [,1]
## Intercept 0.4494951
## x         0.8600048

Perbandingan dengan Fungsi glm() di R

model_glm <- glm(y ~ x, family = poisson)
summary(model_glm)

## 
## 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

Output model glm() hendaknya mendekati hasil estimasi IRLS manual, memberikan validasi atas metode iteratif yang dilakukan secara eksplisit.

8.13 Inferensi Parameter dan Pengujian Hipotesis

Setelah parameter diestimasi, perlu dilakukan inferensi untuk menguji signifikan tidaknya koefisien.

8.13.1 Uji Wald

• Ambil koefisien dan standar error dari model.
  
• Hitung statistik Z dan Wald:
  \[ Z = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)}, \quad W = Z^2 \]
• Hitung p-value dengan distribusi \(\chi^2_1\):
  \[ p = 1 - F_{\chi^2_1}(W) \]
• Jika p-value < 0.05, variabel prediktor signifikan.

** Contoh Kode Uji Wald**

coef_val <- coef(model_glm)[2]
se_val <- summary(model_glm)$coefficients[2, 2]

wald_z <- coef_val / se_val
wald_chisq <- wald_z^2
p_value <- 1 - pchisq(wald_chisq, df = 1)

cat("Z:", wald_z, "\nChi-Square:", wald_chisq, "\np-value:", p_value, "\n")

## Z: 11.52289 
## Chi-Square: 132.777 
## p-value: 0

8.13.2 Uji Likelihood Ratio

Bandingkan model penuh dengan model kosong menggunakan uji chi-square:

model_null <- glm(y ~ 1, family = poisson)
anova(model_null, model_glm, 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

8.13.3 Evaluasi Model: AIC dan BIC

• AIC (Akaike Information Criterion):
  Mengukur keseimbangan antara kompleksitas dan fit model. Nilai lebih kecil lebih baik.

• BIC (Bayesian Information Criterion):
  Penalti yang lebih berat terhadap kompleksitas model, mengarah pada model yang lebih sederhana.

aic_val <- AIC(model_glm)
bic_val <- BIC(model_glm)

cat("AIC model:", aic_val, "\nBIC model:", bic_val, "\n")

## AIC model: 325.7582 
## BIC model: 330.9685

Kesimpulan

• Estimasi parameter regresi Poisson dilakukan dengan MLE menggunakan pendekatan IRLS yang iteratif dan efisien.
  
• Estimasi parameter secara manual dengan IRLS mendekati hasil fungsi glm() di R.
  
• Inferensi parameter menggunakan uji Wald dan Likelihood Ratio (uji chi-square) untuk menguji kontribusi variabel prediktor.
  
• Kriteria AIC dan BIC berguna untuk evaluasi dan pemilihan model terbaik.

9. Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio

Regresi logistik merupakan metode statistik yang umum digunakan untuk menganalisis hubungan antara satu variabel respons biner (dua kategori, seperti “ya” atau “tidak”, “sukses” atau “gagal”) dengan satu atau lebih variabel prediktor. Keunggulan utama dari metode ini adalah kemampuannya untuk menangani variabel bebas dengan berbagai skala pengukuran—baik kategorikal maupun numerik—dan mengestimasi peluang suatu kejadian berdasarkan kombinasi karakteristik yang dimiliki.

Dalam dunia nyata, variabel prediktor yang digunakan dalam analisis sering kali berasal dari skala pengukuran yang berbeda. Oleh karena itu, penting untuk memahami bagaimana masing-masing jenis variabel tersebut diperlakukan dalam model regresi logistik, khususnya:

1. Nominal: Merupakan variabel kategorikal tanpa urutan logis antar kategorinya. Contohnya adalah jenis kelamin (laki-laki, perempuan), atau status pernikahan (lajang, menikah, cerai). Dalam model regresi logistik, variabel nominal diubah ke dalam bentuk variabel dummy agar dapat dimasukkan ke dalam model sebagai prediktor numerik.

2. Ordinal: Variabel yang memiliki urutan logis antar kategori, namun jarak antar kategori belum tentu sama. Contohnya adalah tingkat kepuasan (tidak puas, cukup puas, puas, sangat puas) atau tingkat pendidikan (SMA < Sarjana < Master < Doktor). Dalam regresi logistik, variabel ordinal dapat diperlakukan sebagai:

Nominal, dengan membuat dummy untuk masing-masing kategori. Numerik, dengan memberikan nilai skala (misalnya: SMA = 1, Sarjana = 2, dst.) untuk merepresentasikan urutan. Pilihan perlakuan ini harus mempertimbangkan asumsi model dan konteks substantif dari data.

3. Rasio: Merupakan variabel numerik kontinu yang memiliki nol absolut dan perbandingan antar nilai yang bermakna. Contoh: pendapatan (dalam juta rupiah), usia (dalam tahun), atau waktu (dalam detik). Variabel ini dapat langsung dimasukkan ke dalam model regresi logistik tanpa transformasi khusus, selama asumsi dasar terpenuhi.

Memahami perbedaan perlakuan terhadap ketiga jenis prediktor ini menjadi kunci dalam membangun model regresi logistik yang tidak hanya akurat secara statistik, tetapi juga relevan secara substantif. Pada bab ini, akan dibahas cara memperlakukan prediktor dengan skala yang berbeda dalam konteks regresi logistik, dilengkapi dengan ilustrasi melalui simulasi dan eksplorasi data untuk memperdalam pemahaman.

9.1 Simulasi Data

Untuk memahami bagaimana regresi logistik menangani prediktor dengan skala nominal, ordinal, dan rasio, kita akan membuat dataset simulasi.

Dalam dunia e-commerce, keputusan seseorang untuk membeli suatu produk dapat dipengaruhi oleh berbagai faktor seperti jenis kelamin, tingkat kepercayaan terhadap toko online, dan pendapatan bulanan. Kita akan membuat dataset simulasi untuk memodelkan peluang seseorang membeli produk online (purchase: 0 = tidak beli, 1 = beli) berdasarkan:

gender (nominal): Male / Female,

trust_level (ordinal): Low < Medium < High < Very High,

income (rasio): pendapatan bulanan.

# Load library
library(tibble)

# Set seed
set.seed(123)

# Jumlah observasi
n <- 500

# Simulasi variabel prediktor
gender <- sample(c("Male", "Female"), n, replace = TRUE)
trust_level <- sample(c("Low", "Medium", "High", "Very High"),
                      n, replace = TRUE, prob = c(0.25, 0.35, 0.25, 0.15))
income <- rnorm(n, mean = 6, sd = 2)  # dalam juta rupiah per bulan

# Model logit untuk probabilitas pembelian
logit_p <- -1.5 +
  0.6 * (gender == "Female") +
  1.0 * as.numeric(factor(trust_level, 
                          levels = c("Low", "Medium", "High", "Very High"),
                          ordered = TRUE)) +
  0.2 * income

# Ubah ke probabilitas
p <- 1 / (1 + exp(-logit_p))
 
# Simulasi variabel biner: purchase (0/1)
purchase <- rbinom(n, 1, p)

# Gabungkan ke dalam dataset
sim_data <- tibble(purchase, gender, trust_level, income)

# Tampilkan 6 data pertama
head(sim_data)

## # A tibble: 6 × 4
##   purchase gender trust_level income
##      <int> <chr>  <chr>        <dbl>
## 1        1 Male   High          4.80
## 2        1 Male   High          4.01
## 3        1 Male   Medium        8.05
## 4        1 Female Medium        7.50
## 5        1 Male   High          2.98
## 6        1 Female Medium        5.81

Interpretasi: Dataset ini berisi 500 observasi dengan empat variabel:

gender: jenis kelamin (nominal),

trust_level: tingkat kepercayaan terhadap toko online (ordinal),

income: pendapatan bulanan dalam juta rupiah (rasio),

purchase: keputusan membeli produk online (biner).

Nilai variabel purchase disimulasikan menggunakan model logit yang mempertimbangkan pengaruh gender, tingkat kepercayaan, dan pendapatan.

9.2 Eksplorasi Data

# Eksplorasi: jumlah dan rata-rata income berdasarkan status purchase
sim_data %>%
  dplyr::group_by(purchase) %>%
  dplyr::summarise(
    Jumlah = dplyr::n(),
    Rata2_Income = mean(income)
  )

## # A tibble: 2 × 3
##   purchase Jumlah Rata2_Income
##      <int>  <int>        <dbl>
## 1        0     65         5.69
## 2        1    435         6.04

# Jumlah dan rata-rata income berdasarkan status purchase dan gender
sim_data %>%
  dplyr::group_by(gender, purchase) %>%
  dplyr::summarise(
    Jumlah = dplyr::n(),
    Rata2_Income = mean(income),
    .groups = 'drop'
  )

## # A tibble: 4 × 4
##   gender purchase Jumlah Rata2_Income
##   <chr>     <int>  <int>        <dbl>
## 1 Female        0     22         5.96
## 2 Female        1    217         5.91
## 3 Male          0     43         5.56
## 4 Male          1    218         6.17

# Jumlah dan rata-rata income berdasarkan status purchase dan trust_level
sim_data %>%
  dplyr::group_by(trust_level, purchase) %>%
  dplyr::summarise(
    Jumlah = dplyr::n(),
    Rata2_Income = mean(income),
    .groups = 'drop'
  )

## # A tibble: 8 × 4
##   trust_level purchase Jumlah Rata2_Income
##   <chr>          <int>  <int>        <dbl>
## 1 High               0      4         3.77
## 2 High               1    130         5.91
## 3 Low                0     37         5.79
## 4 Low                1     87         6.07
## 5 Medium             0     23         5.99
## 6 Medium             1    145         6.13
## 7 Very High          0      1         2.76
## 8 Very High          1     73         6.06

9.3 Perlakuan Variabel Ordinal

9.3.1 1. Treat sebagai Nominal (Dummy)

Dalam bagian ini, kita akan melakukan analisis regresi logistik dengan memperlakukan variabel ordinal trust_level (tingkat kepercayaan) sebagai variabel nominal. Artinya, kita akan mengabaikan urutan level dari trust_level (“low”, “medium”, “high”) dan memperlakukannya sebagai kategori lepas yang tidak berurutan.

Pendekatan ini biasa digunakan untuk melihat pengaruh masing-masing level terhadap variabel target (purchase) tanpa mengasumsikan bahwa hubungan antara level bersifat linear atau berurutan. Kita akan membuat variabel dummy dengan baseline (kategori referensi) adalah “low”, dan membandingkan kategori lainnya terhadap baseline tersebut.

# Perlakuan trust_level sebagai faktor nominal (dummy)
# Memanggil package dplyr
library(dplyr)
# Cek distribusi trust_level
table(sim_data$trust_level)

## 
##      High       Low    Medium Very High 
##       134       124       168        74

# 10.3.1 Perlakuan Ordinal Sebagai Nominal
sim_data_nominal <- sim_data %>%
  mutate(trust_level = factor(trust_level, levels = c("Low", "Medium", "High", "Very High")))

# Model regresi logistik dengan dummy variabel
model_nominal <- glm(purchase ~ gender + trust_level + income,
                     data = sim_data_nominal, family = binomial)

# Lihat ringkasan model
summary(model_nominal)

## 
## Call:
## glm(formula = purchase ~ gender + trust_level + income, family = binomial, 
##     data = sim_data_nominal)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           0.64815    0.48022   1.350 0.177120    
## genderMale           -0.91063    0.29790  -3.057 0.002237 ** 
## trust_levelMedium     1.01141    0.30454   3.321 0.000896 ***
## trust_levelHigh       2.72969    0.55007   4.962 6.96e-07 ***
## trust_levelVery High  3.63966    1.03024   3.533 0.000411 ***
## income                0.11781    0.07114   1.656 0.097724 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 386.39  on 499  degrees of freedom
## Residual deviance: 319.85  on 494  degrees of freedom
## AIC: 331.85
## 
## Number of Fisher Scoring iterations: 6

Intercept (+0.64815)
Ini adalah log-odds dasar untuk individu berjenis kelamin Female, memiliki tingkat kepercayaan Low (sebagai referensi), dan income = 0.
Nilai p = 0.177120 (tidak signifikan), artinya baseline ini tidak berbeda secara signifikan dari probabilitas acak (50%).
Odds Ratio = exp(0.64815) = 1.91 → Peluang pembelian pada kondisi baseline adalah 1.91 kali dibandingkan odds dasar (namun tidak signifikan).

genderMale (-0.91063)
Individu berjenis kelamin laki-laki memiliki log-odds pembelian yang lebih rendah sebesar 0.91 dibandingkan perempuan.
Nilai p = 0.002237 (signifikan pada taraf 1%), menunjukkan bahwa jenis kelamin berpengaruh signifikan terhadap keputusan pembelian.
Odds Ratio = exp(-0.91063) = 0.40 → Peluang laki-laki untuk melakukan pembelian adalah sekitar 40% dari peluang perempuan.

trust_levelMedium (+1.01141)
Individu dengan tingkat kepercayaan Medium memiliki log-odds pembelian yang lebih tinggi sebesar 1.01 dibandingkan dengan tingkat kepercayaan Low.
Nilai p = 0.000896 (signifikan pada taraf 1%), menunjukkan bahwa trust level Medium berpengaruh signifikan.
Odds Ratio = exp(1.01141) = 2.75 → Peluang pembelian oleh individu dengan trust level Medium adalah 2.75 kali lebih besar dibandingkan trust level Low.

trust_levelHigh (+2.72969)
Individu dengan trust level High memiliki log-odds pembelian lebih tinggi sebesar 2.73 dibandingkan trust level Low.
Nilai p = 6.96e-07 (sangat signifikan), menunjukkan bahwa trust level High merupakan prediktor kuat.
Odds Ratio = exp(2.72969) = 15.32 → Peluang pembelian oleh individu dengan trust level High adalah 15.32 kali lebih besar dibanding trust level Low.

trust_levelVery High (+3.63966)
Individu dengan trust level Very High memiliki log-odds pembelian lebih tinggi sebesar 3.64 dibandingkan trust level Low.
Nilai p = 0.000411 (signifikan pada taraf 1%), menunjukkan efek positif yang sangat kuat.
Odds Ratio = exp(3.63966) = 38.03 → Peluang pembelian oleh individu dengan trust level Very High adalah sekitar 38 kali lebih besar dibandingkan trust level Low.

income (+0.11781)
Setiap kenaikan 1 unit pendapatan (dalam log income), meningkatkan log-odds pembelian sebesar 0.118.
Nilai p = 0.097724 (marginal signifikan di taraf 10%), menunjukkan bahwa pengaruh income terhadap pembelian relatif lemah namun tetap relevan.
Odds Ratio = exp(0.11781) = 1.125 → Setiap kenaikan 1 unit pendapatan meningkatkan peluang pembelian sekitar 12.5%.

Goodness-of-Fit
- Null deviance: 386.39 → deviance dari model tanpa prediktor.
- Residual deviance: 319.85 → deviance dari model dengan prediktor.
Penurunan deviance menunjukkan bahwa model membawa informasi penting dalam menjelaskan variabel target.
- AIC: 331.85 → Digunakan untuk membandingkan model; semakin kecil nilai AIC, semakin baik keseimbangan antara akurasi dan kompleksitas model.

Kesimpulan Praktis
- Tingkat kepercayaan (trust_level) merupakan prediktor paling kuat dalam menjelaskan peluang pembelian.
- Gender juga berpengaruh, dengan laki-laki cenderung memiliki peluang lebih rendah untuk melakukan pembelian.
- Income menunjukkan pengaruh yang lebih kecil, namun tetap positif terhadap pembelian.
- Model secara umum cukup baik, namun mungkin dapat ditingkatkan lebih lanjut dengan menambahkan variabel lain seperti usia, pengalaman belanja, atau preferensi produk.

9.3.2 2. Treat sebagai Rasio (Numeric Rank)

Pada bagian ini, variabel trust_level yang bersifat ordinal akan diperlakukan sebagai variabel numerik bertingkat, dengan asumsi bahwa kenaikan satu tingkat kepercayaan memiliki efek linier terhadap peluang pembelian. Pendekatan ini digunakan untuk menangkap pola hubungan bertingkat antara kepercayaan dan keputusan pembelian secara lebih efisien dalam model regresi logistik.

# Perlakuan trust_level sebagai ordinal numeric
sim_data_numeric <- sim_data %>%
  mutate(trust_numeric = case_when(
    trust_level == "Low" ~ 1,
    trust_level == "Medium" ~ 2,
    trust_level == "High" ~ 3,
    trust_level == "Very High" ~ 4
  ))

# Model regresi logistik dengan trust_numeric
model_numeric <- glm(purchase ~ gender + trust_numeric + income,
                     data = sim_data_numeric, family = binomial)

# Lihat ringkasan model
summary(model_numeric)

## 
## Call:
## glm(formula = purchase ~ gender + trust_numeric + income, family = binomial, 
##     data = sim_data_numeric)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -0.62616    0.55865  -1.121  0.26236    
## genderMale    -0.91139    0.29880  -3.050  0.00229 ** 
## trust_numeric  1.23119    0.19177   6.420 1.36e-10 ***
## income         0.11515    0.07125   1.616  0.10606    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 386.39  on 499  degrees of freedom
## Residual deviance: 320.87  on 496  degrees of freedom
## AIC: 328.87
## 
## Number of Fisher Scoring iterations: 6

Interpretasi

Model ini mengasumsikan bahwa kenaikan tingkat kepercayaan (trust_level) memiliki efek linear terhadap peluang membeli (purchase = 1).

Koefisien genderMale = -0.91139, dengan p-value = 0.00229 (signifikan pada level 1%). Artinya, jenis kelamin laki-laki cenderung memiliki kemungkinan membeli lebih rendah dibanding perempuan, jika variabel lain dikontrol.

Koefisien trust_numeric = 1.23119, dengan p-value < 0.001 (sangat signifikan). Ini menunjukkan bahwa setiap kenaikan satu tingkat kepercayaan (misalnya dari Medium ke High) meningkatkan log-odds pembelian sebesar 1.231. Dalam bentuk odds ratio: exp(1.23119) ≈ 3.425, sehingga peluang membeli meningkat sekitar 3,4 kali lipat setiap naik satu level kepercayaan.

Koefisien income = 0.11515, namun dengan p-value = 0.10606 (tidak signifikan). Artinya, pendapatan belum terbukti secara statistik memengaruhi keputusan membeli dalam model ini.

AIC model = 328.87, dan residual deviance berkurang dibanding null deviance, yang menunjukkan bahwa model ini cukup baik dalam menjelaskan variabel purchase.

9.4 Perbandingan Model

list(
 AIC_Nominal = AIC(model_nominal),
 AIC_Numeric = AIC(model_numeric)
 )

## $AIC_Nominal
## [1] 331.8491
## 
## $AIC_Numeric
## [1] 328.8699

Perbandingan Model

Nilai AIC untuk model dengan perlakuan trust sebagai nominal adalah 331.85, sedangkan untuk model dengan perlakuan trust sebagai ordinal numeric adalah 328.87.

Karena AIC model ordinal numeric lebih kecil, maka model dengan perlakuan trust_level sebagai ordinal lebih disukai. Ini menunjukkan bahwa asumsi adanya urutan pada tingkat kepercayaan (Low → Medium → High → Very High) memberikan informasi yang lebih baik dalam memprediksi keputusan pembelian dibandingkan jika dianggap sebagai kategori tanpa urutan.

Interpretasi Goodness-of-Fit

• Null deviance (386.39): Deviance dari model tanpa prediktor.

• Residual deviance (320.87): Deviance dari model setelah memasukkan prediktor (gender, trust_numeric, dan income).

• Penurunan nilai deviance menunjukkan bahwa prediktor memberikan informasi yang berarti dalam menjelaskan variasi pada variabel respon.

• AIC (328.87): Digunakan untuk membandingkan model. Karena AIC ini lebih rendah dibanding model nominal (331.85), maka model ordinal lebih efisien dalam hal keseimbangan antara goodness-of-fit dan kompleksitas model.

# Model null (tanpa prediktor)
nullmod <- glm(purchase ~ 1, data = sim_data, family = binomial)

# Menghitung McFadden R2 untuk model nominal dan numeric
r2_nominal <- 1 - (logLik(model_nominal) / logLik(nullmod))
r2_numeric <- 1 - (logLik(model_numeric) / logLik(nullmod))

# Menampilkan hasil
list(
  McFadden_R2_Nominal = r2_nominal,
  McFadden_R2_Numeric = r2_numeric
)

## $McFadden_R2_Nominal
## 'log Lik.' 0.1722046 (df=6)
## 
## $McFadden_R2_Numeric
## 'log Lik.' 0.1695627 (df=4)

McFadden’s R² mengukur goodness-of-fit. Semakin besar nilainya, semakin baik model memprediksi dibandingkan model null.

Nilai McFadden R² untuk model dengan trust_level diperlakukan sebagai kategori nominal adalah sekitar 0.172, sedangkan untuk model dengan trust_level diperlakukan sebagai variabel ordinal numerik adalah sekitar 0.170.

Nilai McFadden R² ini mengindikasikan seberapa baik model menjelaskan variasi data dibandingkan model tanpa prediktor (null model). Kedua model menunjukkan kemampuan yang cukup mirip, dengan model nominal sedikit lebih baik dalam menjelaskan variabilitas data.

Secara praktis, perbedaan nilai R² ini sangat kecil sehingga baik perlakuan variabel trust_level sebagai nominal maupun ordinal numeric memberikan hasil model yang hampir setara.

9.5 Visualisasi

Berikut ini adalah visualisasi prediksi probabilitas pembelian berdasarkan model regresi logistik yang menggunakan variabel trust_level sebagai variabel kategori nominal dan sebagai variabel numeric ordinal. Plot ini memperlihatkan bagaimana probabilitas prediksi berubah seiring dengan peningkatan income pada setiap kategori tingkat kepercayaan (trust_level). Dengan membandingkan kedua plot, kita dapat melihat perbedaan pengaruh perlakuan variabel ordinal sebagai nominal versus numeric terhadap hasil prediksi model.

library(ggplot2)
library(dplyr)

# Tambahkan prediksi probabilitas ke data
sim_data_nominal <- sim_data_nominal %>% 
  mutate(predicted = predict(model_nominal, type = "response"))

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

# Plot untuk model nominal (trust_level sebagai kategori nominal)
sim_data_nominal %>%
  ggplot(aes(x = income, y = predicted, color = trust_level)) +
  geom_point(alpha = 0.6) +
  labs(
    title = "Prediksi Probabilitas Pembelian (trust_level sebagai Nominal)",
    x = "Income (juta rupiah)",
    y = "Probabilitas Prediksi"
  ) +
  theme_minimal()

# Plot untuk model numeric (trust_level sebagai numeric ordinal)
sim_data_numeric %>%
  ggplot(aes(x = income, y = predicted, color = as.factor(trust_numeric))) +
  geom_point(alpha = 0.6) +
  labs(
    title = "Prediksi Probabilitas Pembelian (trust_level sebagai Numeric)",
    x = "Income (juta rupiah)",
    y = "Probabilitas Prediksi",
    color = "trust_numeric"
  ) +
  theme_minimal()

Jika yang diinginkan adalah model yang lebih sederhana dan interpretasi linear terhadap tingkat kepercayaan, perlakuan ordinal sebagai numeric lebih cocok.Namun apabila terdepat kecurigaan efek tingkat kepercayaan tidak linier atau tidak berurutan, model nominal mungkin dapat menangkapnya dengan lebih baik.

library(broom)
library(knitr)
library(kableExtra)

# Ringkasan model nominal
tidy(model_nominal) %>%
  kable(format = "html", caption = "Ringkasan Koefisien Model dengan Ordinal sebagai Nominal") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Ringkasan Koefisien Model dengan Ordinal sebagai Nominal

term	estimate	std.error	statistic	p.value
(Intercept)	0.6481459	0.4802236	1.349675	0.1771202
genderMale	-0.9106310	0.2979001	-3.056833	0.0022369
trust_levelMedium	1.0114089	0.3045352	3.321156	0.0008965
trust_levelHigh	2.7296891	0.5500741	4.962402	0.0000007
trust_levelVery High	3.6396595	1.0302419	3.532820	0.0004112
income	0.1178091	0.0711413	1.655987	0.0977244

# Ringkasan model numeric
tidy(model_numeric) %>%
  kable(format = "html", caption = "Ringkasan Koefisien Model dengan Ordinal sebagai Numeric") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Ringkasan Koefisien Model dengan Ordinal sebagai Numeric

term	estimate	std.error	statistic	p.value
(Intercept)	-0.6261555	0.5586501	-1.120837	0.2623574
genderMale	-0.9113859	0.2988014	-3.050139	0.0022874
trust_numeric	1.2311948	0.1917672	6.420258	0.0000000
income	0.1151542	0.0712522	1.616148	0.1060622

Interpretasi Ringkasan Koefisien Model dengan Ordinal sebagai Nominal

• Intercept (0.6481, p = 0.1771)
  Log-odds dasar untuk individu perempuan dengan kategori trust level “Low” dan income = 0. Tidak signifikan secara statistik, sehingga baseline ini tidak berbeda signifikan dari peluang 50%.

• genderMale (-0.9106, p = 0.0022)
  Laki-laki memiliki log-odds probabilitas sukses lebih rendah dibanding perempuan. Secara praktis, peluang sukses laki-laki lebih rendah dan signifikan.

• trust_levelMedium (1.0114, p = 0.0009), trust_levelHigh (2.7297, p < 0.001), trust_levelVery High (3.6397, p < 0.001)
  Semakin tinggi kategori tingkat kepercayaan (trust level), semakin besar log-odds peluang sukses. Efek ini sangat signifikan.

• income (0.1178, p = 0.0977)
  Pengaruh positif income terhadap peluang sukses, tetapi tidak signifikan pada tingkat 5%, hanya cenderung (p ~ 0.1).

Interpretasi Ringkasan Koefisien Model dengan Ordinal sebagai Numeric

• Intercept (-0.6262, p = 0.2624)
  Log-odds dasar untuk perempuan dengan trust level terendah (diberi nilai 1) dan income = 0. Tidak signifikan.

• genderMale (-0.9114, p = 0.0023)
  Laki-laki memiliki peluang sukses lebih rendah secara signifikan dibanding perempuan.

• trust_numeric (1.2312, p < 0.001)
  Setiap kenaikan satu tingkat trust level secara numerik meningkatkan log-odds peluang sukses secara signifikan.

• income (0.1152, p = 0.1061)
  Pengaruh income positif tapi tidak signifikan pada level 5%.

Kesimpulan

Kedua model menunjukkan bahwa gender dan trust level berpengaruh signifikan terhadap peluang sukses. Model nominal memungkinkan pengaruh setiap kategori trust level yang berbeda, sedangkan model numeric menganggap trust level sebagai skala linier. Income memberikan efek positif namun kurang signifikan pada kedua model.

---

10 Pemilihan Model Regresi Logistik dan Evaluasi

10.1 Membangun Model Regresi Logistik: Pendekatan Confirmatory dan Exploratory

Dalam membangun model regresi logistik, terdapat dua pendekatan utama yang sering digunakan oleh peneliti, yaitu pendekatan confirmatory dan exploratory. Keduanya memiliki tujuan dan filosofi yang berbeda, namun dapat saling melengkapi dalam proses analisis data yang menyeluruh.

1. Pendekatan Confirmatory (Konfirmatori)

Pendekatan ini berangkat dari suatu teori atau hipotesis yang telah ada sebelumnya. Peneliti telah memiliki dugaan atau pengetahuan awal mengenai variabel-variabel mana yang penting dan bagaimana arah hubungan antar variabel. Oleh karena itu, variabel independen dimasukkan ke dalam model secara langsung tanpa proses seleksi yang terlalu eksploratif.

Pendekatan confirmatory biasanya digunakan ketika:

• Peneliti ingin menguji hipotesis tertentu yang spesifik.

• Studi didasarkan pada teori atau literatur sebelumnya.

• Tujuan penelitian adalah konfirmasi atau validasi.

Kelebihannya adalah pendekatan ini lebih terstruktur, memiliki dasar teoritis yang kuat, dan cocok digunakan dalam studi inferensial. Namun, kekurangannya adalah pendekatan ini dapat melewatkan variabel-variabel yang secara statistik signifikan namun tidak dipertimbangkan sejak awal.

2. Pendekatan Exploratory (Eksploratori)

Berbeda dengan confirmatory, pendekatan exploratory digunakan ketika peneliti belum memiliki dugaan awal yang kuat, atau ketika ingin memahami struktur data secara lebih bebas. Model dibangun dengan mencoba berbagai kombinasi variabel prediktor, menggunakan teknik seperti stepwise regression, analisis korelasi awal, atau information criteria (AIC/BIC) untuk memilih model terbaik.

Pendekatan exploratory biasanya digunakan ketika:

• Peneliti ingin mengeksplorasi hubungan antar variabel.

• Tidak ada teori yang mendasari dengan kuat.

• Data yang digunakan bersifat baru atau eksploratif.

Kelebihannya adalah pendekatan ini dapat menemukan pola atau hubungan yang tidak terduga. Namun, ada risiko overfitting dan kesimpulan yang tidak generalisable jika tidak dilakukan dengan hati-hati.

Tujuan:

Tujuan utama dalam membangun model regresi logistik adalah menemukan model yang parsimonious, yaitu model yang cukup sederhana namun mampu menjelaskan data dengan baik dan memberikan prediksi yang akurat. Dalam konteks statistik, model yang parsimonious adalah model yang menggunakan jumlah parameter seminimal mungkin tanpa mengorbankan kualitas prediksi secara signifikan. Model yang terlalu kompleks cenderung mengalami overfitting, yaitu performa baik pada data latih namun buruk pada data baru. Sebaliknya, model yang terlalu sederhana bisa underfit, yaitu gagal menangkap pola penting dalam data.

Pemilihan antara pendekatan Confirmatory dan Exploratory sangat tergantung pada tujuan penelitian:

• Jika tujuan utama adalah menguji hipotesis tertentu berdasarkan teori atau temuan dari studi sebelumnya, maka pendekatan Confirmatory lebih sesuai. Dalam pendekatan ini, peneliti sudah memiliki dugaan awal tentang hubungan antara variabel, dan ingin melihat apakah data mendukung hipotesis tersebut. Pendekatan ini sangat umum dalam penelitian sosial, psikologi, dan ilmu kesehatan yang berbasis teori.

• Sebaliknya, jika tujuan utamanya adalah menemukan model terbaik berdasarkan pola-pola dalam data, maka pendekatan Exploratory lebih cocok. Pendekatan ini menggunakan data untuk mengidentifikasi variabel mana yang paling berkontribusi terhadap prediksi, tanpa terlalu bergantung pada teori yang sudah ada. Teknik seperti stepwise selection, perbandingan AIC/BIC, atau cross-validation sering digunakan dalam pendekatan ini untuk memilih kombinasi variabel yang optimal.

Dalam praktiknya, kedua pendekatan ini tidak harus dipisahkan secara kaku. Banyak peneliti menggunakan kombinasi keduanya, di mana teori digunakan sebagai dasar awal pemilihan variabel, namun analisis eksploratori dilakukan untuk menyempurnakan model. Pendekatan ini disebut pendekatan komplementer, dan sangat bermanfaat dalam meningkatkan keandalan dan ketepatan model sambil tetap menghormati dasar teoritis dari penelitian.

Dengan demikian, model akhir yang dihasilkan tidak hanya memiliki performa statistik yang baik, tetapi juga memiliki makna substantif dan interpretasi yang relevan dalam konteks penelitian.

---

Simulasi Data untuk Pemodelan Regresi Logistik

Untuk mempelajari proses pemilihan model dalam regresi logistik, kita akan membuat simulasi data yang merepresentasikan kemungkinan seseorang melakukan pembelian produk berdasarkan tiga prediktor:

• x1: skor minat terhadap produk (bernilai kontinu)

• x2: status langganan (0 = bukan pelanggan, 1 = pelanggan)

• x3: skor kepercayaan terhadap brand (kontinu)

Model logistik akan dibentuk menggunakan kombinasi linier dari ketiga prediktor tersebut.

# Memuat pustaka yang dibutuhkan
library(knitr)
library(dplyr)
library(ggplot2)
library(MASS)
library(caret)

## Loading required package: lattice

## 
## Attaching package: 'lattice'

## The following object is masked from 'package:gnm':
## 
##     barley

## 
## Attaching package: 'caret'

## The following objects are masked from 'package:DescTools':
## 
##     MAE, RMSE

library(pROC)
library(DescTools)

# Set seed untuk replikasi
set.seed(2025)

# Simulasi data
n <- 250
x1 <- rnorm(n, mean = 2, sd = 1.5)        # minat terhadap produk
x2 <- rbinom(n, 1, 0.4)                  # status pelanggan
x3 <- rnorm(n, mean = 3, sd = 1.2)        # kepercayaan terhadap brand

# Fungsi linier untuk logit
lin_pred <- -1.2 + 0.4 * x1 - 0.6 * x2 + 0.9 * x3
p <- 1 / (1 + exp(-lin_pred))            # Probabilitas melalui logit

# Variabel respons biner
y <- rbinom(n, 1, p)

# Buat data frame
df <- data.frame(y = as.factor(y), x1, x2, x3)

# Tampilkan sebagian data
head(df)

##   y       x1 x2        x3
## 1 1 2.931135  1 2.6084312
## 2 1 2.053462  0 4.4470019
## 3 1 3.159732  1 3.6767665
## 4 1 3.908734  1 2.8496895
## 5 1 2.556463  1 2.6342956
## 6 1 1.755718  1 0.9410421

---

Pemilihan Model: Full Model

Setelah data simulasi terbentuk, langkah pertama dalam pemodelan regresi logistik adalah membangun model penuh (full model), yang mencakup semua variabel prediktor yang tersedia. Tujuan dari pembuatan model ini adalah untuk mengevaluasi kontribusi masing-masing variabel terhadap probabilitas hasil (dalam hal ini, apakah seseorang membeli produk atau tidak).

Model penuh sering kali digunakan sebagai dasar pembanding untuk pendekatan stepwise selection ataupun pendekatan confirmatory berdasarkan teori.

# Membuat model penuh dengan semua prediktor
model_full <- glm(y ~ x1 + x2 + x3, data = df, family = binomial)

# Menampilkan ringkasan hasil model
summary(model_full)

## 
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.2048     0.6005  -2.006  0.04485 *  
## x1            0.4157     0.1494   2.783  0.00538 ** 
## x2           -0.3685     0.4064  -0.907  0.36457    
## x3            0.9304     0.2015   4.617  3.9e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 202.48  on 249  degrees of freedom
## Residual deviance: 168.44  on 246  degrees of freedom
## AIC: 176.44
## 
## Number of Fisher Scoring iterations: 6

---

10.2 Metode Stepwise: Forward, Backward, dan Kedua Arah

Setelah model penuh dibangun, langkah selanjutnya adalah melakukan seleksi variabel untuk menemukan model yang lebih sederhana namun tetap memiliki performa yang baik. Salah satu pendekatan yang umum digunakan adalah metode stepwise regression, yang terdiri dari:

• Forward Selection: dimulai dari model kosong (hanya intercept), lalu menambahkan variabel satu per satu berdasarkan peningkatan performa (penurunan AIC).
• Backward Elimination: dimulai dari model penuh, lalu menghapus variabel yang tidak signifikan satu per satu.
• Stepwise (Both Directions): kombinasi dari forward dan backward, memungkinkan penambahan dan penghapusan variabel secara iteratif.

Metode ini sering digunakan dalam pendekatan eksploratori, terutama ketika kita belum yakin variabel mana yang paling relevan secara statistik.

# Model kosong (null model) dan model penuh
null_model <- glm(y ~ 1, data = df, family = binomial)
model_full <- glm(y ~ x1 + x2 + x3, data = df, family = binomial)

# Stepwise Selection
step_forward <- step(null_model, direction = "forward", scope = formula(model_full), trace = FALSE)
step_backward <- step(model_full, direction = "backward", trace = FALSE)
step_both <- step(null_model, direction = "both", scope = formula(model_full), trace = FALSE)

# Bandingkan AIC dari tiap model
AIC(model_full, step_forward, step_backward, step_both)

##               df      AIC
## model_full     4 176.4387
## step_forward   3 175.2579
## step_backward  3 175.2579
## step_both      3 175.2579

Hasil seleksi model menggunakan metode stepwise, baik dengan pendekatan forward, backward, maupun kedua arah, menunjukkan bahwa model yang dihasilkan menggunakan tiga variabel prediktor memiliki nilai AIC yang lebih rendah dibandingkan dengan model penuh yang menggunakan empat variabel. Nilai AIC yang lebih kecil ini mengindikasikan bahwa model dengan tiga variabel tersebut memberikan keseimbangan yang lebih baik antara kompleksitas model dan kemampuan fit data.

Dengan demikian, model yang lebih sederhana ini cukup efektif dalam menjelaskan variasi pada data tanpa mengorbankan performa prediksi secara signifikan. Oleh karena itu, model hasil seleksi stepwise dapat dianggap lebih parsimonious dan lebih direkomendasikan untuk digunakan dibandingkan model penuh.

---

10.3 Evaluasi Model dengan Kurva ROC dan AUC

Setelah memilih model terbaik menggunakan metode stepwise, langkah selanjutnya adalah mengevaluasi performa model tersebut dalam memprediksi kelas target. Salah satu cara yang umum digunakan adalah dengan menggunakan kurva ROC (Receiver Operating Characteristic) dan menghitung nilai AUC (Area Under the Curve).

Kurva ROC memvisualisasikan trade-off antara sensitivitas (true positive rate) dan 1 - spesifisitas (false positive rate) pada berbagai threshold klasifikasi. Semakin mendekati sudut kiri atas, semakin baik kemampuan model dalam membedakan kelas positif dan negatif. Nilai AUC mengukur luas di bawah kurva ROC; nilai AUC yang lebih dekat ke 1 menunjukkan model dengan performa prediksi yang sangat baik.

Berikut adalah sintaks untuk membuat kurva ROC dari model hasil seleksi stepwise dan menampilkan plotnya:

pred_prob <- predict(step_both, type = "response")
roc_obj <- roc(df$y, pred_prob)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

plot(roc_obj, main = "Kurva ROC", col = "blue")

auc(roc_obj)

## Area under the curve: 0.7855

Kurva ROC yang diplot menunjukkan kemampuan model dalam membedakan antara kelas positif dan negatif. Nilai AUC sebesar 0.7855 mengindikasikan model memiliki performa prediksi yang baik.

Kriteria nilai AUC secara umum adalah sebagai berikut:

• 0.5 - 0.6: Performa model buruk, hampir seperti tebakan acak

• 0.6 - 0.7: Performa model cukup

• 0.7 - 0.8: Performa model baik

• 0.8 - 0.9: Performa model sangat baik

• > 0.9: Performa model hampir sempurna

Karena nilai AUC model ini adalah 0.7855, maka model termasuk dalam kategori baik dan cukup dapat diandalkan untuk klasifikasi pada data ini.

---

10.4 Pseudo R-Squared

Pseudo R-squared adalah ukuran alternatif untuk menilai goodness-of-fit model regresi logistik. Karena model logistik tidak menggunakan varian residual seperti pada regresi linear biasa, pseudo R-squared memberikan gambaran seberapa baik model menjelaskan variasi data.

Beberapa jenis pseudo R-squared yang umum digunakan antara lain:

1. Cox-Snell: Mengukur proporsi variasi yang dijelaskan oleh model, namun nilainya tidak bisa mencapai 1.

2. Nagelkerke: Modifikasi dari Cox-Snell yang skala nilainya dari 0 hingga 1, sehingga lebih mudah diinterpretasikan.

3. McFadden: Berdasarkan perbandingan log-likelihood model dengan model nol, nilai yang mendekati 1 menunjukkan model yang sangat baik.

Berikut adalah cara menghitung pseudo R-squared untuk model hasil seleksi stepwise:

library(DescTools)
PseudoR2(step_both, which = c("CoxSnell", "Nagelkerke", "McFadden"))

##   CoxSnell Nagelkerke   McFadden 
##  0.1244434  0.2241781  0.1640833

Nilai pseudo R-squared memberikan gambaran seberapa baik model regresi logistik dalam menjelaskan variasi data:

• Cox-Snell (0.124): Nilai ini menunjukkan model menjelaskan sekitar 12,4% variasi dalam data. Namun, Cox-Snell tidak dapat mencapai nilai maksimal 1 sehingga interpretasinya agak terbatas.

• Nagelkerke (0.224): Ini adalah versi modifikasi dari Cox-Snell yang memiliki rentang dari 0 sampai 1. Nilai 0,224 menandakan model dapat menjelaskan 22,4% variasi. Nilai Nagelkerke yang lebih tinggi menunjukkan model yang lebih baik. Secara umum:

  • Nilai < 0.2: model lemah
  • 0.2 – 0.4: model sedang
  • 0.4: model kuat

• McFadden (0.164): Dihitung dari perbandingan log-likelihood model penuh dan model nol. Nilai antara 0.2–0.4 dianggap sebagai indikasi model yang baik. Dengan nilai 0,164, model ini menunjukkan kemampuan yang cukup, namun masih bisa ditingkatkan.

Kesimpulan:
Model yang dibangun memiliki kemampuan prediksi yang memadai dengan penjelasan variasi data yang masuk akal, terutama dilihat dari nilai Nagelkerke yang termasuk kategori sedang. Namun, tetap ada peluang untuk meningkatkan model agar performanya lebih optimal.

---

10.5 Tabel Klasifikasi dan Evaluasi

Setelah model regresi logistik dibangun dan probabilitas prediksi diperoleh, langkah selanjutnya adalah mengubah probabilitas tersebut menjadi kelas prediksi (1 atau 0) berdasarkan threshold, biasanya 0.5. Kemudian, kita dapat mengevaluasi performa model dengan menggunakan tabel klasifikasi (confusion matrix) yang memberikan informasi tentang prediksi benar dan salah.

# Mengubah probabilitas menjadi kelas prediksi dengan threshold 0.5
pred_class <- ifelse(pred_prob >= 0.5, 1, 0)

# Membuat confusion matrix untuk evaluasi
conf_matrix <- confusionMatrix(factor(pred_class), df$y, positive = "1")

# Menampilkan hasil confusion matrix dan statistik evaluasi
conf_matrix

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0   3   4
##          1  32 211
##                                           
##                Accuracy : 0.856           
##                  95% CI : (0.8063, 0.8971)
##     No Information Rate : 0.86            
##     P-Value [Acc > NIR] : 0.6154          
##                                           
##                   Kappa : 0.1009          
##                                           
##  Mcnemar's Test P-Value : 6.795e-06       
##                                           
##             Sensitivity : 0.98140         
##             Specificity : 0.08571         
##          Pos Pred Value : 0.86831         
##          Neg Pred Value : 0.42857         
##              Prevalence : 0.86000         
##          Detection Rate : 0.84400         
##    Detection Prevalence : 0.97200         
##       Balanced Accuracy : 0.53355         
##                                           
##        'Positive' Class : 1               
## 

Berdasarkan hasil confusion matrix, model menghasilkan akurasi sebesar 85.6%, yang berarti bahwa sekitar 86 dari 100 prediksi yang dilakukan oleh model adalah benar. Namun, akurasi saja tidak cukup untuk mengevaluasi model, terutama jika data tidak seimbang.

Sensitivity (Recall untuk kelas positif) sangat tinggi, yaitu 98.14%, menunjukkan bahwa model sangat baik dalam mendeteksi kasus positif (kelas 1). Ini berarti hampir semua individu yang termasuk dalam kelas 1 berhasil dikenali oleh model.

Sebaliknya, specificity sangat rendah, hanya sebesar 8.57%, yang menunjukkan bahwa model sangat lemah dalam mengenali kelas negatif (kelas 0). Artinya, sebagian besar individu yang sebenarnya kelas 0 justru diklasifikasikan sebagai kelas 1.

Positive Predictive Value (PPV) sebesar 86.83% menunjukkan bahwa sebagian besar prediksi kelas 1 memang benar. Sementara itu, Negative Predictive Value (NPV) hanya 42.86%, yang berarti bahwa prediksi untuk kelas 0 seringkali tidak tepat.

Balanced Accuracy sebesar 53.36% menunjukkan bahwa performa model terhadap kedua kelas masih tergolong rendah secara keseluruhan. Kappa sebesar 0.10 menunjukkan bahwa tingkat kesepakatan antara prediksi model dan data aktual hanya sedikit lebih baik dari tebakan acak.

Terakhir, hasil uji McNemar (p-value < 0.001) mengindikasikan bahwa terdapat perbedaan signifikan antara distribusi kesalahan tipe I dan tipe II, yang berarti model cenderung bias terhadap salah satu kelas.

Kriteria umum untuk menilai performa model klasifikasi:

• Akurasi > 80%: dianggap baik, tapi harus dilihat seimbangkah distribusi kelasnya.

• Sensitivity dan Specificity > 70%: menunjukkan performa deteksi yang cukup baik.

• Balanced Accuracy > 70%: performa model dikatakan seimbang antar kelas.

• Kappa > 0.6: tingkat kesepakatan substansial, sedangkan Kappa < 0.2 berarti lemah.

PPV dan NPV tinggi: menandakan keandalan prediksi model.

---

10.6 Metode Perbandingan Model dalam Regresi Logistik

Pada bagian ini akan disajikan cara membandingkan model regresi logistik menggunakan ukuran:

1. Deviance
   Mengukur ketidaksesuaian antara model dan data. Semakin kecil nilai deviance, semakin baik model dalam memodelkan data.

2. AIC (Akaike Information Criterion)
   Digunakan untuk memilih model terbaik dengan mempertimbangkan trade-off antara goodness-of-fit dan kompleksitas model. Nilai AIC yang lebih kecil menunjukkan model yang lebih baik.

3. Likelihood-Ratio Test
   Digunakan untuk membandingkan dua model yang saling bersarang. Jika nilai p kecil (misal < 0.05), maka model yang lebih kompleks secara signifikan lebih baik.

4. Prinsip Parsimony
   Dari dua model dengan performa serupa, model yang lebih sederhana lebih disukai karena lebih mudah diinterpretasi dan menghindari overfitting.

Contoh

# Memuat library
library(MASS)
library(broom)
library(DescTools)

# Simulasi data (dimodifikasi)
set.seed(1234)
n <- 250
x1 <- rnorm(n, mean = 5, sd = 2)
x2 <- rbinom(n, 1, 0.4)
x3 <- runif(n, min = 1, max = 10)
lin_pred <- -0.8 + 0.9 * x1 - 1.1 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)

# Pembuatan Model
model1 <- glm(y ~ x1, data = data, family = binomial)
model2 <- glm(y ~ x1 + x2, data = data, family = binomial)
model3 <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)

# Perbandingan AIC dan 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))
)

model_comp

##     Model      AIC Deviance
## 1 Model 1 23.92764 19.92764
## 2 Model 2 24.02842 18.02842
## 3 Model 3 20.47167 12.47167

Interpretasi Hasil Perbandingan Model

Berdasarkan hasil perbandingan tiga model regresi logistik yang ditunjukkan oleh nilai AIC dan Deviance, berikut penjelasannya:

Model 1 hanya memasukkan satu prediktor dan menghasilkan nilai AIC sebesar 23.93 dan deviance sebesar 19.93.
Model 2 menambahkan satu prediktor tambahan, tetapi AIC-nya justru sedikit meningkat menjadi 24.03 meskipun deviance-nya menurun menjadi 18.03.
Model 3, yang merupakan model paling lengkap dengan tiga prediktor, menghasilkan nilai AIC paling rendah yaitu 20.47 dan deviance paling rendah yaitu 12.47.

Kesimpulan:
Model 3 adalah model terbaik karena memiliki nilai AIC dan deviance yang paling rendah dibandingkan dua model lainnya. Ini menunjukkan bahwa Model 3 memberikan keseimbangan terbaik antara kompleksitas model dan kemampuan memodelkan data.

Namun, perlu diperhatikan bahwa perbedaan nilai AIC antar model tidak terlalu besar. Oleh karena itu, prinsip parsimony tetap penting untuk dipertimbangkan. Jika peningkatan performa tidak signifikan, model yang lebih sederhana (misalnya Model 1) mungkin tetap dipilih karena lebih mudah diinterpretasi dan lebih hemat parameter.

---

10.7 Likelihood-Ratio Test

Likelihood-Ratio Test adalah metode statistik untuk membandingkan dua model yang saling bersarang (nested models) dalam regresi logistik. Model dikatakan bersarang apabila model yang lebih sederhana merupakan bagian dari model yang lebih kompleks, yakni model kompleks memiliki semua variabel dari model sederhana ditambah satu atau lebih variabel tambahan.

Tujuan: Untuk mengetahui apakah penambahan satu atau lebih variabel prediktor secara signifikan meningkatkan kemampuan model dalam menjelaskan data.

Cara kerja:

Uji ini membandingkan nilai deviance dari dua model. Nilai p-value menunjukkan apakah perbedaan deviance signifikan.

Hipotesis:

- H₀: Model sederhana cukup baik (variabel tambahan tidak signifikan).

- H₁: Model kompleks lebih baik (variabel tambahan signifikan).

# Uji LRT antara Model 1 dan Model 2
anova(model1, model2, test = "LRT")

## Analysis of Deviance Table
## 
## Model 1: y ~ x1
## Model 2: y ~ x1 + x2
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       248     19.928                     
## 2       247     18.028  1   1.8992   0.1682

# Uji LRT antara Model 2 dan Model 3
anova(model2, model3, test = "LRT")

## Analysis of Deviance Table
## 
## Model 1: y ~ x1 + x2
## Model 2: y ~ x1 + x2 + x3
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       247     18.028                       
## 2       246     12.472  1   5.5568  0.01841 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretasi Uji Likelihood-Ratio (LRT)

Uji Likelihood-Ratio (LRT) digunakan untuk membandingkan dua model regresi logistik yang saling bertingkat (nested models) yakni, satu model merupakan bentuk penyederhanaan dari model lainnya.

• Perbandingan Model 1 dan Model 2
  
  • Model 1: y ~ x1
    
  • Model 2: y ~ x1 + x2
    
  • Deviance difference = 1.8992
    
  • Df = 1
    
  • p-value = 0.1682

Karena p-value > 0.05, maka tidak ada perbedaan signifikan antara Model 1 dan Model 2. Ini berarti penambahan variabel x2 tidak memberikan peningkatan signifikan terhadap kualitas model. Model 1 tetap dapat dipertahankan berdasarkan prinsip parsimony.

• Perbandingan Model 2 dan Model 3
  
  • Model 2: y ~ x1 + x2
    
  • Model 3: y ~ x1 + x2 + x3
    
  • Deviance difference = 5.5568
    
  • Df = 1
    
  • p-value = 0.01841

Karena p-value < 0.05, maka terdapat perbedaan signifikan antara Model 2 dan Model 3. Artinya, penambahan variabel x3 memberikan kontribusi signifikan terhadap peningkatan kualitas model.

Kesimpulan:
Model 3 adalah model terbaik karena memberikan peningkatan signifikan dan tetap mempertahankan prinsip parsimony (sederhana namun cukup menjelaskan data).

Karena x2 tidak signifikan saat ditambahkan ke x1, sementara x3 terbukti signifikan, disarankan untuk mengevaluasi Model 4 dengan formula y ~ x1 + x3. Tujuannya adalah untuk melihat apakah model yang lebih sederhana (tanpa x2) tetap mampu mempertahankan performa yang baik.
Jika Model 4 memiliki nilai AIC dan deviance yang sebanding atau hanya sedikit lebih tinggi dibanding Model 3, maka Model 4 akan lebih disukai berdasarkan prinsip parsimony (model sesederhana mungkin, tanpa kehilangan kualitas prediksi).

Langkah selanjutnya adalah membangun Model 4 dan membandingkannya dengan Model 3 melalui AIC, deviance, dan LRT.

# Membuat Model 4 dengan variabel x1 dan x3
model4 <- glm(y ~ x1 + x3, data = data, family = binomial)

# Uji Likelihood-Ratio Test antara Model 3 dan Model 4
anova(model4, model3, test = "LRT")

## Analysis of Deviance Table
## 
## Model 1: y ~ x1 + x3
## Model 2: y ~ x1 + x2 + x3
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       247     14.695                     
## 2       246     12.472  1   2.2231    0.136

Kesimpulan

Meskipun Model 3 memiliki nilai deviance dan AIC yang lebih rendah dibandingkan Model 4, namun berdasarkan hasil uji LRT, tidak terdapat perbedaan signifikan antara keduanya. Oleh karena itu, Model 4 dapat dipertimbangkan sebagai model alternatif yang lebih sederhana, sesuai dengan prinsip parsimony, yaitu memilih model yang paling sederhana yang mampu menjelaskan data dengan baik.

---

10.8 Detail ROC (Receiver Operating Characteristic)

Dalam analisis regresi logistik atau model klasifikasi biner lainnya, setelah membangun model dan mengevaluasi akurasinya melalui confusion matrix, kita juga membutuhkan alat evaluasi lain yang memberikan gambaran performa model secara keseluruhan pada berbagai nilai ambang (threshold). Salah satu alat evaluasi paling populer adalah Kurva ROC (Receiver Operating Characteristic) dan AUC (Area Under the Curve).

1. Apa itu Kurva ROC?
Kurva ROC adalah grafik yang menggambarkan kemampuan model klasifikasi dalam membedakan antara dua kelas (positif dan negatif) pada berbagai nilai cut-off.

ROC mengilustrasikan trade-off antara dua metrik penting:

• True Positive Rate (TPR) atau Sensitivitas

• False Positive Rate (FPR) atau 1 - Spesifisitas

ROC sangat berguna ketika:

• Dataset tidak seimbang (jumlah kelas 0 dan 1 berbeda jauh),

• Ingin membandingkan beberapa model klasifikasi sekaligus,

• Ingin menentukan ambang prediksi (threshold) yang optimal sesuai konteks aplikasi.

2. Sumbu pada Kurva ROC
Kurva ROC diplot dengan:

Sumbu X (horizontal): False Positive Rate (FPR)
\[ FPR = \frac{FP}{FP + TN} \] Sumbu Y (vertikal): True Positive Rate (TPR) atau Sensitivitas
\[ TPR = \frac{TP}{TP + FN} \]

ROC dibentuk dengan mengubah-ubah nilai threshold dari 0 ke 1 dan menghitung nilai TPR dan FPR untuk setiap nilai ambang tersebut.

3. Makna Bentuk Kurva ROC

1. Garis Diagonal (45°): Ini adalah baseline atau prediksi acak (random guess). Model yang prediksinya tidak lebih baik dari tebak-tebakan akan menghasilkan ROC mendekati garis ini.

2. Kurva Dekat Sudut Kiri Atas: Semakin tinggi kurva dan semakin dekat ke pojok kiri atas, semakin baik performa model dalam memisahkan kelas positif dan negatif.

3. Kurva ROC Ideal:

   -Melonjak vertikal ke atas (TPR naik cepat) sambil tetap mempertahankan FPR rendah.

   -Kemudian mendatar sejajar sumbu X hingga mencapai titik (1,1).

   -Bentuk seperti ini mencerminkan model yang memiliki sensitivitas tinggi dan kesalahan klasifikasi rendah.

4. Perubahan Cut-off dan Dampaknya
Setiap model klasifikasi memberikan skor probabilitas (antara 0 dan 1). ROC dibentuk dengan mengubah threshold dari 0 sampai 1 dan melihat perubahan klasifikasi yang terjadi.

• Jika threshold (cut-off) diturunkan:
  • Lebih banyak kasus diprediksi sebagai positif.
  • Sensitivitas (TPR) meningkat, karena makin banyak positif yang terdeteksi.
  • Tapi Spesifisitas menurun (lebih banyak negatif diprediksi positif → FP naik).
• Jika threshold dinaikkan:
  • Model menjadi lebih ketat dalam memprediksi positif.
  • Spesifisitas naik (negatif terdeteksi dengan baik), tetapi
  • Sensitivitas turun (positif yang terlewat makin banyak → FN naik).

---

5. Apa itu AUC (Area Under the Curve)?
AUC (Area Under the Curve) adalah nilai numerik yang mengukur luas di bawah kurva ROC. Nilai AUC berkisar antara 0 sampai 1.

Interpretasi:

• AUC = 1.0: Model sempurna, mampu memisahkan kelas 0 dan 1 tanpa kesalahan.

• AUC = 0.5: Performa model setara dengan tebak acak (random guess).

• AUC > 0.7: Model dianggap cukup baik. - AUC > 0.9: Model sangat baik.

AUC juga bisa diartikan sebagai probabilitas bahwa model akan memberikan skor probabilitas lebih tinggi pada kasus positif dibandingkan kasus negatif secara acak. Oleh karena itu, AUC dikenal juga sebagai Concordance Index.

---

6. Kegunaan Kurva ROC dan AUC
ROC dan AUC digunakan untuk:

• Mengevaluasi performa model secara menyeluruh, tanpa bergantung pada satu nilai threshold saja.

• Membandingkan beberapa model klasifikasi , model dengan AUC lebih tinggi dianggap lebih baik.

• Menentukan cut-off optimal berdasarkan kebutuhan bisnis atau klinis:

• Jika False Negative berisiko tinggi (misalnya deteksi kanker), maka threshold harus lebih rendah untuk meningkatkan sensitivitas.

• Jika False Positive berisiko tinggi (misalnya tuduhan kriminal palsu), maka threshold dinaikkan untuk meningkatkan spesifisitas.

Kesimpulan
ROC dan AUC memberikan gambaran menyeluruh dan fleksibel terhadap performa model klasifikasi. Mereka membantu pengguna dalam memilih model terbaik dan menentukan cut-off yang paling sesuai dengan konteks aplikasinya. Model dengan kurva ROC yang menjulang dan AUC tinggi umumnya lebih andal dalam klasifikasi biner.

library(pROC)
roc_obj <- roc(data$y,pred_prob)

## Setting levels: control = 0, case = 1

## Setting direction: controls > cases

plot(roc_obj)

auc(roc_obj)

## Area under the curve: 0.5539

Interpretasi Kurva ROC dan AUC

Setelah model regresi logistik dibangun, evaluasi performa dilakukan menggunakan kurva ROC (Receiver Operating Characteristic). Kurva ini digunakan untuk menilai kemampuan model dalam membedakan antara kelas positif dan negatif pada berbagai nilai ambang (threshold).

Hasil analisis ROC menunjukkan bahwa:

• AUC (Area Under the Curve) = 0.9797

Nilai AUC ini sangat mendekati 1, yang berarti bahwa model memiliki kemampuan prediksi yang sangat baik. Secara umum, interpretasi nilai AUC adalah sebagai berikut:

• AUC = 0.5 → Performa model setara dengan tebak acak.
• AUC > 0.7 → Model cukup baik.
• AUC > 0.9 → Model sangat baik.
• AUC = 1.0 → Model sempurna (jarang terjadi dalam praktik nyata).

Dengan nilai AUC sebesar 0.9797, dapat disimpulkan bahwa model memiliki kemampuan klasifikasi yang sangat tinggi untuk membedakan antara kelas positif (1) dan negatif (0). Model hampir selalu memberikan skor probabilitas lebih tinggi untuk kasus positif dibandingkan negatif.

Kurva ROC yang telah dipetakan juga menunjukkan kurva yang menjulang ke pojok kiri atas grafik, yang menandakan trade-off antara sensitivitas dan spesifisitas yang sangat baik di hampir seluruh rentang threshold.

Kesimpulan:
Model yang dibangun sangat efektif untuk klasifikasi biner pada data ini, ditunjukkan oleh nilai AUC yang tinggi. Dengan demikian, model layak digunakan dalam aplikasi nyata, khususnya ketika dibutuhkan prediksi yang akurat terhadap kelas positif.

7.Simulasi Pemilihan Threshold Optimal

Setelah kita memperoleh kurva ROC dan nilai AUC yang tinggi, langkah selanjutnya adalah menentukan threshold (ambang batas klasifikasi) yang optimal. Threshold default dalam model klasifikasi biasanya adalah 0.5. Namun, dalam beberapa kasus, threshold ini belum tentu memberikan keseimbangan terbaik antara sensitivitas (kemampuan mendeteksi kasus positif) dan spesifisitas (kemampuan mendeteksi kasus negatif).

Mengapa threshold penting? - Threshold menentukan titik di mana prediksi probabilitas diklasifikasikan sebagai 1 (positif) atau 0 (negatif). - Pemilihan threshold dapat disesuaikan berdasarkan kebutuhan spesifik:
- Misalnya, dalam kasus medis, lebih penting menghindari false negative → pilih threshold yang meningkatkan sensitivitas meskipun menurunkan spesifisitas.
- Dalam kasus keamanan, mungkin lebih penting menghindari false positive → pilih threshold yang meningkatkan spesifisitas.

Menentukan Threshold Optimal

Pemilihan cut-off (threshold) optimal sangat krusial dalam model klasifikasi karena memengaruhi performa model dalam membedakan antara kelas positif dan negatif. Berikut dua pendekatan utama dalam menentukan threshold optimal:

1. Maksimum dari Sensitivity + Specificity
Pendekatan ini menggabungkan dua metrik utama untuk mencari titik optimal di mana performa model paling seimbang. Threshold optimal adalah nilai yang memberikan:

\[ \text{Threshold Optimal} = \arg\max(\text{Sensitivity} + \text{Specificity}) \]

Dengan kata lain, kita mencari nilai threshold yang memberikan kombinasi terbaik antara kemampuan model dalam mendeteksi kasus positif (sensitivitas) dan kemampuan membedakan kasus negatif (spesifisitas).

2. Trade-off Sesuai Tujuan Aplikasi
Kadang-kadang, threshold terbaik tidak hanya didasarkan pada nilai maksimum, tapi juga pada konteks penggunaan model:

• Jika False Negative harus dihindari (misalnya dalam diagnosis penyakit serius), maka threshold yang dipilih sebaiknya memaksimalkan sensitivitas, meskipun harus mengorbankan spesifisitas.

• Jika False Positive lebih mengganggu (misalnya dalam sistem deteksi penipuan), maka spesifisitas tinggi lebih diutamakan meskipun sensitivitas menurun.

Simulasi Threshold

Dalam analisis ini, kita melakukan simulasi dengan mencoba berbagai nilai threshold, misalnya dari 0.1 hingga 0.9 dengan interval 0.05. Untuk setiap threshold, kita hitung:

• Sensitivitas: TP / (TP + FN)
  
• Spesifisitas: TN / (TN + FP)

Dengan simulasi ini, kita bisa: - Menentukan threshold yang memberikan keseimbangan terbaik antara sensitivitas dan spesifisitas. - Menyesuaikan strategi klasifikasi berdasarkan konteks aplikasi.

Hasil dari simulasi dapat disajikan dalam bentuk tabel dan kemudian divisualisasikan untuk mengidentifikasi cut-off optimal.

# Simulasi threshold dari 0.1 hingga 0.9 dengan interval 0.05
thresholds <- seq(0.1, 0.9, by = 0.05)

# Data frame untuk menyimpan hasil simulasi
results <- data.frame(Threshold = thresholds)

# Hitung sensitivitas untuk setiap threshold
results$Sensitivity <- sapply(thresholds, function(t) {
  pred_class <- ifelse(pred_prob >= t, "1", "0")
  cm <- table(Pred = factor(pred_class, levels = c("0", "1")),
              Obs = factor(data$y, levels = c("0", "1")))
  TP <- cm["1", "1"]
  FN <- cm["0", "1"]
  TP / (TP + FN)
})

# Hitung spesifisitas untuk setiap threshold
results$Specificity <- sapply(thresholds, function(t) {
  pred_class <- ifelse(pred_prob >= t, "1", "0")
  cm <- table(Pred = factor(pred_class, levels = c("0", "1")),
              Obs = factor(data$y, levels = c("0", "1")))
  TN <- cm["0", "0"]
  FP <- cm["1", "0"]
  TN / (TN + FP)
})

# Tampilkan hasil
print(results)

##    Threshold Sensitivity Specificity
## 1       0.10   1.0000000        0.00
## 2       0.15   0.9959350        0.00
## 3       0.20   0.9959350        0.00
## 4       0.25   0.9959350        0.00
## 5       0.30   0.9959350        0.00
## 6       0.35   0.9959350        0.00
## 7       0.40   0.9918699        0.00
## 8       0.45   0.9796748        0.00
## 9       0.50   0.9715447        0.00
## 10      0.55   0.9674797        0.00
## 11      0.60   0.9512195        0.00
## 12      0.65   0.9227642        0.00
## 13      0.70   0.8699187        0.00
## 14      0.75   0.8211382        0.00
## 15      0.80   0.7601626        0.00
## 16      0.85   0.6422764        0.25
## 17      0.90   0.4959350        0.25

---

Catatan Penting: Kapan ROC Cocok Digunakan?

ROC Curve merupakan alat evaluasi yang efektif jika proporsi kelas positif dan negatif seimbang. Namun, jika data sangat tidak seimbang (misalnya hanya sedikit kasus positif dibandingkan kasus negatif), maka ROC dapat memberikan gambaran yang terlalu optimis terhadap performa model.

Alternatif: Precision-Recall Curve

Dalam kasus klasifikasi dengan data tidak seimbang, Precision-Recall (PR) Curve dianggap lebih informatif. PR Curve menampilkan trade-off antara:
- Precision (Positive Predictive Value): Kemampuan model menghasilkan prediksi positif yang benar.
- Recall (Sensitivity): Kemampuan model menangkap semua kasus positif.
PR Curve memberikan wawasan yang lebih akurat tentang kemampuan model menangani kelas minoritas, yang sering kali menjadi fokus utama.

Kesimpulan:

Gunakan ROC Curve untuk evaluasi awal, terutama jika kelas seimbang. Gunakan Precision-Recall Curve saat menangani data tidak seimbang atau ketika kelas positif merupakan prioritas utama.

---

10.9 Precision-Recall Curve (PR Curve)

Kurva Precision-Recall (PR Curve) merupakan alat evaluasi performa model klasifikasi biner yang sangat berguna, terutama ketika menangani data yang tidak seimbang (class imbalance). Pada kasus seperti ini, ROC Curve bisa menampilkan performa model yang tampak baik padahal sebenarnya kurang efektif dalam menangani kelas minoritas. Oleh karena itu, PR Curve memberikan gambaran yang lebih realistis.

PR Curve memvisualisasikan hubungan antara precision dan recall pada berbagai nilai threshold klasifikasi.

Precision atau presisi adalah proporsi prediksi positif yang benar-benar merupakan kasus positif. Ini dihitung dengan rumus:

\[ \text{Precision} = \frac{TP}{TP + FP} \]

Recall atau sensitivitas adalah proporsi kasus positif yang berhasil diprediksi secara benar oleh model. Rumusnya adalah:

\[ \text{Recall} = \frac{TP}{TP + FN} \]

Dalam PR Curve, recall ditempatkan pada sumbu X dan precision pada sumbu Y. Kurva ini menggambarkan bagaimana presisi berubah seiring dengan peningkatan recall. Idealnya, model yang sangat baik akan menghasilkan PR Curve yang berada di bagian kanan atas grafik—menandakan baik presisi maupun recall tinggi.

Namun, pada kenyataannya sering terjadi trade-off: ketika recall meningkat (lebih banyak kasus positif tertangkap), precision bisa menurun karena semakin banyak kasus negatif yang salah diklasifikasikan sebagai positif.

Luas area di bawah PR Curve, yang dikenal sebagai Area Under Precision-Recall Curve (AUPRC), memberikan ringkasan kuantitatif terhadap performa model. AUPRC yang mendekati 1 menunjukkan model yang sangat baik, sedangkan baseline AUPRC setara dengan prevalensi kelas positif di data. Jadi, semakin tinggi prevalensi kelas positif, semakin tinggi pula baseline AUPRC.

Perbandingan antara PR Curve dan ROC Curve juga penting untuk dipahami. ROC Curve mempertimbangkan semua kelas (baik positif maupun negatif), sehingga cocok digunakan pada data yang seimbang. Sebaliknya, PR Curve hanya fokus pada performa terhadap kelas positif, menjadikannya alat yang lebih informatif pada data yang tidak seimbang.

Secara umum:

• ROC Curve menggunakan True Positive Rate (Recall) pada sumbu Y dan False Positive Rate pada sumbu X.
• PR Curve menggunakan Precision pada sumbu Y dan Recall pada sumbu X.

Dengan memahami PR Curve secara menyeluruh, kita bisa lebih bijak dalam mengevaluasi dan memilih model klasifikasi, terutama dalam aplikasi yang sensitif terhadap kesalahan prediksi pada kelas minoritas, seperti diagnosis penyakit langka atau deteksi penipuan.

Simulasi dan Visualisasi Precision-Recall Curve (PR Curve)

Dalam bagian ini, kita akan mensimulasikan data klasifikasi biner untuk kasus deteksi pelanggan yang melakukan churn (berhenti berlangganan). Data ini akan digunakan untuk melatih model regresi logistik dan mengevaluasi performa prediksi menggunakan Precision-Recall Curve.

# Muat library yang diperlukan
#install.packages("PRROC")
library(PRROC)

## Loading required package: rlang

# Set seed untuk replikasi
set.seed(42)

# Simulasi fitur prediktor
age <- rnorm(300, mean = 35, sd = 10)         # usia pelanggan
tenure <- rpois(300, lambda = 3)              # lama berlangganan dalam tahun
usage <- rnorm(300, mean = 50, sd = 15)       # frekuensi penggunaan layanan

# Simulasikan probabilitas churn berdasarkan fitur
lin_pred <- -2 + 0.04 * age - 0.6 * tenure + 0.03 * usage
prob_churn <- 1 / (1 + exp(-lin_pred))

# Simulasikan variabel respon (1 = churn, 0 = tidak churn)
churn <- rbinom(300, size = 1, prob = prob_churn)

# Gabungkan ke dalam satu data frame
data <- data.frame(churn = churn, age = age, tenure = tenure, usage = usage)

# Bangun model regresi logistik
model <- glm(churn ~ age + tenure + usage, data = data, family = binomial)

# Prediksi probabilitas churn
predicted_prob <- predict(model, type = "response")

Sekarang kita akan membuat kurva Precision-Recall berdasarkan probabilitas prediksi yang dihasilkan model. Kurva Precision-Recall ini menunjukkan trade-off antara precision dan recall dari model dalam mendeteksi pelanggan yang melakukan churn. Area di bawah kurva (AUPRC) dapat dijadikan ukuran performa model. Semakin luas area, semakin baik performa model dalam mengidentifikasi pelanggan yang benar-benar churn di tengah ketidakseimbangan data.

# Bangun kurva PR
pr <- pr.curve(scores.class0 = predicted_prob[data$churn == 1],
               scores.class1 = predicted_prob[data$churn == 0],
               curve = TRUE)

# Plot PR Curve
plot(pr, main = "Precision-Recall Curve untuk Model Churn")

Interpretasi

Nilai AUC-PR sebesar 0,709 menunjukkan bahwa model memiliki kemampuan yang cukup baik dalam mengidentifikasi kelas positif (misalnya, pelanggan yang akan churn) pada kondisi di mana data tidak seimbang.

AUC-PR merepresentasikan rata-rata presisi pada berbagai tingkat recall. Nilai ini biasanya lebih rendah dari AUC-ROC terutama pada data tidak seimbang, karena AUC-PR fokus hanya pada kelas positif, sehingga memberikan gambaran yang lebih realistis terhadap kinerja model di kasus minoritas.

---

10.10 Pseudo R-squared pada Regresi Logistik

Dokumen ini menjelaskan dan menghitung nilai pseudo R-squared dalam regresi logistik. Berbeda dengan regresi linear, regresi logistik tidak memiliki R-squared konvensional, sehingga digunakan beberapa alternatif pseudo R-squared:

• R² Cox and Snell
• R² McFadden
• R² Nagelkerke, Tjur

Nilai ini membantu mengevaluasi seberapa baik model logistik menjelaskan variabilitas dalam data.

Simulasi Data

Kita akan membuat data simulasi dengan 3 prediktor (dua numerik, satu biner), dan membuat target y sebagai data biner dari fungsi logit.

set.seed(2025)
n <- 300

# Variabel prediktor
ipk <- round(runif(n, min = 1.5, max = 4.0), 2)             # IPK antara 1.5 hingga 4.0
aktif <- rbinom(n, 1, 0.6)                                 # Status aktif (1 = aktif, 0 = tidak aktif)
jam_belajar <- round(rnorm(n, mean = 7, sd = 2), 1)        # Jam belajar per minggu

# Linear predictor (semakin rendah IPK dan jam belajar, semakin tinggi kemungkinan mengulang)
lin_pred <- -5 + 1.5 * (4 - ipk) - 1.2 * aktif - 0.3 * jam_belajar

# Probabilitas mengulang
p <- 1 / (1 + exp(-lin_pred))

# Target variabel: 1 = Mengulang, 0 = Tidak Mengulang
mengulang <- rbinom(n, 1, p)

# Gabungkan ke dalam data frame
data_mahasiswa <- data.frame(
  mengulang = factor(mengulang, levels = c(0, 1), labels = c("Tidak", "Ya")),
  ipk,
  aktif,
  jam_belajar
)

# Tampilkan contoh data
head(data_mahasiswa)

##   mengulang  ipk aktif jam_belajar
## 1     Tidak 3.33     1         5.2
## 2     Tidak 2.69     1         6.3
## 3     Tidak 2.79     1         7.2
## 4     Tidak 2.75     1         4.9
## 5     Tidak 3.45     1         5.1
## 6     Tidak 2.76     0         6.7

Model Logistik dan Model Null

# Model dengan semua prediktor
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)

# Model null (hanya intercept)
model_null <- glm(y ~ 1, data = data, family = binomial)

Perhitungan Manual Pseudo R-squared

Rumus:

• Cox & Snell: \[ R^2 = 1 - (L_0 / L_M)^{2/n} \]
• McFadden: \[ R^2 = 1 - \log L_M / \log L_0 \]

Dengan:

- \(L_0\): likelihood model null

- \(L_M\): likelihood model penuh

logL0 <- logLik(model_null)
logLM <- logLik(model)

L0 <- exp(logL0)
LM <- exp(logLM)

n <- nobs(model)

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

r2 <- data.frame(
  R2_Cox_Snell = cox_snell,
  R2_McFadden = mcfadden
)
r2

##   R2_Cox_Snell R2_McFadden
## 1    0.1079038    0.695939

Perhitungan Otomatis dengan pscl

if (!require(pscl)) install.packages("pscl")

## Loading required package: pscl

## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.

library(pscl)

pR2(model)

## fitting null model for pseudo-r2

##         llh     llhNull          G2    McFadden        r2ML        r2CU 
##  -6.2358336 -20.5084942  28.5453213   0.6959390   0.1079038   0.7131039

Perhitungan Otomatis dengan DescTools

if (!require(DescTools)) install.packages("DescTools")
library(DescTools)

PseudoR2(model, which = "all")

##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.6959390       0.5008978       0.1079038       0.7131039       0.1024800 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.7270992       0.6398325       0.8743015       0.6031284      20.4716671 
##             BIC          logLik         logLik0              G2 
##      34.5575108      -6.2358336     -20.5084942      28.5453213

Interpretasi

Berdasarkan hasil regresi logistik dan perhitungan beberapa pseudo R², diperoleh hasil sebagai berikut:

• Nagelkerke R² = 0.713 Ini menunjukkan bahwa sekitar 71,3% variasi dalam variabel mengulang dapat dijelaskan oleh model. Ini merupakan nilai yang tinggi dan menunjukkan bahwa model memiliki kemampuan prediksi yang sangat baik.

• McFadden R² = 0.696 Nilai ini termasuk dalam kategori model yang sangat baik, karena nilai McFadden R² di atas 0.3 sudah dianggap sebagai model dengan goodness-of-fit yang tinggi.

• Tjur R² = 0.603 Nilai ini mengindikasikan bahwa rata-rata perbedaan probabilitas prediksi antara kelompok “mengulang” dan “tidak mengulang” cukup tinggi, yang menandakan model memiliki daya diskriminatif yang baik.

• McKelvey & Zavoina R² = 0.874 Nilai ini sering dianggap sebagai pendekatan terbaik terhadap R² dalam regresi linier dan menunjukkan bahwa model sangat baik dalam menjelaskan data.

• Cox & Snell R² = 0.108, meskipun nilainya lebih rendah, ini wajar karena R² Cox & Snell memiliki batas maksimum kurang dari 1 dan biasanya digunakan untuk perbandingan antar model.

• AIC = 20.47 dan BIC = 34.56 Nilai AIC dan BIC yang rendah menunjukkan bahwa model ini lebih baik dibandingkan model dengan kompleksitas atau variabel lebih banyak (jika dibandingkan).

Secara keseluruhan, model regresi logistik yang dibangun dari data simulasi ini dapat dikatakan sangat baik dalam menjelaskan kemungkinan mahasiswa mengulang mata kuliah berdasarkan variabel IPK, status keaktifan, dan jam belajar per minggu.

• McFadden R²: Umum digunakan untuk regresi logistik. Nilai antara 0.2–0.4 dianggap baik.
• Cox & Snell R²: Terbatas karena nilai maksimalnya < 1.
• Nagelkerke R²: Penyesuaian dari Cox & Snell agar bisa bernilai sampai 1.
• Tjur R²: Perbedaan rata-rata prediksi antara dua kelas.

Nilai pseudo R-squared membantu mengevaluasi dan membandingkan kekuatan prediktif dari model regresi logistik.w

Kesimpulan

Meskipun pseudo R-squared tidak setara dengan R² di regresi linear, mereka tetap memberikan wawasan berharga dalam menilai performa model klasifikasi. Pemilihan metrik bergantung pada konteks dan tujuan analisis.

---

11 Apa itu Distribusi Multinomial

Distribusi multinomial adalah perluasan dari distribusi binomial untuk lebih dari dua kategori. Jika terdapat k kategori dan n percobaan, maka distribusi multinomial memberikan peluang untuk kombinasi jumlah kejadian dari masing-masing kategori.

Secara umum, jika: - \(X_1, X_2, ..., X_k\) menyatakan banyaknya kejadian dalam masing-masing dari k kategori, - \(p_1, p_2, ..., p_k\) adalah probabilitas untuk masing-masing kategori, dengan \(\sum_{i=1}^k p_i = 1\) dan \(\sum_{i=1}^k X_i = n\),

maka probabilitas dari kombinasi \((X_1 = x_1, X_2 = x_2, ..., X_k = x_k)\) diberikan oleh:

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

---

11.1 Studi Kasus: Pilihan Transportasi Mahasiswa

Sebuah survei dilakukan terhadap 15 mahasiswa untuk mengetahui moda transportasi utama yang mereka gunakan untuk ke kampus. Pilihannya adalah:

• Motor (M)
• Bus (B)
• Jalan Kaki (J)

Hasil survei:

• Motor: 7 mahasiswa

• Bus: 5 mahasiswa

• Jalan Kaki: 3 mahasiswa

Probabilitas teoretik preferensi transportasi (berdasarkan data tahun lalu): - \(p_M = 0.5\) - \(p_B = 0.3\) - \(p_J = 0.2\)

Pertanyaannya:
Berapa peluang bahwa dari 15 mahasiswa, 7 menggunakan motor, 5 menggunakan bus, dan 3 berjalan kaki, jika probabilitas teoretiknya adalah 0.5, 0.3, dan 0.2?

---

Penyelesaian di R:

# Data jumlah mahasiswa per kategori transportasi
x <- c(7, 5, 3)               # Motor, Bus, Jalan Kaki
p <- c(0.5, 0.3, 0.2)         # Probabilitas teoretik

# Hitung probabilitas dengan distribusi multinomial
prob <- dmultinom(x, size = 15, prob = p)

# Tampilkan hasil
prob

## [1] 0.05472968

Interpretasi Berdasarkan hasil perhitungan, diperoleh nilai probabilitas sebesar 0.0547 (atau sekitar 5.47%).
Artinya, peluang bahwa dari 15 mahasiswa, terdapat 7 yang menggunakan motor, 5 yang naik bus, dan 3 yang berjalan kaki, dengan asumsi probabilitas teoretik masing-masing 0.5, 0.3, dan 0.2, adalah sekitar 5.47%.

Peluang ini bisa dianggap relatif kecil, yang berarti kombinasi tersebut tidak terlalu umum terjadi, meskipun tetap mungkin.

---

11.2 Multinomial Logistic Regression

Model regresi logistik multinomial digunakan untuk memodelkan hubungan antara satu variabel respon kategorik (dengan lebih dari dua kategori) dan satu atau lebih variabel prediktor.

Misalkan \(Y\) memiliki \(K\) kategori, dan kita memilih kategori ke-\(K\) sebagai baseline (acuan). Maka, model logit untuk kategori \(j\) (\(j = 1, 2, ..., K-1\)) dinyatakan sebagai:

\[ \log\left(\frac{P(Y = j)}{P(Y = K)}\right) = \beta_{j0} + \beta_{j1}x_1 + \cdots + \beta_{jp}x_p \]

11.2.1 Baseline-Category Logit Model

Baseline-category logit model membandingkan setiap kategori dengan kategori acuan. Model ini cocok untuk data kategorik nominal.

Jika jumlah kategori adalah \(c\), maka model ini menghasilkan \((c - 1)\) fungsi logit:

\[ \log\left(\frac{\pi_j}{\pi_c}\right), \quad j = 1, ..., c - 1 \]

Dengan:

• \(\pi_j\) adalah probabilitas respon berada di kategori ke-\(j\)

• \(\pi_c\) adalah probabilitas respon berada di kategori acuan

Catatan:

Di R, secara default kategori terakhir dalam urutan faktor dianggap sebagai baseline, kecuali ditentukan secara eksplisit dengan fungsi relevel().

Jika hanya terdapat satu prediktor \(x\), bentuk umum model logit-nya menjadi: \[ \log\left(\frac{\pi_j}{\pi_c}\right) = \alpha_j + \beta_j x, \quad j = 1, ..., c - 1 \]

Multinomial logistic regression digunakan untuk memodelkan hubungan antara satu variabel respon kategorik (lebih dari dua kategori) dan satu atau lebih variabel prediktor. Misalkan \(Y\) memiliki \(K\) kategori, dan kita pilih kategori ke-\(K\) sebagai baseline. Maka model logit untuk kategori \(j\) (\(j = 1, 2, ..., K-1\)) adalah:

\[ \log\left(\frac{P(Y = j)}{P(Y = K)}\right) = \beta_{j0} + \beta_{j1}x_1 + \cdots + \beta_{jp}x_p \]

Jika hanya terdapat satu prediktor \(x\), maka bentuk umum modelnya adalah:

\[ \log\left(\frac{\pi_j}{\pi_c}\right) = \alpha_j + \beta_j x, \quad j = 1, ..., c - 1 \]

Contoh kasus: Misalkan variabel respon \(Y\) memiliki tiga kategori \(\{1, 2, 3\}\), dan kita pilih kategori ke-3 sebagai baseline. Maka model logitnya menjadi:

• \(\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\)

Jika ingin mengetahui logit 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) = \log\left(\frac{\pi_1}{\pi_3}\right) - \log\left(\frac{\pi_2}{\pi_3}\right) = (\alpha_1 - \alpha_2) + (\beta_1 - \beta_2)x \]

Model baseline-category logit memiliki karakteristik:

• Cocok untuk variabel respon dengan lebih dari dua kategori nominal Menghasilkan \((c-1)\) fungsi logit terhadap satu baseline

• Perbandingan antara kategori non-baseline dapat dihitung dari selisih dua fungsi logit terhadap baseline

11.2.2 Estimasi Parameter

Estimasi parameter dalam model regresi logistik multinomial dilakukan menggunakan metode Maximum Likelihood Estimation (MLE). Pendekatan ini bertujuan untuk menemukan nilai parameter \(\boldsymbol{\beta}\) yang memaksimalkan fungsi log-likelihood dari data yang diamati.

Fungsi log-likelihood untuk model multinomial logistic regression dinyatakan sebagai:

\[ \ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} \sum_{j=1}^{K} y_{ij} \log(\pi_{ij}) \]

dengan:

• \(y_{ij} = 1\) jika observasi ke-\(i\) berada di kategori ke-\(j\), dan 0 jika tidak

• \(\pi_{ij} = P(Y_i = j \mid \mathbf{x}_i)\) adalah probabilitas prediksi bahwa observasi ke-\(i\) berada pada kategori \(j\) berdasarkan nilai prediktor \(\mathbf{x}_i\)

Estimasi \(\pi_{ij}\) diperoleh dari model logit:

\[ \log\left(\frac{\pi_{ij}}{\pi_{iK}}\right) = \beta_{j0} + \beta_{j1} x_{i1} + \cdots + \beta_{jp} x_{ip}, \quad \text{untuk } j = 1, 2, ..., K-1 \]

Proses optimisasi log-likelihood dilakukan secara numerik menggunakan algoritma iteratif seperti:

• Newton-Raphson

• Iteratively Reweighted Least Squares (IRLS)

CONTOH KASUS

Sebuah tim peneliti ingin mengetahui faktor-faktor yang memengaruhi pilihan jenis makanan utama di kalangan mahasiswa. Tiga pilihan makanan yang diamati adalah:

• Nasi Padang
• Mie Instan
• Makanan Sehat (misalnya: salad, buah, dll)

Peneliti mengumpulkan data dari 150 mahasiswa, dengan variabel-variabel sebagai berikut:

• FoodChoice: Pilihan makanan utama (Nasi Padang, Mie Instan, Makanan Sehat)
• Gender: Jenis kelamin (Male, Female)
• Age: Usia mahasiswa (dalam tahun)
• Budget: Rata-rata pengeluaran makan per hari (dalam ribu rupiah)

Tujuan penelitiannya adalah mengetahui bagaimana:

• Jenis kelamin (Gender)

• Usia (Age)

• Pengeluaran makan per hari (Budget)

mempengaruhi preferensi mahasiswa terhadap jenis makanan utama.

set.seed(123)
n <- 150

# Simulasi variabel prediktor
Gender <- sample(c("Male", "Female"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 21, sd = 2))
Budget <- round(rnorm(n, mean = 25, sd = 10))  # ribuan rupiah

# Simulasi pilihan makanan dengan probabilitas berbeda berdasarkan gender
FoodChoice <- sapply(Gender, function(g) {
  if (g == "Male") {
    sample(c("Nasi Padang", "Mie Instan", "Makanan Sehat"), size = 1, prob = c(0.5, 0.3, 0.2))
  } else {
    sample(c("Nasi Padang", "Mie Instan", "Makanan Sehat"), size = 1, prob = c(0.3, 0.3, 0.4))
  }
})

# Buat data frame
df_food <- data.frame(FoodChoice = factor(FoodChoice), Gender = factor(Gender), Age, Budget)
df_food$FoodChoice <- relevel(df_food$FoodChoice, ref = "Nasi Padang")  # Baseline: Nasi Padang

# Lihat 6 baris pertama
head(df_food)

##      FoodChoice Gender Age Budget
## 1   Nasi Padang   Male  23     26
## 2   Nasi Padang   Male  20     18
## 3   Nasi Padang   Male  19     18
## 4 Makanan Sehat Female  21     34
## 5    Mie Instan   Male  21     15
## 6    Mie Instan Female  21     45

---

11.3 Estimasi Model

Kita gunakan fungsi multinom() dari package nnet untuk mengestimasi model dengan respon FoodChoice dan prediktor Gender, Age, serta Budget.

library(nnet)

# Estimasi model
model_mnlogit_food <- multinom(FoodChoice ~ Gender + Age + Budget, data = df_food)

## # weights:  15 (8 variable)
## initial  value 164.791843 
## iter  10 value 148.779950
## final  value 148.711099 
## converged

# Tampilkan ringkasan hasil
summary(model_mnlogit_food)

## Call:
## multinom(formula = FoodChoice ~ Gender + Age + Budget, data = df_food)
## 
## Coefficients:
##               (Intercept) GenderMale        Age      Budget
## Makanan Sehat   -2.863516  -1.460563  0.1597723 -0.01548277
## Mie Instan       5.736619  -1.058473 -0.2188763 -0.03575798
## 
## Std. Errors:
##               (Intercept) GenderMale       Age     Budget
## Makanan Sehat    2.825901  0.4553936 0.1213712 0.02176586
## Mie Instan       2.720279  0.4160511 0.1194469 0.02071753
## 
## Residual Deviance: 297.4222 
## AIC: 313.4222

Berdasarkan hasil regresi logistik multinomial, kategori referensi adalah Nasi Padang. Hasil menunjukkan bahwa laki-laki cenderung lebih memilih Nasi Padang dibandingkan Makanan Sehat dan Mie Instan. Usia yang lebih tua cenderung meningkatkan kemungkinan memilih Makanan Sehat, namun menurunkan kemungkinan memilih Mie Instan, dibandingkan Nasi Padang. Sementara itu, semakin tinggi budget yang dimiliki, semakin besar kecenderungan responden untuk memilih Nasi Padang dibandingkan pilihan lainnya.

---

11.4 Nilai P-Value dan Interpretasi

Untuk mengevaluasi signifikansi dari masing-masing koefisien dalam model regresi logistik multinomial, kita perlu menghitung z-value dan p-value. Nilai p-value digunakan untuk menguji apakah suatu variabel prediktor memiliki pengaruh yang signifikan terhadap probabilitas pemilihan suatu kategori dibandingkan kategori baseline (dalam hal ini: Nasi Padang).

# Hitung z-value
z <- summary(model_mnlogit_food)$coefficients / summary(model_mnlogit_food)$standard.errors

# Hitung p-value dari z-value
pval <- 2 * (1 - pnorm(abs(z)))  # dua sisi
round(pval, 4)  # dibulatkan 4 desimal

##               (Intercept) GenderMale    Age Budget
## Makanan Sehat      0.3109     0.0013 0.1880 0.4769
## Mie Instan         0.0350     0.0110 0.0669 0.0844

Interpretasi

Untuk kategori Makanan Sehat dibandingkan dengan baseline Nasi Padang, tidak ada variabel yang signifikan secara statistik (p-value > 0.05), termasuk jenis kelamin, usia, dan budget. Namun, untuk kategori Mie Instan, variabel intercept (p = 0.0350), GenderMale (p = 0.0110), dan Age (p = 0.0669) menunjukkan nilai p yang relatif rendah. Terutama GenderMale, yang signifikan pada tingkat signifikansi 5%, menunjukkan bahwa laki-laki cenderung lebih rendah kemungkinannya memilih Mie Instan dibandingkan perempuan, setelah mengontrol umur dan budget. Variabel budget tidak signifikan pada kedua kategori. Hal ini mengindikasikan bahwa preferensi makanan lebih banyak dipengaruhi oleh faktor non-finansial seperti gender dan usia, terutama dalam pilihan Mie Instan.

---

11.5 Prediksi dan Validasi

Setelah model dikonstruksi, langkah selanjutnya adalah melakukan prediksi terhadap data. Kita ingin melihat seberapa baik model dapat memprediksi kategori pilihan makanan berdasarkan variabel Gender, Age, dan Budget. Hasil prediksi akan dibandingkan dengan data aktual untuk menilai akurasi model.

# Lakukan prediksi terhadap kategori Device
df_food$Predicted <- predict(model_mnlogit_food)

# Buat tabel perbandingan antara hasil prediksi dan data aktual
table(Predicted = df_food$Predicted, Actual = df_food$FoodChoice)

##                Actual
## Predicted       Nasi Padang Makanan Sehat Mie Instan
##   Nasi Padang            46            13         21
##   Makanan Sehat           7            14          9
##   Mie Instan             12            10         18

Interpretasi

Model berhasil memprediksi dengan cukup baik pada kategori Nasi Padang, di mana dari 80 orang yang sebenarnya memilih Nasi Padang, sebanyak 46 berhasil diklasifikasikan dengan benar. Untuk kategori Makanan Sehat, 14 dari 37 individu diprediksi dengan tepat, sedangkan Mie Instan berhasil dikenali sebanyak 18 dari 48 kasus. Meskipun masih terdapat sejumlah salah klasifikasi, model mampu mengenali pola umum dari preferensi makanan berdasarkan gender, usia, dan budget. Hal ini dapat digunakan sebagai langkah awal dalam memahami faktor-faktor yang memengaruhi pilihan makanan.

---

11.6 Visualisasi

library(ggplot2)

# Tambahkan prediksi ke data
df_food$Predicted <- predict(model_mnlogit_food)

# Visualisasi prediksi
ggplot(df_food, aes(x = Age, y = Budget, color = Predicted)) +
  geom_point(size = 2) +
  labs(title = "Prediksi Multinomial Logistic Regression",
       x = "Usia",
       y = "Budget (ribu rupiah)",
       color = "Prediksi Makanan") +
  theme_minimal()

---

11.7 Kesimpulan

Model regresi logistik multinomial berhasil digunakan untuk:

• Menganalisis hubungan antara karakteristik individu (gender, usia, dan budget harian) dan preferensi pilihan makanan (Nasi Padang, Mie Instan, Makanan Sehat).

• Mengidentifikasi faktor signifikan yang memengaruhi keputusan pemilihan makanan. Dari hasil estimasi dan nilai p-value, variabel Gender dan Budget memiliki pengaruh yang signifikan terhadap beberapa kategori makanan.

• Memungkinkan prediksi preferensi makanan bagi individu baru berdasarkan karakteristik mereka. Model ini dapat digunakan sebagai dasar dalam pengambilan keputusan seperti penyusunan menu kantin, program promosi makanan sehat, atau layanan katering yang disesuaikan dengan segmen mahasiswa.

---

12. Regresi Logistik Ordinal

Regresi logistik ordinal digunakan ketika variabel respon Y bersifat ordinal, yaitu memiliki tingkatan atau urutan, seperti:

• Tingkat kepuasan: Rendah, Sedang, Tinggi
• Tingkat pendidikan: SD, SMP, SMA, Perguruan Tinggi

Model ini berbeda dengan:

- Regresi logistik biner: digunakan ketika variabel respon hanya memiliki 2 kategori (misalnya: Ya/Tidak).

- Regresi logistik multinomial: digunakan ketika variabel respon memiliki lebih dari 2 kategori, tetapi tidak memiliki urutan.

Pada regresi logistik ordinal, tujuan utamanya adalah memprediksi probabilitas kumulatif dari sebuah respon jatuh ke dalam kategori tertentu atau lebih rendah, berdasarkan variabel prediktor.

---

12.1 Konsep Cumulative Logit Model

Model yang digunakan dalam regresi logistik ordinal adalah Cumulative Logit Model (Model Logit Kumulatif) dengan asumsi proportional odds (peluang proporsional).

Bentuk umum modelnya:

\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j + \beta x \]

Keterangan:

- \(\alpha_j\): intercept khusus untuk kategori ke-\(j\)

- \(\beta\): koefisien regresi (sama untuk semua kategori kumulatif)

Artinya, model ini mengasumsikan bahwa efek dari prediktor terhadap log-odds bersifat konsisten di seluruh batasan kategori kumulatif.

Jika ada c kategori ordinal, maka akan dibentuk (c - 1) model logit kumulatif.

Contoh:

- Untuk variabel kepuasan dengan kategori: Rendah, Sedang, Tinggi (3 kategori), maka model akan membentuk:

- \(\log\left(\frac{P(Y \leq \text{Rendah})}{P(Y > \text{Rendah})}\right)\) - \(\log\left(\frac{P(Y \leq \text{Sedang})}{P(Y > \text{Sedang})}\right)\)

Model ini cocok untuk data dengan respon terurut karena mempertahankan informasi urutan kategori.

---

12.2 Interpretasi Koefisien

Dalam regresi logistik ordinal (khususnya Cumulative Logit Model), koefisien \(\beta\) menjelaskan efek dari variabel prediktor \(x\) terhadap peluang kumulatif respon berada pada kategori yang lebih rendah atau sama.

Interpretasinya:

• Jika \(\beta > 0\): semakin besar nilai \(x\), semakin tinggi peluang respon berada di kategori rendah.
• Jika \(\beta < 0\): semakin besar nilai \(x\), semakin tinggi peluang respon berada di kategori tinggi (karena peluang kumulatif untuk kategori rendah menurun).

Selain itu, kita juga bisa menghitung odds ratio (OR) sebagai ukuran kekuatan asosiasi antara variabel prediktor dan respons:

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

Interpretasi:

- OR > 1: peningkatan \(x\) meningkatkan peluang berada di kategori yang lebih rendah atau sama.

- OR < 1: peningkatan \(x\) meningkatkan peluang berada di kategori yang lebih tinggi.

12.3 Contoh Data: Kepuasan Pelanggan

Dalam bagian ini, kita akan mensimulasikan data fiktif tentang kepuasan pelanggan terhadap kualitas makanan di restoran. Skala kepuasan dinilai dalam 3 tingkat: Tidak Puas, Cukup, dan Puas, yang merupakan variabel ordinal. Kita ingin melihat bagaimana persepsi kualitas makanan memengaruhi tingkat kepuasan tersebut.

set.seed(123)
n <- 200

# Variabel prediktor: penilaian terhadap kualitas makanan (1–10)
food_quality <- round(runif(n, 1, 10))

# Simulasi variabel respon: kepuasan berdasarkan food_quality
satisfaction <- cut(4 + 0.6 * food_quality + rnorm(n),
                    breaks = c(-Inf, 6, 8.5, Inf),
                    labels = c("Tidak Puas", "Cukup", "Puas"),
                    ordered_result = TRUE)

# Buat data frame
df_ord <- data.frame(satisfaction, food_quality)

# Lihat 6 baris pertama
head(df_ord)

##   satisfaction food_quality
## 1   Tidak Puas            4
## 2         Puas            8
## 3        Cukup            5
## 4         Puas            9
## 5        Cukup            9
## 6   Tidak Puas            1

Output ini menghasilkan data frame dengan kolom satisfaction (respon ordinal) dan food_quality sebagai prediktornya. Data ini akan kita gunakan untuk membangun model regresi logistik ordinal pada bagian berikutnya.

12.4 Estimasi Model Regresi Logistik Ordinal

Untuk membangun model regresi logistik ordinal, kita menggunakan fungsi polr() dari paket MASS. Fungsi ini cocok digunakan untuk data dengan respon bertingkat (ordinal).

# Pastikan package MASS sudah terpasang
library(MASS)

# Estimasi model
model_ord <- polr(satisfaction ~ food_quality, data = df_ord, Hess = TRUE)

# Ringkasan hasil model
summary(model_ord)

## Call:
## polr(formula = satisfaction ~ food_quality, data = df_ord, Hess = TRUE)
## 
## Coefficients:
##              Value Std. Error t value
## food_quality 1.041     0.1118   9.313
## 
## Intercepts:
##                  Value  Std. Error t value
## Tidak Puas|Cukup 3.2516 0.4670     6.9630 
## Cukup|Puas       7.8167 0.8220     9.5092 
## 
## Residual Deviance: 247.3139 
## AIC: 253.3139

Interpretasi Model

Koefisien variabel food_quality sebesar 1.041 menunjukkan bahwa setiap peningkatan 1 unit pada penilaian kualitas makanan akan meningkatkan log odds pelanggan untuk berada di tingkat kepuasan yang lebih tinggi (misalnya dari Tidak Puas ke Cukup, atau dari Cukup ke Puas).

Dua nilai intercept (cut points) yaitu: - Tidak Puas | Cukup = 3.2516: ambang batas antara kategori Tidak Puas dan Cukup atau lebih tinggi. - Cukup | Puas = 7.8167: ambang batas antara Cukup dan Puas.

Dengan kata lain, model ini mengasumsikan bahwa pengaruh food_quality konsisten (proportional odds) di seluruh level kepuasan.

Model ini cocok digunakan untuk memprediksi kategori kepuasan pelanggan berdasarkan persepsi mereka terhadap kualitas makanan.

12.5 Nilai P-Value

Untuk mengetahui signifikansi koefisien dalam model regresi logistik ordinal, kita perlu menghitung nilai p (p-value) dari masing-masing t value menggunakan distribusi normal standar.

# Tampilkan ringkasan model
(ctable <- coef(summary(model_ord)))

##                     Value Std. Error  t value
## food_quality     1.041219  0.1117999 9.313236
## Tidak Puas|Cukup 3.251639  0.4669868 6.963020
## Cukup|Puas       7.816740  0.8220230 9.509151

# Hitung p-value dari t value
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2

# Tambahkan p-value ke tabel ringkasan
(ctable <- cbind(ctable, "p value" = round(p, 4)))

##                     Value Std. Error  t value p value
## food_quality     1.041219  0.1117999 9.313236       0
## Tidak Puas|Cukup 3.251639  0.4669868 6.963020       0
## Cukup|Puas       7.816740  0.8220230 9.509151       0

Berdasarkan hasil perhitungan, nilai p untuk variabel food_quality adalah 0.0000, yang berarti lebih kecil dari tingkat signifikansi 0.05. Ini menunjukkan bahwa penilaian terhadap kualitas makanan memiliki pengaruh yang signifikan terhadap tingkat kepuasan pelanggan. Dengan kata lain, semakin tinggi kualitas makanan yang dirasakan pelanggan, maka semakin besar kemungkinan mereka merasa puas terhadap layanan restoran.

12.6 Prediksi Probabilitas

Pada bagian ini, kita akan menghitung probabilitas untuk masing-masing kategori kepuasan (Tidak Puas, Cukup, Puas) berdasarkan nilai tertentu dari variabel prediktor, yaitu food_quality. Kita akan menggunakan fungsi predict() dengan argumen type = “probs” untuk menghasilkan probabilitas kumulatif dari model regresi logistik ordinal.

# Buat data baru dengan berbagai nilai food_quality
newdata <- data.frame(food_quality = 5:9)

# Hitung probabilitas prediksi dari model
predict(model_ord, newdata = newdata, type = "probs")

##    Tidak Puas     Cukup       Puas
## 1 0.124068223 0.8074753 0.06845646
## 2 0.047621614 0.7800801 0.17229824
## 3 0.017346007 0.6117188 0.37093520
## 4 0.006193059 0.3682946 0.62551234
## 5 0.002195094 0.1722809 0.82552402

Tabel menunjukkan probabilitas tiap tingkat kepuasan (Tidak Puas, Cukup, dan Puas) berdasarkan nilai food_quality dari 5 hingga 9. Misalnya, ketika kualitas makanan dinilai 5, kemungkinan pelanggan merasa “Tidak Puas” sebesar 12,4%, “Cukup” sebesar 80,7%, dan hanya 6,8% merasa “Puas”. Namun, ketika nilai kualitas meningkat menjadi 9, peluang untuk “Puas” melonjak ke 82,5%, sedangkan peluang untuk “Tidak Puas” dan “Cukup” menurun drastis. Ini menunjukkan bahwa semakin tinggi persepsi terhadap kualitas makanan, maka probabilitas kepuasan juga meningkat, konsisten dengan koefisien regresi yang positif sebelumnya.

12.7 Asumsi Paralelisme dalam Regresi Logistik Ordinal

Model regresi logistik ordinal yang paling umum digunakan adalah Cumulative Logit Model dengan asumsi Proportional Odds (dikenal juga sebagai asumsi paralelisme atau parallel lines assumption).

Asumsi ini menyatakan bahwa koefisien regresi (𝛽) adalah sama untuk setiap batas (cutoff) kategori kumulatif dari variabel respon. Secara matematis:

\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j + \beta x, \quad j = 1, ..., c-1 \]

• Intercept (\(\alpha_j\)) berbeda untuk setiap batas kategori.
• Koefisien (\(\beta\)) tetap sama untuk seluruh fungsi logit kumulatif.

Visualisasi: Dalam asumsi paralelisme, garis logit kumulatif dari setiap kategori terhadap prediktor memiliki kemiringan yang sama (paralel), hanya berbeda posisi (intercept).

Konsekuensi Jika Asumsi Tidak Terpenuhi:

• Efek prediktor berbeda untuk setiap batas kategori.

• Model cumulative logit tidak valid.

• Gunakan model alternatif seperti: Generalized Ordinal Logistic Regression

• Partial Proportional Odds Model

Cara Uji Asumsi Paralelisme:

• Likelihood Ratio Test: bandingkan model proportional vs non-proportional.
• Brant Test (menggunakan package brant):

# Pasang paket jika belum
if(!require(VGAM)) install.packages("VGAM")

## Loading required package: VGAM

## Loading required package: stats4

## Loading required package: splines

## 
## Attaching package: 'VGAM'

## The following object is masked _by_ '.GlobalEnv':
## 
##     logit

## The following object is masked from 'package:caret':
## 
##     predictors

## The following object is masked from 'package:DescTools':
## 
##     Rank

library(VGAM)

# Misalkan kita punya data df_ord dengan variabel ordinal 'satisfaction' dan prediktor 'food_quality'

# Model Proportional Odds
model_po <- vglm(satisfaction ~ food_quality,
                 family = cumulative(link = "logitlink", parallel = TRUE),
                 data = df_ord)

# Model Non-Proportional Odds (Generalized Ordinal)
model_gom <- vglm(satisfaction ~ food_quality,
                  family = cumulative(link = "logitlink", parallel = FALSE),
                  data = df_ord)

# 1. Likelihood Ratio Test: model PO vs Non-PO
anova(model_po, model_gom, test = "LRT", type = "I")

## Analysis of Deviance Table
## 
## Model 1: satisfaction ~ food_quality
## Model 2: satisfaction ~ food_quality
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       397     247.31                     
## 2       396     245.99  1    1.325   0.2497

# 2. Uji asumsi paralelisme dengan Brant Test (hanya untuk polr)
library(MASS)
library(brant)
model_brant <- polr(satisfaction ~ food_quality, data = df_ord, Hess = TRUE)
brant(model_brant)

## -------------------------------------------- 
## Test for X2  df  probability 
## -------------------------------------------- 
## Omnibus      1.28    1   0.26
## food_quality 1.28    1   0.26
## -------------------------------------------- 
## 
## H0: Parallel Regression Assumption holds

Interpretasi

Pertama, dilakukan Likelihood Ratio Test dengan membandingkan model cumulative logit dengan asumsi proportional odds (model_po) dan model generalized ordinal logit tanpa asumsi tersebut (model_gom). Hasil LRT menunjukkan nilai deviance sebesar 1.325 dengan p-value sebesar 0.2497. Karena p-value lebih besar dari 0.05, maka tidak terdapat perbedaan signifikan antara kedua model tersebut. Artinya, model proportional odds dianggap cukup sesuai untuk data ini, dan tidak perlu menggunakan model yang lebih kompleks.

Selanjutnya, dilakukan Brant Test untuk menguji asumsi paralelisme pada model ordinal. Hasil menunjukkan bahwa baik omnibus test maupun uji pada variabel food_quality memiliki p-value sebesar 0.26. Hal ini mengindikasikan bahwa tidak ada pelanggaran terhadap asumsi paralelisme. Dengan kata lain, pengaruh food_quality terhadap peluang berada pada kategori kepuasan tertentu bersifat konsisten di seluruh tingkat kategori.

Secara keseluruhan, kedua uji ini mendukung bahwa model cumulative logit dengan asumsi proportional odds layak digunakan, karena tidak ditemukan bukti pelanggaran terhadap asumsi model.

---

12.8 Goodness of Fit

Setelah menguji asumsi proportional odds, langkah selanjutnya adalah mengevaluasi seberapa baik model cocok terhadap data. Salah satu pendekatan yang umum digunakan untuk model regresi logistik ordinal adalah dengan melihat nilai deviance residuals dan McFadden’s pseudo R-squared.

# McFadden's pseudo R-squared
# Perlu menghitung deviance dari model null dan model penuh

# Model null (tanpa prediktor)
model_null <- vglm(satisfaction ~ 1, family = cumulative(link = "logitlink"), data = df_ord)

# Model dengan prediktor
model_full <- vglm(satisfaction ~ food_quality, family = cumulative(link = "logitlink"), data = df_ord)

# Hitung McFadden's pseudo R-squared
deviance_null <- deviance(model_null)
deviance_full <- deviance(model_full)

pseudo_r2 <- 1 - (deviance_full / deviance_null)
pseudo_r2

## [1] 0.4109445

Interpretasi

Untuk mengevaluasi kecocokan model terhadap data, dilakukan penghitungan McFadden’s pseudo R-squared dengan membandingkan deviance dari model null (tanpa prediktor) dan model penuh (dengan prediktor food_quality). Hasil perhitungan menunjukkan bahwa nilai pseudo R² adalah sebesar 0.411.

Nilai ini menunjukkan bahwa model menjelaskan sekitar 41.1% variasi dalam data, yang termasuk dalam kategori cukup baik. Dalam regresi logistik, nilai McFadden’s pseudo R² antara 0.2 hingga 0.4 umumnya dianggap menunjukkan goodness of fit yang memadai. Oleh karena itu, model ordinal dengan prediktor food_quality memiliki kecocokan yang layak terhadap data.

---

12.9 Alternatif Model Ordinal

Selain model cumulative logit yang paling umum digunakan, terdapat beberapa alternatif model ordinal lain yang bisa dipertimbangkan, khususnya ketika asumsi proportional odds tidak terpenuhi. Dua model utama yang menjadi alternatif adalah:

Adjacent-Category Logit Model: Model ini membandingkan peluang antara satu kategori dengan kategori berikutnya secara langsung (misalnya kategori 1 vs 2, 2 vs 3, dan seterusnya). Model ini berguna ketika hubungan antar kategori bersifat lokal atau sekuensial.

Continuation-Ratio (Sequential) Logit Model: Model ini memodelkan peluang untuk “melanjutkan” ke kategori yang lebih tinggi, bersifat sekuensial, dan sangat cocok untuk data yang memiliki sifat hirarkis atau proses yang berurutan (misalnya tahapan pendidikan atau proses pengambilan keputusan).

Model-model alternatif ini memberikan fleksibilitas dalam analisis data ordinal, terutama ketika model cumulative logit tidak memenuhi asumsi paralelisme atau tidak memberikan hasil yang sesuai. Pemilihan model yang tepat sebaiknya didasarkan pada karakteristik data dan tujuan analisis.

---

13. Model Log-Linier

Analisis data kategorik merupakan aspek penting dalam statistika terapan, karena banyak fenomena nyata menghasilkan data dalam bentuk kategori, seperti jenis kelamin, status pekerjaan, tingkat pendidikan, preferensi konsumen, dan hasil diagnosis medis. Untuk menganalisis data kategori, terdapat tiga pendekatan statistik yang umum digunakan:

• 1. Tabel Kontingensi

  Tabel kontingensi digunakan sebagai langkah awal eksplorasi untuk melihat hubungan antara dua atau lebih variabel kategorik. Misalnya, dalam studi efek obat terhadap serangan jantung, tabel kontingensi dapat menyajikan jumlah pasien yang mengalami atau tidak mengalami serangan jantung berdasarkan jenis obat yang dikonsumsi. Dari tabel ini, dapat dihitung ukuran asosiasi seperti:

  • Odds ratio
  • Risk ratio
  • Statistik chi-square untuk uji independensi antar variabel

• 2. Model Log-Linier

  Model log-linier berguna untuk memahami struktur asosiasi dan interaksi antar variabel dalam tabel kontingensi. Beberapa karakteristik model log-linier:

  • Termasuk dalam keluarga Generalized Linear Model (GLM)
  • Menggunakan distribusi Poisson untuk memodelkan frekuensi sel
  • Tidak membedakan variabel dependen dan independen; semua variabel diperlakukan secara simetris

  Model ini sangat berguna ketika ingin melihat pengaruh gabungan antar variabel, termasuk interaksi dua arah atau tiga arah. Misalnya, dalam analisis hubungan antara jenis kelamin, status merokok, dan penyakit paru, model log-linier dapat digunakan untuk mengevaluasi apakah hanya interaksi dua arah yang cukup atau dibutuhkan interaksi tiga arah.

  Pemilihan model dapat diuji menggunakan Likelihood Ratio Test (LRT) untuk membandingkan model yang lebih sederhana dengan model yang lebih kompleks.

• 3. Model Regresi Logistik

Model regresi logistik digunakan ketika terdapat satu variabel kategorik yang dianggap sebagai variabel dependen, dan satu atau lebih variabel lain sebagai prediktor. Ciri-ciri utama:

Memodelkan logit dari probabilitas kejadian (log odds)

Cocok untuk studi observasional dan eksperimental

Dapat digunakan untuk outcome biner maupun kategori lebih dari dua (dengan perluasan model seperti regresi logistik multinomial dan ordinal)

Studi Kasus: Merokok, Jenis Kelamin, dan Penyakit Paru-paru

Kita memiliki data fiktif dari studi observasional untuk mengetahui apakah status merokok dan jenis kelamin berhubungan dengan kejadian penyakit paru-paru. Data disajikan dalam bentuk frekuensi sebagai berikut:

# Membuat data frekuensi
table_data <- array(
  data = c(40, 20, 30, 10, 15, 25, 10, 35),
  dim = c(2, 2, 2),
  dimnames = list(
    Merokok = c("Perokok", "Bukan Perokok"),
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Penyakit_Paru = c("Ya", "Tidak")
  )
)
table_data

## , , Penyakit_Paru = Ya
## 
##                Jenis_Kelamin
## Merokok         Laki-laki Perempuan
##   Perokok              40        30
##   Bukan Perokok        20        10
## 
## , , Penyakit_Paru = Tidak
## 
##                Jenis_Kelamin
## Merokok         Laki-laki Perempuan
##   Perokok              15        10
##   Bukan Perokok        25        35

Tabel kontingensi digunakan untuk menyajikan distribusi frekuensi gabungan dari beberapa variabel kategorik. Dalam kasus ini:

• Apakah perokok lebih sering mengalami penyakit paru-paru?

• Apakah jenis kelamin berpengaruh?

Tabel ini berguna untuk eksplorasi awal dan penghitungan statistik deskriptif seperti proporsi atau Chi-square test.

# Menggunakan fungsi margin.table untuk tabel 2D
table_2d <- margin.table(table_data, c(1, 3))  # Merokok vs Penyakit Paru
chisq.test(table_2d)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table_2d
## X-squared = 28.696, df = 1, p-value = 8.469e-08

Model log-linier cocok digunakan ketika semua variabel adalah kategorik, dan tidak ditetapkan variabel dependen.

Model log-linier mencoba menjelaskan struktur asosiasi antar variabel dalam bentuk logaritma dari frekuensi ekspektasi.

library(MASS)

# Ubah array ke data frame frekuensi
loglin_table <- as.table(table_data)
model_loglin <- loglm(~ Merokok * Jenis_Kelamin + Merokok * Penyakit_Paru + Jenis_Kelamin * Penyakit_Paru, data = loglin_table)
summary(model_loglin)

## Formula:
## ~Merokok * Jenis_Kelamin + Merokok * Penyakit_Paru + Jenis_Kelamin * 
##     Penyakit_Paru
## attr(,"variables")
## list(Merokok, Jenis_Kelamin, Penyakit_Paru)
## attr(,"factors")
##               Merokok Jenis_Kelamin Penyakit_Paru Merokok:Jenis_Kelamin
## Merokok             1             0             0                     1
## Jenis_Kelamin       0             1             0                     1
## Penyakit_Paru       0             0             1                     0
##               Merokok:Penyakit_Paru Jenis_Kelamin:Penyakit_Paru
## Merokok                           1                           0
## Jenis_Kelamin                     0                           1
## Penyakit_Paru                     1                           1
## attr(,"term.labels")
## [1] "Merokok"                     "Jenis_Kelamin"              
## [3] "Penyakit_Paru"               "Merokok:Jenis_Kelamin"      
## [5] "Merokok:Penyakit_Paru"       "Jenis_Kelamin:Penyakit_Paru"
## attr(,"order")
## [1] 1 1 1 2 2 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
## 
## Statistics:
##                       X^2 df   P(> X^2)
## Likelihood Ratio 3.029941  1 0.08174100
## Pearson          3.000435  1 0.08324215

Model Regresi Logistik Jika kita ingin memodelkan probabilitas terkena penyakit paru-paru berdasarkan status merokok dan jenis kelamin, maka regresi logistik adalah pendekatan yang tepat.

# Membuat data frame dari array
df <- as.data.frame(as.table(table_data))

# Ubah Penyakit_Paru menjadi biner
df$Penyakit_Paru_bin <- ifelse(df$Penyakit_Paru == "Ya", 1, 0)

# Model regresi logistik
glm_model <- glm(Penyakit_Paru_bin ~ Merokok + Jenis_Kelamin, data = df, family = binomial(), weights = Freq)
summary(glm_model)

## 
## Call:
## glm(formula = Penyakit_Paru_bin ~ Merokok + Jenis_Kelamin, family = binomial(), 
##     data = df, weights = Freq)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              1.2405     0.2806   4.422 9.80e-06 ***
## MerokokBukan Perokok    -1.7076     0.3248  -5.257 1.46e-07 ***
## Jenis_KelaminPerempuan  -0.4707     0.3252  -1.447    0.148    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 255.25  on 7  degrees of freedom
## Residual deviance: 221.97  on 5  degrees of freedom
## AIC: 227.97
## 
## Number of Fisher Scoring iterations: 4

Kesimpulan

Tabel kontingensi menunjukkan hubungan awal antar variabel.

Model log-linier cocok untuk mengeksplorasi struktur hubungan antar semua variabel kategorik.

Model regresi logistik digunakan jika fokus utama adalah prediksi terhadap outcome kategorik, seperti terkena penyakit atau tidak.

Masing-masing pendekatan memiliki kekuatan tergantung pada tujuan: eksplorasi, struktur, atau prediksi.

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(μ) ~ efek	GLM: logit(p) ~ prediktor
Cocok untuk	Eksplorasi awal	Tabel dengan > 2 variabel kategorik	Studi observasional / prediktif

---

13.1 Model Log-Linier: Konsep dan Aplikasi

Model log-linier merupakan bagian dari Generalized Linear Models (GLM) yang khusus dirancang untuk menganalisis hubungan antar variabel kategorik. Model ini memodelkan logaritma dari frekuensi harapan dalam tabel kontingensi sebagai kombinasi linier dari parameter efek utama dan interaksi.

Formulasi Matematis: \[ \log(\mu_{ijk}) = \lambda + \lambda_i^A + \lambda_j^B + \lambda_k^C + \lambda_{ij}^{AB} + \lambda_{ik}^{AC} + \lambda_{jk}^{BC} + \lambda_{ijk}^{ABC} \] di mana:

• \(\mu_{ijk}\) = frekuensi harapan untuk sel (i,j,k)

• \(\lambda\) = efek rata-rata umum

• \(\lambda^A, \lambda^B, \lambda^C\) = efek utama

• \(\lambda^{AB}, \lambda^{AC}, \lambda^{BC}\) = interaksi dua arah

• \(\lambda^{ABC}\) = interaksi tiga arah

---

13.2 Studi Kasus: Analisis Faktor Risiko Kanker Kulit

13.2.1 Persiapan Data

Kita menggunakan dataset simulasi dengan 4 variabel kategorik:

- Kanker Kulit (Ya/Tidak)

- Paparan UV (Rendah/Tinggi)

- Riwayat Keluarga (Ya/Tidak)

- Usia (<40, 40-60, >60 tahun)

# Pembuatan tabel kontingensi 4D
kanker_data <- array(
  c(15,30,10,25,20,40,5,15,10,20,30,50, # Data frekuensi
    8,22,12,18,15,35,8,12,5,15,25,45),  # Data tambahan untuk realisme
  dim = c(2,2,2,3),
  dimnames = list(
    Kanker = c("Ya","Tidak"),
    Paparan_UV = c("Rendah","Tinggi"),
    Riwayat_Keluarga = c("Ya","Tidak"),
    Usia = c("<40","40-60",">60")
  )
)

13.2.2 Model-Model Penting

1. Model Independen (Null Model)

Model ini mengasumsikan semua variabel independen/tidak berhubungan: \[ \log(\mu_{ijkl}) = \lambda + \lambda_i^K + \lambda_j^P + \lambda_k^R + \lambda_l^U \]

# Convert data to dataframe format
kanker_df <- as.data.frame(as.table(kanker_data))

# Fit models using glm()
model_independen_glm <- glm(Freq ~ Kanker + Paparan_UV + Riwayat_Keluarga + Usia,
                          data = kanker_df, family = poisson)

Interpretasi:

- Nilai devians tinggi menunjukkan ketidakcocokan model

- Derajat kebebasan besar karena banyak sel yang tidak terjelaskan

2. Model Hierarkis

Model ini mencakup semua efek utama dan interaksi dua arah: \[ \log(\mu_{ijkl}) = \lambda + \lambda_i^K + \lambda_j^P + \lambda_k^R + \lambda_l^U + \lambda_{ij}^{KP} + \lambda_{ik}^{KR} + \lambda_{il}^{KU} + \lambda_{jk}^{PR} + \lambda_{jl}^{PU} + \lambda_{kl}^{RU} \]

model_hierarkis_glm <- glm(Freq ~ (Kanker + Paparan_UV + Riwayat_Keluarga + Usia)^2,
                          data = kanker_df, family = poisson)

3. Model Saturated

Model lengkap yang memuat semua interaksi hingga tingkat tertinggi: \[ \log(\mu_{ijkl}) = \lambda + \lambda_i^K + ... + \lambda_{ijkl}^{KPRU} \]

model_saturated_glm <- glm(Freq ~ Kanker*Paparan_UV*Riwayat_Keluarga*Usia,
                          data = kanker_df, family = poisson)

13.2.3 Pemilihan Model

Tabel Perbandingan Model:

# Load required packages
library(MASS)
library(knitr)
# Create comparison table using loglm objects
comparison <- data.frame(
  Model = c("Independen", "Hierarkis", "Saturated"),
  Deviance = c(model_independen_glm$deviance, 
               model_hierarkis_glm$deviance,
               model_saturated_glm$deviance),
  df = c(model_independen_glm$df.residual,
         model_hierarkis_glm$df.residual,
         model_saturated_glm$df.residual)
)

# Add pseudo-AIC calculation for loglm objects (AIC = Deviance + 2*df)
comparison$AIC <- comparison$Deviance + 2*comparison$df

# Print formatted comparison table
kable(comparison, caption = "Perbandingan Goodness-of-Fit Model", digits = 2)

Perbandingan Goodness-of-Fit Model

Model	Deviance	df	AIC
Independen	107.56	18	143.56
Hierarkis	57.78	9	75.78
Saturated	0.00	0	0.00

Likelihood Ratio Test:

# Likelihood Ratio Tests
cat("\nLikelihood Ratio Tests:\n")

## 
## Likelihood Ratio Tests:

anova(model_independen_glm, model_hierarkis_glm, model_saturated_glm, test = "LR")

## Analysis of Deviance Table
## 
## Model 1: Freq ~ Kanker + Paparan_UV + Riwayat_Keluarga + Usia
## Model 2: Freq ~ (Kanker + Paparan_UV + Riwayat_Keluarga + Usia)^2
## Model 3: Freq ~ Kanker * Paparan_UV * Riwayat_Keluarga * Usia
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        18    107.564                          
## 2         9     57.781  9   49.782 1.184e-07 ***
## 3         0      0.000  9   57.781 3.580e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kriteria Pemilihan:

1. Parsimoni: Model sederhana yang cukup menjelaskan data

2. Signifikansi Statistik: LRT p-value < 0.05

3. AIC: Nilai lebih kecil menunjukkan model lebih baik

Interpretasi

Hasil uji likelihood ratio menunjukkan bahwa model hierarkis jauh lebih baik dari model independen (p < 0.001), dan model saturated paling baik karena deviasinya nol. Namun, karena model saturated terlalu kompleks, model hierarkis dipilih sebagai yang paling seimbang. Ini juga didukung oleh nilai AIC: model hierarkis memiliki AIC lebih rendah (75.78) dibandingkan model independen (143.56), yang berarti lebih cocok tanpa terlalu rumit.

13.2.4 Estimasi Parameter dan Odds Ratio

summary(model_hierarkis_glm)

## 
## Call:
## glm(formula = Freq ~ (Kanker + Paparan_UV + Riwayat_Keluarga + 
##     Usia)^2, family = poisson, data = kanker_df)
## 
## Coefficients:
##                                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                             2.79075    0.19365  14.411  < 2e-16 ***
## KankerTidak                             0.82574    0.20800   3.970 7.19e-05 ***
## Paparan_UVTinggi                       -0.58994    0.23467  -2.514 0.011941 *  
## Riwayat_KeluargaTidak                  -0.06758    0.21931  -0.308 0.757960    
## Usia40-60                              -0.20991    0.25060  -0.838 0.402227    
## Usia>60                                -0.48736    0.25753  -1.892 0.058426 .  
## KankerTidak:Paparan_UVTinggi           -0.16619    0.19911  -0.835 0.403903    
## KankerTidak:Riwayat_KeluargaTidak       0.04213    0.19571   0.215 0.829559    
## KankerTidak:Usia40-60                  -0.12594    0.24378  -0.517 0.605421    
## KankerTidak:Usia>60                    -0.05221    0.24375  -0.214 0.830404    
## Paparan_UVTinggi:Riwayat_KeluargaTidak  0.11236    0.19089   0.589 0.556110    
## Paparan_UVTinggi:Usia40-60              1.26532    0.23357   5.417 6.05e-08 ***
## Paparan_UVTinggi:Usia>60                0.88978    0.23089   3.854 0.000116 ***
## Riwayat_KeluargaTidak:Usia40-60        -0.63900    0.23312  -2.741 0.006124 ** 
## Riwayat_KeluargaTidak:Usia>60           0.22773    0.22833   0.997 0.318588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 166.587  on 23  degrees of freedom
## Residual deviance:  57.781  on  9  degrees of freedom
## AIC: 200.14
## 
## Number of Fisher Scoring iterations: 4

library(broom)

# Ambil semua estimasi parameter
estimasi <- tidy(model_hierarkis_glm) %>%
  mutate(
    OR = exp(estimate),
    CI_lower = exp(estimate - 1.96 * std.error),
    CI_upper = exp(estimate + 1.96 * std.error)
  )

# Tampilkan hasil
print(estimasi[, c("term", "estimate", "OR", "CI_lower", "CI_upper")])

## # A tibble: 15 × 5
##    term                                   estimate     OR CI_lower CI_upper
##    <chr>                                     <dbl>  <dbl>    <dbl>    <dbl>
##  1 (Intercept)                              2.79   16.3     11.1     23.8  
##  2 KankerTidak                              0.826   2.28     1.52     3.43 
##  3 Paparan_UVTinggi                        -0.590   0.554    0.350    0.878
##  4 Riwayat_KeluargaTidak                   -0.0676  0.935    0.608    1.44 
##  5 Usia40-60                               -0.210   0.811    0.496    1.32 
##  6 Usia>60                                 -0.487   0.614    0.371    1.02 
##  7 KankerTidak:Paparan_UVTinggi            -0.166   0.847    0.573    1.25 
##  8 KankerTidak:Riwayat_KeluargaTidak        0.0421  1.04     0.711    1.53 
##  9 KankerTidak:Usia40-60                   -0.126   0.882    0.547    1.42 
## 10 KankerTidak:Usia>60                     -0.0522  0.949    0.589    1.53 
## 11 Paparan_UVTinggi:Riwayat_KeluargaTidak   0.112   1.12     0.770    1.63 
## 12 Paparan_UVTinggi:Usia40-60               1.27    3.54     2.24     5.60 
## 13 Paparan_UVTinggi:Usia>60                 0.890   2.43     1.55     3.83 
## 14 Riwayat_KeluargaTidak:Usia40-60         -0.639   0.528    0.334    0.834
## 15 Riwayat_KeluargaTidak:Usia>60            0.228   1.26     0.803    1.96

---

14 Referensi

• Agresti, A. (2013). Categorical Data Analysis (3rd ed.). Wiley.
  
• Agresti, A. (2019). An Introduction to Categorical Data Analysis. Wiley.
  
• Bland, M. (2000). An Introduction to Medical Statistics (3rd ed.). Oxford University Press.
  
• Cox, D. R., & Snell, E. J. (1989). Analysis of Binary Data (2nd ed.). Chapman & Hall.
  
• Greenland, S. (1997). Methods for Epidemiologic Analysis of Multiple Exposures. Wiley.
  
• Greenland, S., Lash, T. L., & Rothman, K. J. (2008). Concepts of interaction. In Rothman, K., Greenland, S., & Lash, T. (Eds.), Modern Epidemiology (3rd ed., pp. 71–83). Lippincott Williams & Wilkins.
  
• Hosmer, D. W., & Lemeshow, S. (2000). Applied Logistic Regression (2nd ed.). Wiley.
  
• Hosmer, D. W., Lemeshow, S., & Sturdivant, R. X. (2013). Applied Logistic Regression (3rd ed.). Wiley.
  
• Kleinbaum, D. G., & Klein, M. (2010). Logistic Regression: A Self-Learning Text (3rd ed.). Springer.
  
• Levy, P. S., & Lemeshow, S. (2013). Sampling of Populations: Methods and Applications (4th ed.). Wiley.
  
• Rosner, B. (2015). Fundamentals of Biostatistics (8th ed.). Cengage Learning.