Ebook ADK

1 KATA PENGANTAR

Puji syukur penulis panjatkan kepada Tuhan yang Maha Esa karena atas berkat dan rahmatnya sehingga E-Book Analisis Kategori ini dapat disusun dan diselesaikan dengan baik. Buku ini disusun dalam rangka menyelesaikan pembelajaran analisis kategori pada prodi Statistika FMIPA Unpad dan memberikan pemahaman yang lebih mendalam teori dan penerapan analisis kategori dalam kehidupan.Dalam buku ini penulis mencoba untuk menjelaskan konsep-konsep dasar dan aplikasi dari analisis kategorik pada data ril.

Dalam penyusunan e-book ini, penulis mencoba untuk menuliskan teori-teori dan penjelasan dalam bahasa yang mudah dipahami mengingat analisis kategorik ini akan semakin banyak digunakan dalam kehidupan sehari-hari. Penulis berupaya menyajikan materi secara sistematis dan mencoba untuk memberikan contoh pada masing-masing teori.

Ucapan terima kasih penulis panjatkan kepada semua pihak yang sudah mendukung mulai dari Pak I Gede Nyoman Mindra Jaya selaku dosen pengampu mata kuliah analisis data kategorik dan teman-teman seperjuangan di statistika.

Akhir kata, penulis menyadari bahwa e-book ini masih jauh dari kata sempurna. Oleh karena itu kritik dan saran sangat terbuka untuk pembelajaran penulis kedepannya. Semoga buku ini dapat menjadi referensi yang berguna bagi siapapun yang membaca.

Jatinangor, Juni 2025 Penulis

2 PENDAHULUAN

Dalam statistik terdapat beberapa jenis data, salah satunya adalah data kategorik. Data kategorik adalah data yang bersifat non-numerik dan terdiri dari kategori atau kelompok.

Analisis data kategorik memilki peran penting pada bidang-bidang tersebut tetapi melakukan analisis dengan data yang tidak kontinu memerlukan pendekatan khusus seperti tabel kontingensi, uji chi-square, analisis regresi logistik, dan metode estimasi berbasis model probabilistik. Pendekatan-pendekatan ini memungkinkan untuk menganalisis data yang bersifat kategorik sehingga dapat digunakan pada banyak bidang.

2.1 Tujuan Analisis Data Kategorik

Masing-masing analisis data memilki tujuan tersendiri, analisis data kategorik memilki beberapa tujuan yaitu :

2.1.1 Menyajikan data secara visual

Untuk membantu memahami data, diperlukan pemvisualisan data dengan berbagai cara seperti diagram batang, pie chart, tabel kontingensi, dan masih banyak lagi. Hal ini juga dapat mengindetifikasi pola atau tren dari data.

2.1.2 Menganalisis hubungan antar variabel kategorik

Analisis data kategorik pada beberapa bidang memerlukan untuk melihat apakah ada hubungan dalam dua ketegori, contohnya apakah status pernikahan memiliki hubungan dengan prefensi produk, apakah ada hubungan antara tingkat pendidikan dan pekerjaan.

2.1.3 Membandingkan kelompok

Dalam beberapa bidang dibutuhkan perbandingan antar kelompok, contohnya tingkat kepuasan pelanggan berdasarkan wilayah, dan lainnya. Hal ini juga dapat membantu

2.2 Definisi dan Ruang Lingkup Analisis Data Kategori

Analisis data kategori adalah teknik statistik yang digunakan untuk menganalisis variabel yang bersifat kategorik. Ada beberapa jenis data yaitu :

2.2.1 Nominal vs Ordinal

Nominal adalah data yang tidak memiliki urutan. Contohnya jenis kelamin, warna, dsb

Ordinal adalah data yang memilki urutan. Contohnya tingkat pendidikan, status sosial dsb

2.2.2 Data Biner vs Multikategori

Biner adalah data yang hanya memiliki dua kategori. Contohnya ya/tidak, menang/kalah, dsb

Multikategori adalah data yang memiliki lebih dari dua kategori.

2.3 Perbedaan Data Kategorik dan Kuantitatif

Data kategorik tidak dapat diukur dalam skala numerik yang kontinu sedangkan data kuantitatif bisa

2.4 Manfaat Analisis Data Kategori dalam Berbagai Bidang

Data kategorik memiliki peran penting dalam berbagai bidang terutama di bidang kesehatan dan sosial. Contoh dalam bidang kesehatan data kategoriknya adalah pasien sehat (ya/tidak), tingkat kecelakaan seseorang (tinggi/sedang/rendah) atau pada bidang sosial, status sosial (tinggi/rendah/sedang), status menikah (kawin/tidak), bahkan sampai pekerjaan seseorang. Tidak hanya pada kedua bidang tersebut, masih banyak bidang lain yang juga menggunakan data kategorik.

2.4.1 Ilmu Sosial dan Psikologi

Contoh dalam bidang ini adalah

• Survei kepuasan pelanggan dengan skala Likert

• Analisis hubungan antara faktor sosial ekonomi dan tingkat pendidikan seseorang

• Penelitian efek terapi psikolog terhadap pasien

2.4.2 Kesehatan dan Kedokteran

Contoh dalam bidang kesehatan dan kedokteran

• Mengkategorikan status kesehatan pasien

• Efektivitas terapi

• Hubungan antara penderita sebuah penyakit dengan penyakit lainnya

2.4.3 Pemasaran dan Bisnis

Contoh dalam bidang kesehatan dan kedokteran

• Menentukan preferensi pelanggan terhadap produk tertentu

• Analisis segmentasi pelanggan berdasarkan usia, golongan, jenis kelamin

• Melihat kepuasan pelanggan

2.4.4 Pendidikan

Beberapa contoh dalam bidang pendidikan adalah

• Survei kepuasan pelajar terhadap metode pengajaran sekolah

• Analisis hubungan latar belakang ekonomi dan prestasi akademik

• Studi efektivitas program pembelajaran berbasis teknologi

2.4.5 Kebijakan Publik dan Pemerintahan

Contoh untuk bidang ini akan berkaitan dengan kebutuhan masyarakat :

• Analisis tingkat kepuasan masyarakat terhadap layanan publik

• Studi tentang tingkat partisipasi pemilu

• Evaluasi efektivitas program bantuan berdasarkan kategori penerima manfaat

2.4.6 Keamanan dan Kriminalitas

Contoh dalam bidang ini adalah

• Analisis kategori kejahatan yang paling sering terjadi di suatu wilayah

• STudi tentang faktor demografis yang berkorelasi dengan tingkat kriminalitas

• Evaluasi efektivitas kebijakan penegakan hukum terhadap kategori pelanggaran

3 METODE DALAM ANALISIS DATA KATEGORI

Dalam analisis data kategori terdapat banyak metode tergantung pada tujuan penelitian yang diinginkan. Metode-metode umum yang sering digunakan adalah :

3.1 Tabel Kontingensi dan Uji Chi-Square

• Digunakan untuk menguji hubungan antara dua variabel kategori

• Contohnya untuk menguji apakah ada hubungan latar belakang ekonomi dengan prestasi akademik

3.2 Regresi Logistik

• Digunakan untuk memprediksi probabilitas suatu kejadian berdasarkan variabel kategori

• Contohnya memprediksi kepuasan pelanggan berdasarkan tingkat pelayanannya

3.3 Analisis Correspondence

• Digunakan untuk memahami hubungan antara berbagai kategori dalam satu dataset

• Contohnya analisis tentang preferensi makanan berdasarkan kelompok usia

3.4 Decision Tree dan Random Forest

• Metode machine learning yang sering digunakan untuk klasifikasi berbasis data kategori

Dengan pendekatan yang tepat, analisis data kategori dapat membantu untuk memvisualisasikan data yang bersifat kategori. Selain itu analisis yang tepat memungkinkan pengambilan keputusan yang tepat dan akurat.

4 DISTRIBUSI PROBABILITAS DALAM DATA KATEGORI

Variabel acak kategori adalah variabel yang hanya memiliki beberapa kategori diskrit

4.1 Distribusi Bernoulli

Distribusi Bernoulli digunakan untuk percobaan biner, yaitu percobaan yang memiliki dua kemungkinan hasil: - Sukses (1) dengan probabilitas p - Gagal (0) dengan probabilitas 1 - p

Fungsi probabilitasnya:

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

Keterangan Notasi:

• X: Variabel acak biner (0 atau 1)

• p: Probabilitas sukses (X = 1)

Contoh Variabel Acak Bernoulli – Hasil dari lemparan koin (Kepala = 1, Ekor = 0) – Keberhasilan atau kegagalan dalam suatu percobaan klinis

4.2 Distribusi Binomial

Distribusi Binomial adalah generalisasi dari distribusi Bernoulli yang digunakan untuk n kali percobaan independen yang masing-masing memiliki dua kemungkinan hasil (sukses atau gagal). Jika setiap percobaan memiliki probabilitas sukses p, maka distribusi Binomial memiliki fungsi probabilitas:

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

Keterangan Notasi:

• X: Jumlah keberhasilan dalam n percobaan

• n: Jumlah percobaan

• k: Jumlah keberhasilan yang diamati

• p: Probabilitas keberhasilan dalam satu percobaan

• \(\binom{n}{k}\): Kombinasi “n pilih k”, dihitung sebagai \(\frac{n!}{k!(n-k)!}\)

Contoh Variabel Acak Binomial – Jumlah keberhasilan dalam 10 kali lemparan koin – Jumlah pasien yang sembuh setelah diberikan obat tertentu dalam suatu studi klinis

4.3 Distribusi Multinomial

Distribusi Multinomial adalah generalisasi lebih lanjut dari distribusi Binomial, digunakan ketika setiap percobaan memiliki lebih dari dua kemungkinan hasil. Jika suatu eksperimen dilakukan n kali, dan setiap percobaan dapat menghasilkan salah satu dari k kategori dengan probabilitas \(p_1,p_2,...,p_k\) maka distribusi probabilitasnya:

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

Keterangan Notasi

• \(X_{i}\): Frekuensi kemunculan kategori ke-i
• \(n_{}\): Jumlah total percobaan
• \(x_{i}\) Jumlah kejadian kategori ke-i
• \(p_{i}\): Probabilitas kategori ke-i

Contoh Variabel Acak Multinomial
- Pemilihan kandidat dalam pemilu (beberapa kandidat, satu suara per pemilih)

4.4 Distribusi Poisson

Distribusi Poisson digunakan untuk menghitung jumlah kejadian dalam interval waktu atau ruang tertentu dengan rata-rata kejadian \(\lambda\) per unit waktu/ruang. Fungsi probabilitasnya:

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

Keterangan Notasi

• X: Jumlah kejadian dalam interval tertentu

• \(\lambda\): Rata-rata kejadian dalam interval tersebut

• k: Jumlah kejadian yang diamati

Contoh Variabel Acak Poisson yaitu Jumlah panggilan telepon masuk ke pusat layanan dalam satu jam, Jumlah kecelakaan lalu lintas di satu jalan dalam sehari

5 DESAIN SAMPLING DALAM ANALISIS DATA KATEGORI

Dalam analisis data kategori, desain sampling memiliki peran yang krusial dalam menentukan validitas dan reliabilitas hasil penelitian. Pemilihan desain sampling yang tepat bergantung pada tujuan penelitian dan jenis data yang dikumpulkan.

Secara umum, desain sampling dalam analisis data kategori dapat diklasifikasikan ke dalam dua pendekatan utama, yaitu prospective sampling dan retrospective sampling. Masing-masing pendekatan ini memiliki karakteristik dan metode sampling yang berbeda, yang sering digunakan dalam penelitian eksperimental maupun observasional seperti eksperimen, studi kohort, dan studi kasus-kontrol.

5.1 Prospective Sampling

Prospective sampling adalah metode pengambilan sampel dimana subjek penelitian diidentifikasi dan diikuti dalam periode waktu tertentu untuk mengamati perkembangan variabel yang diteliti. Pendekatan ini memungkinkan peneliti untuk mengontrol variabel sebelum pengukuran hasil dilakukan, sehingga sering digunakan dalam studi kausal dan eksperimental. Beberapa jenis desain sampling dalam metode ini meliputi:

5.1.1 Eksperimen

Dalam studi eksperimental, subjek secara acak dialokasikan ke dalam kelompok perlakuan dan kontrol. 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 sampel diambil secara acak dari setiap strata.
• Cluster Sampling: Populasi dibagi menjadi kelompok-kelompok (cluster), kemudian beberapa cluster dipilih secara acak untuk dianalisis.

5.1.2 Studi Kohort

Studi kohort adalah penelitian observasional dimana kelompok individu dengan karakteristik tertentu diikuti dari waktu ke waktu untuk mengamati kejadian yang dipelajari. Jenis sampling yang umum dalam studi kohort meliputi:

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

5.2 Retrospective Sampling

Retrospective sampling adalah metode dimana data dikumpulkan dari peristiwa yang telah terjadi sebelumnya. Teknik ini sering digunakan dalam studi observasional yang mencari hubungan antara faktor risiko dan hasil tertentu.

5.2.1 Studi Kasus-Kontrol

Dalam studi kasus-kontrol, sekelompok individu dengan kondisi tertentu (kasus) dibandingkan dengan kelompok tanpa kondisi tersebut (kontrol). Teknik sampling yang sering digunakan 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 kemunculan kasus.

5.2.2 Studi Kohort Retrospektif

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

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

6 TABEL KONTINGENSI

Tabel kontingensi adalah tabel yang digunakan untuk menunjukkan hubungan antara dua (atau lebih) variabel kategorik. Bentuk dari tabel kontingensi adalah berbentuk matriks yang berisikan frekuensi dari hubungan antara dua variabel yang ada.

6.1 Tabel Kontingensi 2 × 2

Bentuk paling sederhana dari tabel kontingensi adalah tabel kontingensi 2x2 dimana tabel ini melihat apakah terdapat hubungan/asosiasi pada dua variabel.

Tabel kontingensi \(2 \times 2\) memiliki struktur sebagai berikut:

\[ \begin{array}{c|cc|c} & \text{Kategori 1 ( + )} & \text{Kategori 2 ( – )} & \text{Total} \\ \hline \text{Grup 1} & n_{11} & n_{12} & n_{1\cdot} \\ \text{Grup 2} & n_{21} & n_{22} & n_{2\cdot} \\ \hline \text{Total} & n_{\cdot 1} & n_{\cdot 2} & n \\ \end{array} \]

• \(n_{11}\) = Jumlah kasus dalam kategori positif dari Grup 1
• \(n_{12}\) = Jumlah kasus dalam kategori negatif dari Grup 1
• \(n_{21}\) = Jumlah kasus dalam kategori positif dari Grup 2
• \(n_{22}\) = Jumlah kasus dalam kategori negatif dari Grup 2
• \(n_{1\cdot}\) = Total observasi dalam Grup 1 (\(n_{1\cdot} = n_{11} + n_{12}\))
• \(n_{2\cdot}\) = Total observasi dalam Grup 2 (\(n_{2\cdot} = n_{21} + n_{22}\))
• \(n_{\cdot 1}\) = Total observasi dalam Kategori 1 (\(n_{\cdot 1} = n_{11} + n_{21}\))
• \(n_{\cdot 2}\) = Total observasi dalam Kategori 2 (\(n_{\cdot 2} = n_{12} + n_{22}\))
• \(n\) = Total keseluruhan sampel (\(n = n_{11} + n_{12} + n_{21} + n_{22}\))

Contoh tabel kontingensi 2x2 adalah sebagai berikut :

# Data Tabel Kontingensi 2x2
data <- matrix(c(25, 3, 28, 25, 4, 29, 50, 57 ), 
                nrow = 3, byrow = TRUE)

## Warning in matrix(c(25, 3, 28, 25, 4, 29, 50, 57), nrow = 3, byrow = TRUE):
## data length [8] is not a sub-multiple or multiple of the number of rows [3]

colnames(data) <- c("Tinggi/Sedang","Rendah","Total") 
rownames(data) <- c("Baik", "Cukup/Kurang","Tota")

kable(data, align = "c", 
      caption = "Tabel Kontingensi 2x2 Hasil Belajar Siswa terhadap Pengamalam Ibadah Sehari-Hari ")

Tabel Kontingensi 2x2 Hasil Belajar Siswa terhadap Pengamalam Ibadah Sehari-Hari

	Tinggi/Sedang	Rendah	Total
Baik	25	3	28
Cukup/Kurang	25	4	29
Tota	50	57	25

6.2 Peluang Tabel Kontingensi 2x2

6.2.0.1 Peluang Bersama

Peluang bersama adalah peluang bahwa kedua variabel terjadi secara bersamaan dalam suatu sel tabel kontingensi

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

6.2.0.2 Peluang Marginal

Peluang marginal adalah probabilitas kejadian suatu variabel tanpa mempertimbangkan variabel lainnya. Peluang ini adalah peluang kejadian hanya satu kategorinya:

• Peluang marginal baris:

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

• Peluang marginal kolom:

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

6.2.0.3 Peluang Bersyarat

Peluang bersyarat adalah probabilitas suatu kejadian terjadi dengan syarat kejadian lain telah terjadi:

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

6.2.0.4 Contoh

Dengan contoh tabel kontingensi sebelumnya, akan dihitung distribusi peluang dari tabel tersebut \[ \begin{array}{|l|c|c|c|} \hline \textbf{Tinggi/Sedang} & \textbf{Rendah} & \textbf{Total} \\ \hline \textbf{Baik} & 25 & 3 & 28 \\ \hline \textbf{Cukup/Kurang} & 25 & 4 & 29 \\ \hline \textbf{Total} & 50 & 7 & 57 \\ \hline \end{array} \]

#Tabel Kontingensi 2x2

\[ \begin{array}{|l|c|c|c|} \hline \textbf{Tinggi/Sedang} & \textbf{Rendah} & \textbf{Total} \\ \hline \textbf{Baik} & 25 & 3 & 28 \\ \hline \textbf{Cukup/Kurang} & 25 & 4 & 29 \\ \hline \textbf{Total} & 50 & 7 & 57 \\ \hline \end{array} \]

#Peluang bersama
n <- sum(data)
joint <- data/n

#Peluang marginal
marginalr <- rowSums(data)/n
marginalc <- colSums(data)/n

#Peluang bersyarat
conditional <- data/rowSums(data)

#Hasil <- 
list (Peluang_Bersama <- joint,
      Peluang_marginal_baris <- marginalr,
      Peluang_marginal_kolom <- marginalc,
      Peluang_bersyarat <- conditional)

## [[1]]
##              Tinggi/Sedang     Rendah     Total
## Baik              0.101626 0.01219512 0.1138211
## Cukup/Kurang      0.101626 0.01626016 0.1178862
## Tota              0.203252 0.23170732 0.1016260
## 
## [[2]]
##         Baik Cukup/Kurang         Tota 
##    0.2276423    0.2357724    0.5365854 
## 
## [[3]]
## Tinggi/Sedang        Rendah         Total 
##     0.4065041     0.2601626     0.3333333 
## 
## [[4]]
##              Tinggi/Sedang     Rendah     Total
## Baik             0.4464286 0.05357143 0.5000000
## Cukup/Kurang     0.4310345 0.06896552 0.5000000
## Tota             0.3787879 0.43181818 0.1893939

Interpretasi :

6.3 Ukuran Asosiasi Tabel Kontingensi 2x2

Asosiasi tabel kontingensi 2x2 digunakan untuk menentukan hubungan antara dua variabel kategori. Dengan frekuensi yang ada pada tabel kontingensi memungkinkan kita untuk menilai hubungan statistik yang terjadi. Ukuran asosiasi yang dapat dihitung adalah :

• Risk Difference (RD)

• Relative Risk (RR)

• Odds Ratio (OR)

• Uji Chi-Square & Fishers’s exact Test

6.3.1 Risk Difference (RD)

Risk Difference (RD) atau perbedaan resiko adalah ukuran yang mengambarkan perbedaan antara probabilitas kejadian suatu hasil yag sama dalam dua kelompok kategori yang berbeda. RD dihitung dengan selisih antara risiko kejadian dalam kelompok 1 dan kelompok 2.

\[ RD = \left( \frac{n_{11}}{n_{1.}} \right) - \left( \frac{n_{21}}{n_{2.}} \right) \]

Keterangan :

• Jika RD > 0, maka risiko kejadian pada kelompok 1 lebih tinggi dibandingkan kelompok 2

• Jika RD < 0, maka risiko kejadian pada kelompok 1 lebih rendah dibandingkan kelompok 2

• Jika RD = 0, maka tidak ada perbedaan risiko antara dua kelompok

6.3.2 Relative Risk (RR)

Relative risk (RR) atau risiko relatif adalah ukuran yang membandingkan risiko kejadian dua kelompok yang berbeda.

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

Keterangan :

• Jika RR > 1, maka kejadian kelompok 1 lebih sering terjadi dibandingkan kelompok 2

• Jika RR < 1, maka kejadian kelompok 1 lebih jarang terjadi dibandingkan kelompok 2

• Jika RR = 1, maka kejadian kedua kelompok memiliki risiko yang sama

6.3.3 Odds Ratio (OR)

Odds ratio (OR) atau rasio odds adalah ukuran yang membandingkan peluang (odds) terjadinya suatu kejadian antara dua kelompok

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

Keterangan :

• Jika OR > 1, maka peluang kejadian kelompok 1 lebih besar dibandingkan kelompok 2

• Jika OR < 1, maka peluang kejadian kelompok 1 lebih kecil dibandingkan kelompok 2

• Jika OR = 1, maka tidak ada perbedaan peluang kejadian kelompok 1 dan kelompok 2

6.3.4 Perbedaan RD,RR, dan OR

\[ \begin{array}{|l|p{4.5cm}|p{4.5cm}|p{5cm}|} \hline \textbf{Ukuran} & \textbf{Definisi} & \textbf{Desain Sampling yang Cocok} & \textbf{Interpretasi} \\ \hline \textbf{Risk Difference (RD)} & Selisih probabilitas kejadian antara dua kelompok & Studi kohort atau eksperimen acak & Menunjukkan tambahan atau pengurangan risiko absolut \\ \hline \textbf{Relative Risk (RR)} & Perbandingan risiko kejadian di dua kelompok & Studi kohort atau eksperimen klinis & Mengukur berapa kali lebih besar risiko di satu kelompok dibandingkan kelompok lainnya \\ \hline \textbf{Odds Ratio (OR)} & Perbandingan odds antara dua kelompok & Studi kasus-kontrol atau studi observasional & Mengukur hubungan antara variabel eksposur dan kejadian dalam desain studi kasus-kontrol \\ \hline \end{array} \]

Dari ukuran asosiasi, RD membandingkan selisih, RR membandingkan risiko, OR membandingkan odds.

6.3.5 Contoh

Dari tabel kontingensi yang ada akan dihitung ukuran asosiasi nya

\[ \begin{array}{|l|c|c|c|} \hline & \textbf{Tinggi/Sedang} & \textbf{Rendah} & \textbf{Total} \\ \hline \textbf{Baik} & 25 & 3 & 28 \\ \hline \textbf{Cukup/Kurang} & 25 & 4 & 29 \\ \hline \textbf{Total} & 50 & 7 & 57 \\ \hline \end{array} \]

# Data input
n11 <- 25  # Baik, Tinggi/Sedang
n12 <- 3   # Baik, Rendah
n21 <- 25  # Cukup/Kurang, Tinggi/Sedang
n22 <- 4   # Cukup/Kurang, Rendah

n1. <- n11 + n12  # Total Baik
n2. <- n21 + n22  # Total Cukup/Kurang

# Risk Difference (RD)
RD <- (n11 / n1.) - (n21 / n2.)

# Relative Risk (RR)
RR <- (n11 / n1.) / (n21 / n2.)

# Odds Ratio (OR)
OR <- (n11 * n22) / (n12 * n21)

# Tampilkan hasil
list(
  "Risk Difference (RD)" = RD,
  "Relative Risk (RR)" = RR,
  "Odds Ratio (OR)" = OR
)

## $`Risk Difference (RD)`
## [1] 0.03078818
## 
## $`Relative Risk (RR)`
## [1] 1.035714
## 
## $`Odds Ratio (OR)`
## [1] 1.333333

7 INFERENSI TABEL KONTINGENSI DUA ARAH

Inferensi dalam statistik adalah proses pengambilan keputusan untuk populasi berdasarkan data sampel. Dalam tabel kontingensi dua arah, inferensi digunakan untuk menganalisis hubungan anatra dua variabel kategorik yang disusun dalam tabel kontingensi. Inferensi dilakukan melalui estimasi dan pengujian hipotesis.

7.1 Estimasi

Dalam inferensi estimasi menjadi penting yaitu memperkirakan paramete populasi berdasarkan data sampel. Estimasi dibagi menjadi dua, estimasi titik dan estimasi interval.

7.1.1 Estimasi Titik

Estimasi titik digunakan untuk menentukan satu nilai spesifik sebegai perkiraan terbaik untuk populasi berdasarkan data sampel. Estimasi titik dihitung sebagai berikut :

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

Keterangan :

• \(\hat{p}\) adalah estimasi titik proporsi
  
• \(x\) adalah jumlah individu dalam kategori tertentu
  
• \(n\) adalah total jumlah individu dalam sampel

7.1.2 Estimasi Interval

Selain estimasi titik, ada juga estimasi interval dimana hasil dari perkiraan merupakan rentangg nilai dengan tingkat kepercayaan tertentu. Estimasi interval lebih sering digunakan dibandingkan estimasi titik karena dianggap lebih akurat dalam untuk memprediksi dalam suatu rentang dibandingkan hanya satu titik. Estimasi interval dapat dihitung sebagai berikut :

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

Keterangan :

• \(Z_{\alpha/2}\) adalah nilai dari distribusi normal standar untuk tingkat kepercayaan tertentu,
  
• \(\hat{p}\) adalah estimasi titik proporsi,
  
• \(n\) adalah ukuran sampel.

7.2 Hipotesis

7.2.1 Uji Proporsi

Uji proporsi digunakan untuk membandingkan apakah ada perbedaaan signfikan antara dua kejadian pada dua kelompok yang berbeda.

Hipotesis untuk uji z dua proporsi yaitu :

• Hipotesis Nol (\(H_0\)): Tidak ada perbedaan proporsi antara dua kelompok, yaitu \(p_1 = p_2\)

• Hipotesis Alternatif (\(H_1\)): Terdapat perbedaan proporsi antara dua kelompok, yaitu \(p_1 \ne p_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)}} \]

dimana :

Estimasi proporsi dalam masing-masing kelompok diberikan oleh:

\[ \hat{p}_1 = \frac{n_{11}}{n_{1.}}, \quad \hat{p}_2 = \frac{n_{21}}{n_{2.}} \]

Estimasi proporsi gabungan (pooling proportion):

\[ \hat{p} = \frac{n_{11} + n_{21}}{n_{1.} + n_{2.}} \]

Daerah kritisnya mengikuti tabel Z pada distribusi normal baku N(0,1). Jika |Z| lebih besar dibandingkan nilai kritisnya maka \(H_0\) ditolak, sedangkan jika hal lain maka \(H_0\) diterima.

7.2.2 Uji Asosiasi

Seperti yang sudah dibahas pada bab sebelumnya, ukuran asosiasi mengukur hubungan antara dua variabel kategorik. Terdapat 3 uji dari uji asosiasi, risk difference (RD), relative risk (RR), dan odds ratio (OR).

Untuk hipotesis uji asosiasi untuk ketiga pengujian adalah sama tetapi akan berbeda dalam menghitung statistik uji nya :

• Hipotesis Nol (\(H_0\)): Tidak ada asosiasi antara dua variabel

• Hipotesis Alternatif (\(H_1\)): Terdapat asosiasi antara dua variabel

7.2.2.1 Risk Difference (RD)

Risk Difference (RD) mengukur selisih dari probabilitas dua kejadiaan antara dua kelompok.

\[ RD = \left( \frac{n_{11}}{n_{1.}} \right) - \left( \frac{n_{21}}{n_{2.}} \right) \]

Dimana standard error untuk RD :

\[ SE_{RD} = \sqrt{ \frac{\hat{p}_1 (1 - \hat{p}_1)}{n_1} + \frac{\hat{p}_2 (1 - \hat{p}_2)}{n_2} } \] Untuk statistik uji RD dapat dihitung dengan :

\[ Z = \frac{RD}{SE(RD)} =\frac{ \frac{n_{11}}{n_{1.}} - \frac{n_{21}}{n_{2.}} }{ \sqrt{ \frac{n_{11}(n_{1.} - n_{11})}{n_{1.}^3} + \frac{n_{21}(n_{2.} - n_{21})}{n_{2.}^3} } } \]

7.2.2.2 Relative Risk (RR)

Relative risk membandingkan kemungkinan kejadian antara dua kelompok.

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

Standard error untuk log(RR) :

\[ SE_{\ln(RR)} = \sqrt{ \frac{1}{n_{11}} - \frac{1}{n_{1.}} + \frac{1}{n_{21}} - \frac{1}{n_{2.}} } \] Untuk statistik uji RR dapat dihitung dengan :

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

7.2.2.3 Odds Ratio

Odds Ratio membandingkan peluang kejadian antara dua kelompok

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

Standard Error untuk log(OR)

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

Untuk statistik uji OR dapat dihitung dengan :

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

7.2.3 Contoh Uji Proporsi dan Uji Asosiasi

Tabel kontingensi 2x2 yang menulis hubunagn antara hasil belajar dan amalan ibadah sehari-hari

\[ \begin{array}{|l|c|c|c|} \hline & \textbf{Tinggi/Sedang} & \textbf{Rendah} & \textbf{Total} \\ \hline \textbf{Baik} & 25 & 3 & 28 \\ \hline \textbf{Cukup/Kurang} & 25 & 4 & 29 \\ \hline \textbf{Total} & 50 & 7 & 57 \\ \hline \end{array} \]

Uji proporsi :

# Membuat matriks dari data tabel
data <- matrix(c(25, 3, 25, 4), nrow = 2, byrow = TRUE)
dimnames(data) <- list(
  Kualitas = c("Baik", "Cukup/Kurang"),
  Kecerdasan = c("TinggiSedang", "Rendah")
)

# Tampilkan data
print(data)

##               Kecerdasan
## Kualitas       TinggiSedang Rendah
##   Baik                   25      3
##   Cukup/Kurang           25      4

#Uji Proporsi Dua Sampel
# Lakukan uji proporsi dua sampel
prop_test <- prop.test(
  x = c(data[1, 1], data[2, 1]),  # jumlah sukses di masing-masing grup
  n = c(sum(data[1, ]), sum(data[2, ])),  # total dalam masing-masing grup
  correct = FALSE  # tanpa continuity correction
)

## Warning in prop.test(x = c(data[1, 1], data[2, 1]), n = c(sum(data[1, ]), :
## Chi-squared approximation may be incorrect

# Tampilkan hasil uji
print(prop_test)

## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  c(data[1, 1], data[2, 1]) out of c(sum(data[1, ]), sum(data[2, ]))
## X-squared = 0.12535, df = 1, p-value = 0.7233
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.1391392  0.2007155
## sample estimates:
##    prop 1    prop 2 
## 0.8928571 0.8620690

Didapatkan nilai p-value untuk tabel kontingensi yang adalah sebesar 0,7233. Untuk p-value > alpha (0,05) maka \[H_0\] diterima yang artinya tidak terdapat perbedaan proporsi antara hasil belajar baik dan kurang dengan amalan ibadah tinggi.

Uji Asosiasi :

n11 <- data["Baik", "TinggiSedang"]
n12 <- data["Baik", "Rendah"]
n21 <- data["Cukup/Kurang", "TinggiSedang"]
n22 <- data["Cukup/Kurang", "Rendah"]

n1 <- n11 + n12
n2 <- n21 + n22

p1 <- n11 / n1
p2 <- n21 / n2

#Risk Difference
RD <- p1 - p2
SE_RD <- sqrt((p1 * (1 - p1)) / n1 + (p2 * (1 - p2)) / n2)
Z_RD <- RD / SE_RD
p_RD <- 2 * (1 - pnorm(abs(Z_RD)))

cat("Risk Difference (RD):", round(RD, 4), "\n")

## Risk Difference (RD): 0.0308

cat("SE(RD):", round(SE_RD, 4), "\n")

## SE(RD): 0.0867

cat("Z-statistik RD:", round(Z_RD, 4), "\n")

## Z-statistik RD: 0.3551

cat("p-value RD:", round(p_RD, 4), "\n")

## p-value RD: 0.7225

#Relative Risk 
RR <- p1 / p2
SE_logRR <- sqrt(1/n11 - 1/n1 + 1/n21 - 1/n2)
Z_RR <- log(RR) / SE_logRR
p_RR <- 2 * (1 - pnorm(abs(Z_RR)))

cat("Relative Risk (RR):", round(RR, 4), "\n")

## Relative Risk (RR): 1.0357

cat("SE(log(RR)):", round(SE_logRR, 4), "\n")

## SE(log(RR)): 0.099

cat("Z-statistik RR:", round(Z_RR, 4), "\n")

## Z-statistik RR: 0.3544

cat("p-value RR:", round(p_RR, 4), "\n")

## p-value RR: 0.723

#Odds Ratio 
OR <- (n11 * n22) / (n12 * n21)
SE_logOR <- sqrt(1/n11 + 1/n12 + 1/n21 + 1/n22)
Z_OR <- log(OR) / SE_logOR
p_OR <- 2 * (1 - pnorm(abs(Z_OR)))

cat("Odds Ratio (OR):", round(OR, 4), "\n")

## Odds Ratio (OR): 1.3333

cat("SE(log(OR)):", round(SE_logOR, 4), "\n")

## SE(log(OR)): 0.8145

cat("Z-statistik OR:", round(Z_OR, 4), "\n")

## Z-statistik OR: 0.3532

cat("p-value OR:", round(p_OR, 4), "\n")

## p-value OR: 0.7239

Dari hasil pengujian nilai p-value nya didapatkan bahwa dnegan ketiga uji tidak ada asosisasi antara dua kelompok tersebut.

7.2.4 Uji Independensi

Uji independensi digunakan untuk menguji apakah data yang dimiliki memiliki hubungan antara dua variabel. Pengujian independensi dilakukan dengan uji Chi-Square dan jika tabel kontingensi tidak berukuran 2x2 atau berbentuk persegi dapat menggunakan uji Chi-Square parsial

7.2.4.1 Uji Chi-Square

Uji Chi-Square digunakan untuk menguji apakah ada hubungan (independensi) antara variabel kategorik. Untuk menghitung Chi-Square digunakan rumus sebagai berikut :

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

di mana:

• \(O\) adalah nilai observasi dalam tabel kontingensi.
• \(E\) adalah nilai yang diharapkan, dihitung sebagai:

\[ E_{ij} = \frac{(R_i \times C_j)}{N} \]

dengan:

• \(R_i\) = total baris ke-\(i\)
• \(C_j\) = total kolom ke-\(j\)
• \(N\) = total sampel

7.2.4.2 Partisi Chi-Square

Partisi Chi-Square digunakan untuk mengindetifikasi kategori mana dalam tabel kontingensi yang bertanggung jawab dalam hubungan yang signifikan. Jika dalam keseluruhan uji chi-square sudah signifikan, maka partisi Chi-square dapat menguraikan hubungan dalam subkelompok yang lebih kecil.

7.2.4.3 Uji Likelihood Ratio (G²)

Uji Likelihood Ratio (G²) adalah alternatif dari uji chi-square yang digunakan untuk menguji hipotesis independensi dalam tabel kontingensi \(I \times J\). Statistik uji ini diberikan oleh:

\[ G^2 = 2 \sum_i \sum_j n_{ij} \ln \left( \frac{n_{ij}}{\hat{\mu}_{ij}} \right) \]

Di mana: - \(n_{ij}\): frekuensi observasi - \(\hat{\mu}_{ij}\): nilai ekspektasi = \(\frac{R_i \cdot C_j}{N}\)

Statistik \(G^2\) mengikuti distribusi chi-square dengan derajat bebas \((I-1)(J-1)\).
Tolak \(H_0\) jika \(G^2 \geq \chi^2_{(1-\alpha), (I-1)(J-1)}\)

7.2.4.4 Uji Exact Fisher

Uji Fisher’s Exact digunakan untuk menguji hubungan dua variabel kategorikal dalam tabel kontingensi, dimana asumsi Chi-square tidak berlaku karena sampel yang sangat kecil

Keunggulan :

• Cocok untuk ukuran sampel kecil

• Tidak memerlukan asumsi normalitas atau chi-square

• Memberikan hasil yang lebih akurat pada data yang frekuensinya kecil

Keterbatasan :

• Perhitungan bisa menjadi berat jika ukuran besar

• Jika ukuran tabel kontingensi besar akan menyulitkan sehingga hanya cocok untuk tabel kontingensi kecil (2x2 atau 3x3)

7.2.4.5 Contoh Uji Independensi

Data yang adalah tabel kontingensi antara amalan ibadah dengan hasil pembelajaran

\[ \begin{array}{|l|c|c|c|} \hline & \textbf{Tinggi/Sedang} & \textbf{Rendah} & \textbf{Total} \\ \hline \textbf{Baik} & 25 & 3 & 28 \\ \hline \textbf{Cukup/Kurang} & 25 & 4 & 29 \\ \hline \textbf{Total} & 50 & 7 & 57 \\ \hline \end{array} \]

Uji Chi-Square

• Hipotesis Nol (\(H_0\)): Tidak ada hubungan antara variabel amalan ibadah dan hasil pembelajaran

• Hipotesis Alternatif (\(H_1\)): Ada hubungan antara variabel amalan ibadah dan hasil pembelajaran

{r} # Membuat matriks dari data tabel data <- matrix(c(25, 3, 25, 4), nrow = 2, byrow = TRUE) dimnames(data) <- list( Kualitas = c("Baik", "Cukup/Kurang"), Kecerdasan = c("TinggiSedang", "Rendah") ) # Melakukan uji chi-square uji_chisq <- chisq.test(data, correct = FALSE) uji_chisq}

dengan nilai p-value didapatkan nilai p-value sebesar 0,7233 yang artinya tidak ada hubungan antar variabel amalan ibadah dan hasil belajar.

Uji Likelihood Ratio :

{r} library(DescTools) # Uji G-Test GTest(data, correct = "none")}

Dari hasil uji likelihood ratio didapatkan bahwa nilai p-value lebih besar dibandingkan alpha, maka tidak ada hubungan antara dua variabel

Uji Exact Fisher

fisher.test(data)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  data
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.2014724 10.0012500
## sample estimates:
## odds ratio 
##   1.326645

Dari hasil exact fisher, didapatkan nilai p-value 1 yang artinya tidak ada hubungan antaar variabel hasil belajar dan amalan ibadah.

Dari ketiga uji independensi (Chi-square, Likelihood Ratio, Exact Fisher) didapatkan bahwa kedua variabel sebenarnya tidak memiliki hubungan. Hal ini memungkinkan karena data awalnya berupa data 3x3 yang disatukan untuk salah satu kategorinya.

7.3 Analisi Residual dalam Tabel Kontingensi

Residual dalam tabel kotingensi digunakan untuk mengindentifikasi sel mana yang menyumbang paling banyak terhadap hubungan anatara variabel kategori. Residual mengukur selisih frekuensi yang diamati dengan frekuensi yang diharapkan berdasarkan model independensi. Jika residualnya besar berarti frekuensi observasi pada sel tersebut berbeda dengan yang diharapkan, sedangkan kalau kecil berarti frekuensi palin mendekati nilai ekspektasi.

Ada beberapa jenis residual salah satunya adalah Pearson Residual

7.3.1 Pearson Residual

Rumus Pearson Residual:

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

Dimana:

• \(O_{ij}\): nilai observasi pada sel ke-\(i,j\)

• \(E_{ij}\): nilai ekspektasi pada sel ke-\(i,j\)

Standarized Residual (Adjusted Residual)

\[ r_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}(1 - p_{i+})(1 - p_{+j})}} \]

Dimana:

• \(p_{i+}\) dan \(p_{+j}\): probabilitas marginal dari baris dan kolom

7.3.1.1 Contoh

Dengan tabel kontingensi :

\[ \begin{array}{|l|c|c|c|} \hline & \textbf{Tinggi/Sedang} & \textbf{Rendah} & \textbf{Total} \\ \hline \textbf{Baik} & 25 & 3 & 28 \\ \hline \textbf{Cukup/Kurang} & 25 & 4 & 29 \\ \hline \textbf{Total} & 50 & 7 & 57 \\ \hline \end{array} \]

akan dihitung Pearson Residualnya dan Adjusted Pearson residualnya

obs <- matrix(c(25, 3, 25, 4), nrow = 2, byrow = TRUE) 
colnames(obs) <- c("Tinggi/Sedang", "Rendah")
rownames(obs) <- c("Baik", "Cukup/Kurang")

#Hitung nilai harapan (ekspektasi)

total <- sum(obs) 
row_totals <- rowSums(obs) 
col_totals <- colSums(obs)

exp <- outer(row_totals, col_totals) / total

#Pearson Residual

pearson_res <- (obs - exp) / sqrt(exp) 
round(pearson_res, 3)

##              Tinggi/Sedang Rendah
## Baik                 0.088 -0.237
## Cukup/Kurang        -0.087  0.232

Dari hasil perhitunggan pearson residualnya, didapatkan nilai yang cukup kecil dan mendekati 0, yang artinya perbedaan nilai observasi dan nilai ekspektasi tidak signifikan pada masing-masing sel.

7.3.2 Pendeteksian Outlier dengan Residual

Ketika nilai residual sangat besar, secara positif atau negatif maka hal ini disebut outlier. Outlier menunjukkan bahwa nilai nya jauh lebih tinggi atau lebih rendah dibandingkan nilai ekspektasi.

Sel terindikasi sebagai outlier ketika nilai Pearson Residual lebih besar dari 2. Jika nilai Adjusted Residual lebih besar dari 3 maka sel dianggap outlier signifikan.

8 TABEL KONTINGENSI TIGA ARAH

Tabel kontingensi tiga arah adalah perpanjangan dari tabel kontingensi dua arah yang digunakan untuk menganalisis hubungan antara tiga variabel kategori secara simultan. Misal dalam tabel kontingensi dua arah ada hubungan dua variabel (X dan Y) dipengaruhi dengan variabel ketiga (Z). Variabel ketiga (Z) disebut variabel kontrol atau kovariat.

Kegunaan tabel marginal adalah untuk melihat pola asosiasi secara agregat, tetapi sering kali mengabaikan efek kovariat yang dapat memberikan pemahaman lebih mendalam. Oleh karena itu, analisis yang mempertimbangkan tabel parsial biasanya lebih diutamakan untuk mendapatkan kesimpulan yang lebih akurat dalam penelitian ilmiah.

Tabel kotingensi tiga arah dapat ditulis menjadi dua tabel yaitu tabel parsial dan tabel marginal

8.1 Tabel Parsial dan Marginal

Tabel parsial adalah tabel yang mengelompokkan X dan Y berdasarkan setiap level Z, sedangkan tabel marginal adalah tabel yang mengabaikan Z dengan jumlah data semua level Z

Struktur tabel kontingensi tiga arah :

\[ \begin{array}{|c|c|c|c|c|c|} \hline Z & X, Y=1 & Y=2 & \cdots & Y=J & Y=\cdot \\ \hline \text{Z=1, X=1} & n_{111} & n_{121} & \cdots & n_{1J1} & n_{1+1} \\ \text{Z=1, X=2} & n_{211} & n_{221} & \cdots & n_{2J1} & n_{2+1} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \text{Z=1, X=I} & n_{I11} & n_{I21} & \cdots & n_{IJ1} & n_{I+1} \\ \hline \text{Z=k, X=1} & n_{11k} & n_{12k} & \cdots & n_{1Jk} & n_{1+k} \\ \text{Z=k, X=2} & n_{21k} & n_{22k} & \cdots & n_{2Jk} & n_{2+k} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \text{Z=k, X=I} & n_{I1k} & n_{I2k} & \cdots & n_{IJk} & n_{I+k} \\ \hline X = \cdot & n_{+1+} & n_{+2+} & \cdots & n_{+J+} & n_{+++} \\ \hline \end{array} \]

Keterangan:

• \(n_{ijk}\): Frekuensi observasi untuk \(X = i\), \(Y = j\), dan \(Z = k\)

• \(n_{i++}, n_{+j+}, n_{++k}\): Frekuensi marginal

Tabel Frekuensi marginal menampilkan total observasi untuk setiap variabel dengan mengabaikan variabel lainnya yang ada dalam tabel kontingensi tiga arah. Tabel ini nantinya bersisa frekuensi dari tabel kontingensi tiga arah dengan variabel tersisa.

Contoh data :

Tabel kontingensi berikut merupakan hasil penelitian yang diperoleh dari data Pusat Penelitian Pendapat Nasional yang melakukan survey social umum 1975 untuk melihat perbandingan dari jenis kelamin dan tingkat pendidikan responden dalam menentukan sikap terhadap wanita.

Variabel nya :

X : Masa pendidikan responden ( \(\le\) 8 , 9–12, dan \(\ge\) 13)

Y : Jawaban responden yaitu setuju atau tidak setuju dengan pernyataan

Z : Jenis kelamin responden \[ \begin{array}{|c|c|c|c|} \hline Z, X & Y = \text{Setuju} & Y = \text{Tidak Setuju} & Y = \cdot \\ \hline \text{Pria, } \leq 8 & 72 & 47 & 119 \\ \text{Pria, } 9{-}12 & 110 & 106 & 216 \\ \text{Pria, } \geq 13 & 44 & 187 & 231 \\ \hline \text{Wanita, } \leq 8 & 85 & 23 & 108 \\ \text{Wanita, } 9{-}12 & 173 & 28 & 201 \\ \text{Wanita, } \geq 13 & 210 & 87 & 297 \\ \hline X = \cdot & 694 & 488 & 1182 \\ \hline \end{array} \]

Dari tabel frekuensi tersebut dapat dijadikan tabel parsial sebagai berikut :

Tabel frekuensi parsial untuk Z (Jenis Kelamin) = Pria

\[ \begin{array}{|l|c|c|c|} \hline \text{X (Pendidikan)} & \text{Y = Setuju} & \text{Y = Tidak Setuju} & \text{Total} \\ \hline \leq 8 & 72 & 47 & 119 \\ 9{-}12 & 110 & 106 & 216 \\ \geq 13 & 44 & 187 & 231 \\ \hline \text{Total} & 226 & 340 & 566 \\ \hline \end{array}\]

Tabel frekuensi parsial untuk Z (Jenis Kelamin) = Wanita

\[ \begin{array}{|l|c|c|c|} \hline \text{X (Pendidikan)} & \text{Y = Setuju} & \text{Y = Tidak Setuju} & \text{Total} \\ \hline \leq 8 & 85 & 23 & 108 \\ 9{-}12 & 173 & 28 & 201 \\ \geq 13 & 210 & 87 & 297 \\ \hline \text{Total} & 468 & 138 & 606 \\ \hline \end{array} \]

Dari tabel frekuensi parsial dapat dihidung tabel frekuensi marginal sebagai berikut :

Tabel Frekuensi Marginal untuk X (Pendidikan) dan Y (Jawaban)

\[ \begin{array}{|l|c|c|c|} \hline \text{X (Pendidikan)} & \text{Y = Setuju} & \text{Y = Tidak Setuju} & \text{Total} \\ \hline \leq 8 & 72 + 85 = 157 & 47 + 23 = 70 & 227 \\ 9{-}12 & 110 + 173 = 283 & 106 + 28 = 134 & 417 \\ \geq 13 & 44 + 210 = 254 & 187 + 87 = 274 & 528 \\ \hline \text{Total} & 694 & 478 & 1172 \\ \hline \end{array} \]

Tabel Frekuensi Marginal untuk Z (Jenis Kelamin) dan Y (Jawaban)

\[ \begin{array}{|l|c|c|c|} \hline \text{Z (Jenis Kelamin)} & \text{Y = Setuju} & \text{Y = Tidak Setuju} & \text{Total} \\ \hline \text{Pria} & 226 & 340 & 566 \\ \text{Wanita} & 468 & 138 & 606 \\ \hline \text{Total} & 694 & 478 & 1172 \\ \hline \end{array} \]

8.2 Distribusi Peluang

8.2.1 Peluang Bersama (Joint Probability)

Peluang bersama didefinisikan sebagai :

\[ P(Z, X, Y) = \frac{f(Z, X, Y)}{N} \]

8.2.2 Peluang Marjinal (Marginal Probability)

• Peluang marginal untuk \(Z\):

\[ P(Z) = \sum_{x} \sum_{y} P(Z, X=x, Y=y) \]

• Peluang marginal untuk \(X\):

\[ P(X) = \sum_{z} \sum_{y} P(Z=z, X, Y=y) \]

• Peluang marginal untuk \(Y\):

\[ P(Y) = \sum_{z} \sum_{x} P(Z=z, X=x, Y) \]

---

8.2.3 Peluang Bersyarat (Conditional Probability)

Peluang bersyarat menunjukkan peluang suatu kejadian dengan asumsi kejadian lain telah terjadi.

• Peluang \(X\) dan \(Y\) diberikan \(Z\):

\[ P(X, Y \mid Z) = \frac{P(Z, X, Y)}{P(Z)} \]

• Peluang \(Y\) diberikan \(X\) dan \(Z\):

\[ P(Y \mid X, Z) = \frac{P(Z, X, Y)}{P(Z, X)} \]

• Peluang \(X\) diberikan \(Y\):

\[ P(X \mid Y) = \frac{P(X, Y)}{P(Y)} \]

8.2.4 Contoh

Gunakan data berikut (Tabel 5):

• Responden total: \(N = 1172\)
• Variabel:
  • \(Z\): Jenis Kelamin (Pria, Wanita)
  • \(X\): Pendidikan ( \(\le\) 8 , 9–12, dan \(\ge\) 13)
  • \(Y\): Jawaban (Setuju, Tidak Setuju)

8.2.4.1 Peluang Bersama

Contoh: Peluang responden adalah Pria, pendidikan 9–12 tahun, dan menjawab Setuju.

\[ P(Z = \text{Pria}, X = 9{-}12, Y = \text{Setuju}) = \frac{110}{1172} \approx 0.0939 \]

8.2.4.2 Peluang Marginal

Contoh: Peluang marginal responden menjawab Setuju:

\[ P(Y = \text{Setuju}) = \frac{694}{1172} \approx 0.592 \]

Contoh lain: Peluang marginal responden adalah Wanita:

\[ P(Z = \text{Wanita}) = \frac{606}{1172} \approx 0.517 \]

---

8.2.4.3 Peluang Bersyarat

Contoh: Peluang responden menjawab Setuju, diberikan bahwa responden adalah Pria.

\[ P(Y = \text{Setuju} \mid Z = \text{Pria}) = \frac{226}{566} \approx 0.399 \]

Contoh lain: Peluang responden adalah Pendidikan >=, diberikan bahwa menjawab Tidak Setuju.

\[ P(X = \geq 13 \mid Y = \text{Tidak Setuju}) = \frac{274}{478} \approx 0.573 \]

Implementasi dengan R :

data <- data.frame(
  Jenis_Kelamin = rep(c("Pria", "Wanita"), each = 3),
  Pendidikan = rep(c("<=8", "9-12", "<= 13"), times = 2),
  Setuju = c(72, 110, 44, 85, 173, 210),
  Tidak_Setuju = c(47, 106, 187, 33, 28, 187)
)

# Hitung total per baris dan total keseluruhan
data$Total <- data$Setuju + data$Tidak_Setuju
total_responden <- sum(data$Total)

# Joint probability: P(Z = Pria, X = 9-12, Y = Setuju)
joint_prob <- data$Setuju[data$Jenis_Kelamin == "Pria" & data$Pendidikan == "9-12"] / total_responden
joint_prob

## [1] 0.08580343

# Peluang marginal untuk Setuju
p_setuju <- sum(data$Setuju) / total_responden
# Peluang marginal untuk Tidak Setuju
p_tidak_setuju <- sum(data$Tidak_Setuju) / total_responden

c(Peluang_Setuju = p_setuju, Peluang_Tidak_Setuju = p_tidak_setuju)

##       Peluang_Setuju Peluang_Tidak_Setuju 
##            0.5413417            0.4586583

# Total responden pria
total_pria <- sum(data$Total[data$Jenis_Kelamin == "Pria"])
# Jumlah pria yang setuju
setuju_pria <- sum(data$Setuju[data$Jenis_Kelamin == "Pria"])
# Conditional probability
p_setuju_given_pria <- setuju_pria / total_pria

p_setuju_given_pria

## [1] 0.3992933

8.3 Ukuran Asosiasi

Seperti ukuran asosiasi untuk tabel kontingensi dua arah, ukuran asosiasi untuk tabel kontingensi tiga arah terdiri dari tiga yaitu Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR).

Rumus untuk masing-masing ukuran asosiasi adalah sebagai berikut :

1. Risk Difference antara dua kelompok didefinisikan sebagai:

\[ RD = P(Y|X_1, Z) - P(Y|X_2, Z) \]

2. Risiko Relatif

Risiko relatif membandingkan probabilitas kejadian antara dua kelompok:

\[ RR = \frac{P(Y|X_1, Z)}{P(Y|X_2, Z)} \]

3. Odds Ratio

Odds ratio adalah perbandingan antara odds dua kelompok:

\[ OR = \frac{\frac{P(Y|X_1, Z)}{1 - P(Y|X_1, Z)}}{\frac{P(Y|X_2, Z)}{1 - P(Y|X_2, Z)}} \]

8.3.1 Contoh

# Data input
data <- data.frame(
  Jenis_Kelamin = rep(c("Pria", "Wanita"), each = 3),
  Pendidikan = rep(c("<=8", "9-12", ">=13"), times = 2),
  Setuju = c(72, 110, 44, 153, 173, 28),
  Tidak_Setuju = c(47, 106, 187, 85, 83, 187)
)

data

##   Jenis_Kelamin Pendidikan Setuju Tidak_Setuju
## 1          Pria        <=8     72           47
## 2          Pria       9-12    110          106
## 3          Pria       >=13     44          187
## 4        Wanita        <=8    153           85
## 5        Wanita       9-12    173           83
## 6        Wanita       >=13     28          187

Ukuran Asosisasi Pendidikan terhadap setuju untuk jenis kelamin = pria

# Filter data pria
pria <- subset(data, Jenis_Kelamin == "Pria")

# Ambil data untuk X1 dan X2
x1 <- subset(pria, Pendidikan == "<=8")
x2 <- subset(pria, Pendidikan == ">=13")

# Hitung probabilitas
p_y_x1 <- x1$Setuju / (x1$Setuju + x1$Tidak_Setuju)
p_y_x2 <- x2$Setuju / (x2$Setuju + x2$Tidak_Setuju)

# Ukuran asosiasi
bp <- p_y_x1 - p_y_x2
rr <- p_y_x1 / p_y_x2
or <- (p_y_x1 / (1 - p_y_x1)) / (p_y_x2 / (1 - p_y_x2))

bp; rr; or

## [1] 0.4145658

## [1] 3.176471

## [1] 6.510638

Ukuran Asosisasi Pendidikan terhadap setuju untuk jenis kelamin = wanita

# Filter data wanita
wanita <- subset(data, Jenis_Kelamin == "Wanita")

# Ambil data untuk X1 dan X2
x1_w <- subset(wanita, Pendidikan == "<=8")
x2_w <- subset(wanita, Pendidikan == ">=13")

# Hitung probabilitas
p_y_x1_w <- x1_w$Setuju / (x1_w$Setuju + x1_w$Tidak_Setuju)
p_y_x2_w <- x2_w$Setuju / (x2_w$Setuju + x2_w$Tidak_Setuju)

# Ukuran asosiasi
bp_w <- p_y_x1_w - p_y_x2_w
rr_w <- p_y_x1_w / p_y_x2_w
or_w <- (p_y_x1_w / (1 - p_y_x1_w)) / (p_y_x2_w / (1 - p_y_x2_w))

# Tampilkan hasil
bp_w; rr_w; or_w

## [1] 0.5126246

## [1] 4.936224

## [1] 12.02143

Kesimpulan :

1. Risk Difference

   Pria : pria dengan pendidikan <= 8 tahun memiliki probabilitias 41,46% lebih tinggi untuk setuju dibandingkan pria dengan pendidikan >= 13 tahun

   Wanita : wanita dengan pendidikan <= 8 tahun memiliki probabilitias 51,26% lebih tinggi untuk setuju dibandingkan weanita dengan pendidikan >= 13 tahun

2. Risk Ratio

   Pria : pria dengan pendidikan <=8 tahun memiliki kemungkinan setuju yang 3.18 kali lebih besar dibanding pria dengan pendidikan >=13 tahun.

   Wanita : wanita dengan pendidikan <=8 tahun memiliki kemungkinan setuju yang 4.94 kali lebih besar dibanding wanita dengan pendidikan >=13 tahun.

3. Odds Ratio

   Pria : odds (peluang relatif) untuk setuju pada pria dengan pendidikan <=8 tahun adalah 6.5 kali lebih besar dibanding pria dengan pendidikan >=13 tahun.

   Wanita : odds untuk setuju pada wanita dengan pendidikan <=8 tahun adalah 12 kali lebih besar dibanding wanita dengan pendidikan >=13 tahun.

8.4 Conditional Independence

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

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

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

Atau dalam bentuk frekuensi:

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

Keterangan:

\(n_{ijk}\): Frekuensi pengamatan untuk kategori ke-i, j, dan k dari variabel X, Y, dan Z

\(n_{k++}\): Total frekuensi untuk kategori ke-k dari Z (semua X dan Y)

\(n_{i+k}\): Jumlah frekuensi untuk X = i, Z = k (semua Y)

\(n_{+jk}\): Jumlah frekuensi untuk Y = j, Z = k (semua X)

Conditional independence sangat penting dalam analisis statistik karena membantu mengidentifikasi hubungan yang bersifat tidak langsung antar variabel.

8.4.1 Pengujan Conditional Independence

Pengujian ini bermaksud untuk menguji apakah dua variabel menjadi independen ketika dikendalikan oleh variabel ketiga. Rumus secara teoritis :

Secara matematis:

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

Dalam bentuk frekuensi:

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

Untuk menguji conditional independence digunakan uji Chi-square dengan tabel kontingensi tiga dimensi.

8.4.2 Contoh

# Data asli
data <- data.frame(
  Jenis_Kelamin = rep(c("Pria", "Wanita"), each = 3),
  Pendidikan = rep(c("<=8", "9-12", ">=13"), times = 2),
  Setuju = c(72, 110, 44, 153, 173, 28),
  Tidak_Setuju = c(47, 106, 187, 85, 83, 187)
)

# Ubah ke format array 3D: [Pendidikan, Persetujuan, Jenis Kelamin]
array_data <- array(
  data = c(
    72, 47,   # Pria, <=8
    110, 106, # Pria, 9-12
    44, 187,  # Pria, >=13
    153, 85,  # Wanita, <=8
    173, 83,  # Wanita, 9-12
    28, 187   # Wanita, >=13
  ),
  dim = c(3, 2, 2),
  dimnames = list(
    Pendidikan = c("<=8", "9-12", ">=13"),
    Persetujuan = c("Setuju", "Tidak_Setuju"),
    Jenis_Kelamin = c("Pria", "Wanita")
  )
)

array_data

## , , Jenis_Kelamin = Pria
## 
##           Persetujuan
## Pendidikan Setuju Tidak_Setuju
##       <=8      72          106
##       9-12     47           44
##       >=13    110          187
## 
## , , Jenis_Kelamin = Wanita
## 
##           Persetujuan
## Pendidikan Setuju Tidak_Setuju
##       <=8     153           83
##       9-12     85           28
##       >=13    173          187

mantelhaen.test(array_data)

## 
##  Cochran-Mantel-Haenszel test
## 
## data:  array_data
## Cochran-Mantel-Haenszel M^2 = 33.568, df = 2, p-value = 5.138e-08

Kesimpulan :

Dari hasil pengujian dengan chi-square yang sudah digabungkan dengan CMH Test didapatkan p-value < 0,05 yang artinya terdapat ketergantungan antara pendidikan dan persetujuan meskipun sudah dikontrol jenis kelamin. Dimana ini artinya tidak ada conditional independence.

8.5 Inferensi Tabel Kontingensi Tiga Arah

Jika odds ratio relatif konstan, dapat menghitung odds ratio bersama dengan menggunakan estimasi Mantel Haenszel.

8.5.1 Independensi Bersyarat dalam Tabel Kontingensi Tiga Arah

Independensi bersyarat adalah konsep penting dalam analisis tabel kontingensi tiga arah. Ini merujuk pada kondisi di mana dua variabel, X dan Y, independen dalam setiap level variabel ketiga, Z.

Pengujian independensi bersyarat dilakukan dengan metode statistik seperti uji Cochran–Mantel–Haenszel (CMH).

Definisi Independensi Bersyarat

Dua variabel X dan Y dikatakan independen bersyarat terhadap variabel ketiga Z, jika rasio odds (OR) mereka dalam setiap strata Z sama dengan 1.

Secara matematis, ini dapat dituliskan sebagai:

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

Artinya, setelah mengendalikan pengaruh Z, tidak ada hubungan antara X dan Y dalam setiap strata.

Hal yang Perlu Diperhatikan

• Meskipun X dan Y independen dalam setiap strata Z, ini tidak berarti bahwa mereka juga independen dalam tabel marginal (tanpa mempertimbangkan Z).
  
• Fenomena ini dikenal sebagai Paradoks Simpson.

Paradoks Simpson sangat penting diperhatikan dalam analisis data epidemiologi, ilmu sosial, dan model prediktif karena dapat menyebabkan kesalahan dalam interpretasi hubungan antar variabel.

8.5.2 Pengujian Statistik untuk Independensi Bersyarat

Tujuan Uji CMH adalah untuk menguji hubungan antara dua variabel kategori dengan mempertimbangkan efek dari variabel perancu (confounding variable).

Kegunaan:

• Menguji hubungan yang dikendalikan oleh faktor perancu dalam tabel kontingensi berlapis (stratified contingency tables)
• Mengurangi bias dari variabel ketiga
• Menguji apakah asosiasi antara variabel kategori dengan outcome berbeda dari strata ke strata

Ide Dasar Uji CMH

Uji CMH bekerja dengan konsep tabel kontingensi berlapis (stratified 2 x 2 tables), yang masing-masing lapisan mewakili strata dari variabel ketiga (Z).

Uji CMH akan menguji apakah terdapat hubungan antara treatment (X) dan outcome (Y), namun juga mempertimbangkan efek dari variabel ketiga (Z), sehingga tidak terjadi confounding.

Hipotesis Uji CMH

• Hipotesis nol (\(H_0\)):
  \[ \theta_{XY|k} = 1 \quad \text{untuk setiap } k = 1, 2, ..., K \] (Tidak ada hubungan antara X dan Y di semua strata Z)

• Hipotesis alternatif ($H_1$$):
  \[ \theta_{XY|k} \ne 1 \quad \text{untuk paling sedikit satu } k \] (Ada hubungan antara X dan Y pada setidaknya satu strata Z)

Statistik Uji untuk uji CMH :

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

Keterangan :

• \(n_{11k}\): nilai frekuensi sel baris 1 kolom 1 pada tabel parsial ke-k
• \(\mu_{11k}\): nilai ekspektasi (nilai harapan) dari \(n_{11k}\), dihitung dengan rumus:

\[ \mu_{11k} = E(n_{11k}) = \frac{n_{1.k} \cdot n_{.1k}}{n_{..k}} \]

• \(\text{var}(n_{11k})\): varians dari \(n_{11k}\), diberikan oleh:

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

Keputusan uji :

• Statistik CMH mengikuti distribusi Chi-square dengan derajat kebebasan (df) = 1
• Tolak hipotesis nol \(H_0\) jika:

\[ CMH > \chi^2_{(1)} \quad \text{atau} \quad \text{p-value} < \alpha \]

8.5.3 Odds Ratio Bersama

Penaksir (Khusus Tabel 2 x 2 x K)

• Dalam tabel kontingensi 2 x 2 x K, terdapat K tabel parsial, sehingga ada sebanyak K odds ratio bersyarat.
• Jika nilai-nilai odds ratio bersyarat relatif sama (tidak berbeda secara ekstrem) dan memiliki arah yang sama, maka dapat ditentukan sebuah nilai tunggal odds ratio yang disebut odds ratio bersama.
• Odds ratio bersama ditaksir oleh statistik Mantel-Haenszel.

---

Rumus Odds Ratio Bersama

Odds ratio bersama ditaksir menggunakan rumus:

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

Dengan:

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

Standard Error untuk log Odds Ratio Bersama

Standard error untuk log odds ratio bersama dihitung dengan rumus:

\[ \hat{\sigma}^2[\log(\hat{\theta}_{MH})] = \frac{ \sum_k \left( \frac{n_{11k} + n_{22k}}{n_{11k} \cdot n_{22k}} + \frac{n_{12k} + n_{21k}}{n_{12k} \cdot n_{21k}} \right) }{ 4 \cdot \sum_k \frac{n_{11k} \cdot n_{22k}}{n_{..k}} \cdot \sum_k \frac{n_{12k} \cdot n_{21k}}{n_{..k}} } \]

Interval Kepercayaan log Odds Ratio Bersama

\[ \log(\hat{\theta}_{MH}) \pm Z_{\alpha/2} \cdot SE[\log(\hat{\theta}_{MH})] \]

Contoh Odds Ratio Bersama

Dari data delapan studi yang membahas mengenai hubungan merokok dan kanker paru-paru di 8 kota di China, kita akan menentukan odds ratio bersama, standard error, dan interval kepercayaan 95%.

8.5.3.1 Contoh

Dengan tabel kontingensi tiga arah yang digunakan pada bab ini, akan dicari nilai odds ratio bersama nya

# Input data
data <- data.frame(
  Jenis_Kelamin = rep(c("Pria", "Wanita"), each = 3),
  Pendidikan = rep(c("<=8", "9-12", ">=13"), times = 2),
  Setuju = c(72, 110, 44, 153, 173, 28),
  Tidak_Setuju = c(47, 106, 187, 85, 83, 187)
)

data

##   Jenis_Kelamin Pendidikan Setuju Tidak_Setuju
## 1          Pria        <=8     72           47
## 2          Pria       9-12    110          106
## 3          Pria       >=13     44          187
## 4        Wanita        <=8    153           85
## 5        Wanita       9-12    173           83
## 6        Wanita       >=13     28          187

# Subset untuk pria dan wanita
pria <- subset(data, Jenis_Kelamin == "Pria")
wanita <- subset(data, Jenis_Kelamin == "Wanita")

# Ambil pendidikan ekstrem untuk tabel 2x2
pria_strata <- data.frame(
  Pendidikan = c("<=8", ">=13"),
  Setuju = c(72, 44),
  Tidak_Setuju = c(47, 187)
)

wanita_strata <- data.frame(
  Pendidikan = c("<=8", ">=13"),
  Setuju = c(153, 28),
  Tidak_Setuju = c(85, 187)
)

# Jadikan bentuk tabel
tables <- list(
  t(as.matrix(pria_strata[, -1], rownames.force = NA)),
  t(as.matrix(wanita_strata[, -1], rownames.force = NA))
)

tables

## [[1]]
##              [,1] [,2]
## Setuju         72   44
## Tidak_Setuju   47  187
## 
## [[2]]
##              [,1] [,2]
## Setuju        153   28
## Tidak_Setuju   85  187

mantel_haenszel_or <- function(tables) {
  num <- 0
  denom <- 0
  for (table in tables) {
    a <- table["Setuju", 1]
    b <- table["Setuju", 2]
    c <- table["Tidak_Setuju", 1]
    d <- table["Tidak_Setuju", 2]
    n <- a + b + c + d
    num <- num + (a * d) / n
    denom <- denom + (b * c) / n
  }
  or_mh <- num / denom
  return(or_mh)
}

or_mh <- mantel_haenszel_or(tables)
or_mh

##   Setuju 
## 9.104422

Dari hasil perhitungan didapatkan nilai odds ratio bersama yaitu sebesar 9 dimana hal ini diartikan bahwa respondengn dengan pendidikan \(\le\) 8 tahun memiliki kemungkinan setuju terhadap pernyataan sekitar 9 kali lebih besar dibandingkan yang berpendidikan >=13 tahun

8.5.4 Uji Homogenitas Odds Ratio dengan Statistik Breslow-Day

Definisi Asosiasi Homogen

Asosiasi homogen terjadi jika odds ratio pada setiap tabel parsial bernilai sama:

\[ \theta_{XY(1)} = \theta_{XY(2)} = \cdots = \theta_{XY(K)} \]

Jika odds ratio konstan di semua strata variabel kontrol (\(Z\)), maka tidak ada interaksi antara variabel \(X\) dan \(Y\) dalam pengaruhnya terhadap \(Z\).

Jika odds ratio berbeda-beda antar level \(Z\), maka terdapat interaksi antara \(X\) dan \(Y\) dalam hubungannya dengan \(Z\).

Pengujian Homogenitas dengan Statistik Breslow-Day

Hipotesis :

\[ \begin{aligned} H_0 &: \theta_{XY(1)} = \theta_{XY(2)} = \cdots = \theta_{XY(K)} \\ H_1 &: \text{Setidaknya ada satu } \theta_{XY(k)} \text{ yang berbeda} \end{aligned} \]

Statistik Uji Breslow-Day

\[ X^2_{\text{BD}} = \sum_{j=1}^{K} \frac{(a_j - \tilde{a}_j)^2}{\text{Var}(\tilde{a}_j \mid \hat{\theta}_{MH})} \]

Keterangan :

• \(a_j\) : jumlah kasus terpapar yang diamati dalam strata \(j\)
• \(\tilde{a}_j\) : jumlah kasus terpapar yang diharapkan berdasarkan hipotesis nol
• \(\text{Var}(\tilde{a}_j \mid \hat{\theta}_{MH})\) : varians dari \(\tilde{a}_j\)

Statistik uji ini mengitu distribusi chi square dengan db = k-1. Jika BD lebih besar dari nilai kritis maka \(H_0\) ditolak.

8.5.4.1 Contoh

Dengan tabel kontingensi sebelumnya akan diuji homogenitas odds rationya dnegan Breslow-Day

# Input data
breslow_array <- array(
  c(
    # Strata 1: Pendidikan <=8
    72, 47,
    153, 85,
    # Strata 2: Pendidikan 9–12
    110, 106,
    173, 83,
    # Strata 3: Pendidikan >=13
    44, 187,
    28, 187
  ),
  dim = c(2, 2, 3),
  dimnames = list(
    Jenis_Kelamin = c("Pria", "Wanita"),
    Sikap = c("Setuju", "Tidak Setuju"),
    Pendidikan = c("<=8", "9-12", ">=13")
  )
)

# Tampilkan isi array
breslow_array

## , , Pendidikan = <=8
## 
##              Sikap
## Jenis_Kelamin Setuju Tidak Setuju
##        Pria       72          153
##        Wanita     47           85
## 
## , , Pendidikan = 9-12
## 
##              Sikap
## Jenis_Kelamin Setuju Tidak Setuju
##        Pria      110          173
##        Wanita    106           83
## 
## , , Pendidikan = >=13
## 
##              Sikap
## Jenis_Kelamin Setuju Tidak Setuju
##        Pria       44           28
##        Wanita    187          187

# Gunakan fungsi mantelhaen.test dari base R
bd_result <- mantelhaen.test(breslow_array, correct = FALSE)

bd_result

## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  breslow_array
## Mantel-Haenszel X-squared = 4.0372, df = 1, p-value = 0.04451
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.6063417 0.9947430
## sample estimates:
## common odds ratio 
##           0.77663

library (DescTools)
breslow_test <- BreslowDayTest (breslow_array)
print(breslow_test)

## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  breslow_array
## X-squared = 12.984, df = 2, p-value = 0.001516

Kesimpulan :

Dengan nilai common odds rasio sebesar 0,7766 menunjukkan bahwa pria memiliki kemungkinan setuju yang lebih rendah dibandingkan wanita setelah dikontrol oleh variabel lama pendidikan. Hasil dari Breslow Test yang dilakukan adalah nilai p-value sebesar 0,0015 dimana nilai ini lebih kecil dari alpha yang artinya odds ratio tidak homogen di seluruh strata

9 GENERALIZED LINEAR MODEL (GLM)

Generalized Linear Model adalah perluasan dari model regresi linear klasik. GLM membuat peneliti dapat memodelkan data yang variabelnya tidak berdistribusi normal dan/atau antar variabel prediktor dengan rata-rata dari respons tidak linear.

GLM terdiri dari 3 komponen yaitu :

1. Distribusi dari exponential family untuk variabel respons
2. Fungsi link yang menghubungkan ekspektasi dari variabel respons ke kombinasi linear dari prediktor
3. Fungsi linear prediktor : \(\eta = X\beta\)

9.1 Exponential Family

Distribusi termasuk dalam exponential family jika dapat ditulis dalam bentuk:

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

Contoh distribusi yang termasuk dalam exponential family:

• Distribusi Normal
• Distribusi Binomial
• Distribusi Poisson
• Distribusi Gamma

Contoh Pembuktian: Distribusi Binomial

Fungsi probabilitas distribusi binomial adalah:

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

Kita tuliskan ulang dalam bentuk exponential family:

\[ P(Y = y) = \exp\left\{ \log \binom{n}{y} + y \log \left( \frac{\pi}{1 - \pi} \right) + n \log(1 - \pi) \right\} \]

Dengan substitusi:

• \(\theta = \log \left( \frac{\pi}{1 - \pi} \right)\)
• \(b(\theta) = -n \log(1 - \pi)\)
• \(\phi = 1\)

Maka distribusi binomial termasuk dalam exponential family.

9.2 Model Regresi Logistik

Persamaan regresi logistik menyerupai regresi linear, dimana nilai input dikombinasikan secara linear dengan koefisien untuk menghasilkan prediksi. Hasil prediksi dari regresi logistik menjadi nilai biner, yaitu 0 atau 1 dengan menggunakan fungsi aktivasi sigmoid.

Regresi logistik menganalisis hubungan antara satu atau lebih variabel independen dan mengklasifkasikan data ke dalam kelas-kelas diskrit. Model ini banyak digunakan dalam pemodelan prediktif, di mana model memperkirakan probabilitas matematis apakah suatu entitas termasuk ke dalam kategori tertentu atau tidak.

Sebagai contoh, angka 0 dapat mewakili kelas negatif, dan angka 1 mewakili kelas positif. Regresi logistik biasanya digunakan untuk masalah klasifkasi biner, di mana variabel hasil hanya memiliki dua kemungkinan kategori (0 dan 1). Beberapa contoh penerapan klasifkasi biner di mana respons biner diharapkan atau tersirat antara lain:

• Menentukan probabilitas serangan jantung: Dengan bantuan model logistik, tenaga medis dapat menentukan hubungan antara variabel-variabel seperti berat badan, olahraga, dan sebagainya, untuk memprediksi apakah seseorang berisiko mengalami serangan jantung atau komplikasi medis lainnya.

• Kemungkinan diterima di universitas: Sistem aplikasi dapat memperkirakan probabilitas seorang siswa untuk diterima di universitas tertentu atau program studi tertentu dengan menganalisis hubungan antara variabel-variabel penentu seperti skor GRE, GMAT, atau TOEFL.

• Mengidentifkasi email spam: Kotak masuk email diflter untuk menentukan apakah suatu email merupakan komunikasi promosi atau spam dengan cara memahami variabel prediktor dan menerapkan algoritma regresi logistik untuk memeriksa keasliannya.

Keunggulan Utama Regresi Logistik

Lebih mudah diterapkan dalam metode pembelajaran mesin Modelpembelajaran mesin dapat dibangun secara efektif dengan menggunakan proses training (pelatihan) dan testing (pengujian). Proses pelatihan bertujuan untuk mengenali pola dalam data masukan (misalnya gambar) dan mengaitkannya dengan keluaran tertentu (label). Melatih model logistik dengan algoritma regresi tidak membutuhkan daya komputasi yang tinggi. Oleh karena itu, regresi logistik lebih mudah untuk diimplementasikan, ditafsirkan, dan dilatih dibandingkan metode machine learning lainnya.

1. Cocok untuk data yang dapat dipisahkan secara linear Dataset yang dapat dipisahkan secara linear mengacu pada grafk di mana dua kelas data dapat dipisahkan oleh garis lurus. Dalam regresi logistik, variabel respons (y) hanya memiliki dua nilai. Oleh karena itu, jika data bersifat dapat dipisahkan secara linear, maka klasifkasi ke dalam dua kelas berbeda dapat dilakukan secara efektif.

2. Memberikan wawasan yang berharga Regresi logistik dapat mengukur seberapa relevan atau penting suatu variabel independen/prediktor (melalui ukuran koefsien), serta menunjukkan arah hubungan atau asosiasi antara prediktor dan respons (apakah positif atau negatif).

3. Persamaan dan Asumsi dalam Regresi Logistik Regresi logistik menggunakan fungsi logistik yang disebut fungsi sigmoid untuk memetakan prediksi dan probabilitasnya. Fungsi sigmoid adalah kurva berbentuk huruf S (S-shaped curve) yang mengubah setiap nilai riil menjadi rentang antara 0 dan 1.

Selain itu, jika keluaran dari fungsi sigmoid (yaitu probabilitas yang diperkirakan) lebih besar dari ambang batas yang telah ditentukan dalam grafk, maka model akan memprediksi bahwa suatu observasi termasuk dalam kelas tertentu. Sebaliknya, jika nilai probabilitas tersebut lebih kecil dari ambang batas, model akanmemprediksi bahwa observasi tersebut tidak termasuk ke dalam kelas tersebut.

Sebagai contoh

• Jika hasil fungsi sigmoid lebih dari 0.5, maka output dianggap sebagai 1 (kelas positif).
• Jika hasilnya kurang dari 0.5, maka output diklasifikasikan sebagai 0 (kelas negatif).
• Jika grafik menuju ke arah negatif secara ekstrem, maka nilai prediksi \(y\) akan menjadi 0, dan sebaliknya.

Dengan kata lain, jika keluaran fungsi sigmoid adalah 0.65, maka itu berarti terdapat peluang sebesar 65% bahwa peristiwa tersebut akan terjadi — misalnya dalam kasus pelemparan koin.

Rumus fungsi sigmoid :

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

jika hasil dari fungsi sigmoid lebih besar dari ambang batas (misalnya 0.5), maka diklasifikasikan sebagai 1. Jika kurang dari 0.5, diklasifikasikan sebagai 0.

Simulasi dan Visualisasi Regresi Logistik

# Simulasi data untuk regresi logistik
set.seed(42)
n <- 100
x <- seq(-4, 4, length.out = n)
log_odds <- -0.5 + 1.5 * x
prob <- 1 / (1 + exp(-log_odds))
y <- rbinom(n, 1, prob)

# Buat data frame
data <- data.frame(x = x, y = y, prob = prob)

Plot kurva sigmoid

# Visualisasi menggunakan base R
plot(x, y, pch = 16, col = "gray60",
     xlab = "X", ylab = "Y / Probabilitas",
     main = "Simulasi Regresi Logistik dengan Kurva Sigmoid")

lines(x, prob, col = "blue", lwd = 2)
abline(h = 0.5, col = "red", lty = 2)

legend("topleft",
       legend = c("Data Biner (0/1)", "Kurva Logistik", "Ambang 0.5"),
       col = c("gray60", "blue", "red"),
       pch = c(16, NA, NA),
       lty = c(NA, 1, 2),
       lwd = c(NA, 2, 1),
       pt.cex = 1.5,
       bty = "n")

Kurva sigmoid dalam regresi logistik menunjukkan hubungan non-linear antara variabel prediktor dan probabilitas output. Pendekatan ini efektif untuk klasifikasi biner seperti deteksi penyakit, email spam, dan prediksi ya/tidak. Spesifikasi model :

Fungsi link (logit): Fungsi link function logit dapat dinyatakan dalam bentuk berikut:

\[ g(\mu) = \log\left(\frac{\mu}{1 - \mu}\right) \]

Model regresi logistik :

\[ \log\left(\frac{\mu}{1 - \mu}\right) = X\beta \]

Fungsi inverse link :

\[ \mu = \frac{\exp(X\beta)}{1 + \exp(X\beta)} \]

Estimasi Parameter

Metode estimasi parameter pada GLM umumnya menggunakan Maximum Likelihood Estimation (MLE).

Log-likelihood fungsi untuk regresi logistik:

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

Dengan:

\[ \pi_i = \frac{\exp(x_i^\top \beta)}{1 + \exp(x_i^\top \beta)} \]

Estimasi dilakukan melalui iterasi Newton-Raphson atau algoritma Fisher Scoring.

9.3 Model Regresi Poisson

Regresi Poisson digunakan ketika variabel respons adalah data cacah (count data), yaitu bilangan bulat non-negatif. Model ini merupakan bagian dari Generalized Linear Model (GLM) dengan asumsi bahwa distribusi variabel respons adalah distribusi Poisson.

Distribusi Poisson memilki fungsi probabilitas :

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

Bentuk dalam format exponential family :

\[ f(y; \theta) = \exp \left\{ y \log(\lambda) - \lambda - \log(y!) \right\} \]

dengan:

• \(\theta = \log(\lambda)\)
• \(b(\theta) = e^{\theta} = \lambda\)
• \(\phi = 1\)
• \(c(y, \phi) = -\log(y!)\)

Maka distribusi Poisson termasuk dalam exponential family

Fungsi Link

Fungsi link kanonik untuk distribusi Poisson adalah fungsi logaritma:

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

Sehingga modelnya menjadi:

\[ \log(\mu_i) = x_i^T \beta \]

dan fungsi inverse link:

\[ \mu_i = \exp(x_i^T \beta) \]

Estimasi Parameter

Estimasi parameter \(\beta\) dilakukan dengan metode Maximum Likelihood Estimation (MLE).
Log-likelihood fungsi untuk regresi Poisson:

\[ l(\beta) = \sum_{i=1}^{n} \left[ y_i x_i^T \beta - \exp(x_i^T \beta) - \log(y_i!) \right] \]

Nilai \(\beta\) dapat diperoleh melalui metode numerik seperti iterasi Newton-Raphson.

10 INFERENSI GLM

Dalam Generalized Linear Model (GLM), inferensi statistik membutuhkan pemahaman terhadap ekspektasi dan varians dari estimator model, terutama untuk mengembangkan alat-alat uji seperti Wald test, Likelihood Ratio test, dan interval kepercayaan.

1. Ekspektasi Estimator

Ekspektasi menunjukkan apakah suatu estimator tak bias, yaitu:

\[ \mathbb{E}[\hat{\beta}] = \beta \]

Dalam GLM, MLE dari \(\hat{\beta}\) bersifat asymptotically unbiased.

2. Varians Estimator

Varians menunjukkan presisi dari estimasi parameter:

\[ \text{Var}(\hat{\beta}) \approx (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \]

di mana \(\mathbf{W}\) adalah matriks bobot yang tergantung pada distribusi dan fungsi link.

Distribusi Asimptotik Estimator

dengan ukuran sampel besar :

\[ \hat{\beta} \sim \mathcal{N}(\beta, \text{Var}(\hat{\beta})) \]

Distribusi ini adalah dasar dari:

• Uji Wald
• Confidence interval
• P-value

Varians dalam GLM Tidak Konstan

Tidak seperti regresi linear (OLS) yang mengasumsikan homoskedastisitas:

\[ \text{Var}(Y_i) = \sigma^2 \]

Dalam GLM:

\[ \text{Var}(Y_i) = \phi V(\mu_i) \]

• \(\phi\) = parameter dispersi
  
• \(V(\mu)\) = fungsi varians

Contoh:

• Poisson: \(V(\mu) = \mu\)
• Binomial: \(V(\mu) = \mu(1 - \mu)\)

10.1 Mencari Ekspektasi dan Varians dalam GLM

10.1.1 Ekspektasi

Jika diturunkan berdasarkan fungsi momen:

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

Untuk keluarga eksponensial:

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

atau:

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

Maka ekspektasi turunan pertama:

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

Dan ekspektasi dari turunan pertama:

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

Maka:

\[ \mu = b'(\theta) \]

10.2 Varians

Turunan kedua:

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

Sehingga:

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

10.3 Metode Penaksiran Parameter

10.3.1 Maximum Likelihood Estimation (MLE)

• Prinsip dasar: memaksimumkan fungsi likelihood/log-likelihood.
• Langkah:
  • Turunan pertama = 0
  • Turunan kedua < 0

Namun, karena bentuk GLM tidak eksplisit, digunakan metode numerik.

10.3.2 Metode Optimisasi

10.3.2.1 Newton-Raphson

Menggunakan score vector (gradien) dan Hessian matrix

Iterasi :

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

10.3.2.2 Fisher Scoring

Modifikasi Newton-Raphson, mengganti Hessian dengan matriks informasi Fisher

10.3.2.3 IRLS (Iteratively Reweighted Least Square)

Modifikasi Fisher scoring, hasil estimasi mirip dengan Least Square

Implementasi Newton-Raphson

Statistik score ke-j:

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

Turunan kedua:

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

Taylor expansion:

\[ U(\hat{\beta}) \approx U(\beta) + H(\beta)(\hat{\beta} - \beta) \]

Estimasi parameter:

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

10.4 Diagnostik Model GLM

Diagnostik digunakan untuk mengevaluasi apakah model sudah tepat. Diagnostik ini bisa dengan dua cara yaitu :

• Uji formal
• Grafik antara nilai prediksi vs nilai aktual

10.4.1 Statistik Devians

• Mengukur apakah ada model lain yang lebih baik.
• Nilai devians besar → model tidak cocok.
• Devians adalah:

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

• Devians membandingkan model terhadap saturated model.
• Devians kecil → model lebih cocok.

10.4.2 Statistik Chi-Kuadrat Pearson

• Menguji apakah model lebih baik daripada tidak ada model sama sekali.
• Statistik:

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

• Jika signifikan → model lebih baik daripada tanpa model.

Catatan :

• Untuk data yang dikelompokkan, statistik devians dan chi-kuadrat Pearson mengikuti distribusi Chi-Square.
• Untuk data tidak dikelompokkan, tidak mengikuti distribusi Chi-Square.
• Devians diminimalkan oleh MLE → cocok digunakan untuk evaluasi model.

10.4.3 Analisis Residual

• Residual adalah selisih antara observasi dengan prediksi.
• Dapat digunakan untuk memeriksa penyimpangan sistematis.
• Dapat diplot untuk menilai asumsi model.

10.5 Detail Metode estimasi dan Inferensi Regresi Logistik

Regresi logistik digunakan untuk memodelkan probabilitas dari variabel respons biner (0/1) berdasarkan satu atau lebih variabel prediktor. Estimasi parameter dilakukan menggunakan Maximum Likelihood Estimation (MLE) karena model tidak linear dalam parameternya.

Fungsi model logistik :

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

Log-likelihood untuk n observasi :

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

10.5.1 Estimasi dengan Newton Raphson

Metode Newton-Raphson digunakan untuk mencari nilai parameter \(\beta\) yang memaksimalkan fungsi log-likelihood pada model regresi logistik.

Fungsi model logistik :

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

Log-likelihood untuk n observasi :

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

10.5.2 Langkah-langkah Newton Raphson

1. Turunan Pertama (Score Function)

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

2. Turunan Kedua (Hessian Matrix)

\[ \mathbf{H}(\beta) = -\mathbf{X}^\top \mathbf{W} \mathbf{X}, \quad \text{dengan } \mathbf{W} = \text{diag}(\pi_i (1 - \pi_i)) \]

3. Iterasi Newton-Raphson

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

10.5.3 Inferensi Parameter

1. Uji Wald

Tujuan Uji Wald

Untuk menguji signifikansi parameter \(\beta_j\) dalam model regresi logistik:

• \(\mathbf{H_0} : \beta_j = 0\)
• \(\mathbf{H_1} : \beta_j \neq 0\)

Teori Wald Test

Dari teori estimasi MLE, estimator \(\hat{\beta}_j\) mendekati distribusi normal:

\[ \hat{\beta}_j \sim N(\beta_j, \text{Var}(\hat{\beta}_j)) \]

Jika \(\mathbf{H_0}\) benar (yaitu \(\beta_j = 0\)), maka:

\[ Z = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \sim N(0,1) \]

Dengan statistik Wald:

\[ W = Z^2 = \left( \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \right)^2 \sim \chi^2_1 \]

Kesimpulan Uji Wald didasarkan pada rasio antara estimasi parameter dan standar error-nya. Dengan menaikkan nilai Z menjadi kuadrat (Z²), kita memperoleh distribusi Chi-Square untuk pengujian hipotesis parameter individual dalam model regresi logistik.

2. Uji Likelihood Ratio (Chi-Square)

Membandingkan model penuh dengan model tanpa prediktor.

10.5.4 Evaluasi Kebaikan Model

1. Akaike INformation Criterion (AIC)

Semakin kecil AIC, maka semakin baik model

2. Bayesian Information Criterion (BIC)

Alternatif terhadap AIC, menghukum kompleksitas model

10.6 Detail Metode Estimasi dan Inferensi Regresi Poisson

Model regresi Poisson digunakan untuk memodelkan data count (jumlah kejadian) di mana variabel respons mengikuti distribusi Poisson. Estimasi dilakukan dengan Maximum Likelihood Estimation (MLE), dan inferensi dilakukan dengan uji Wald dan Likelihood Ratio Test.

Model Regresi Poisson :

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

Fungsi distribusi poisson :

\[ \log(\lambda_i) = \mathbf{x}_i^T \beta \]

Estimasi Parameter (MLE)

Log-likelihood fungsi:

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

Dengan:

\[ \lambda_i = \exp(\mathbf{x}_i^T \beta) \]

10.6.1 Estimasi dilakukan dengan metode iterasi (IRLS)

Estimasi dilakukan dengan metode iterasi (IRLS)

Estimasi parameter model regresi Poisson menggunakan pendekatan Maximum Likelihood Estimation (MLE) dengan metode Iteratively Reweighted Least Squares (IRLS) secara manual, tanpa menggunakan glm().

Tahap 1: Definisikan Model regresi Poisson:

\[ \log(\lambda_i) = \mathbf{x}_i^T \beta \quad \text{sehingga} \quad \lambda_i = \exp(\mathbf{x}_i^T \beta) \]

Tahap 2: Mencari Log-likelihood yang dimaksimumkan:

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

Tahap 3: Formulasi iteratif:

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

Dengan:

\[ \mathbf{W} = \text{diag}(\lambda_i) \]

\[ \mathbf{z} = \eta + \frac{y - \lambda}{\lambda} \]

dan

\[ \eta_i = \log(\lambda_i) = \mathbf{x}_i^T \beta \]

untuk mengevaluasi model terbaik tetap digunakan AIC dan BIC

11 REGRESI LOGISTIK

Regresi logistik adalah salah satu metode statistik untuk memodelkan hubungan antara satu atau lebih variabel prediktor dengan variabel respon yang bersifat biner (dua kategori). Variabel prediktor ini sendiri memiliki beberapa skala pengukuran, yaitu :

• Nominal : merupakan variabel yang tidak memiliki urutan antar kategorinya. Pengukuran ini hanya membuat kategori menjadi pembeda. Contoh : jenis kelamin (laki-laki dan perempuan). Dalam konteks regresi logistik, jenis kelamin ini diubah menjadi data nominal dengan variabel dummy (1 dan 0)

• Ordinal : merupakan variabel yang memiliki urutan antar kategorinya, tapi jarak antar kategori belum tentu memiliki jarak yang sama. Contoh : tingkat pendidikan (SD < SMP < SMA < Sarjana). Dalam regresi logistik, variabel ordinal ini dapat diperlakukan :

  • Sebagai nominal dengan membuat dummy untuk tiap kategori.

  • Dijadikan rasio dengan memberi nilai tertentu yang mencerminkan urutan dan jarak antar kategorinya (misal SD = 1, SMP = 2, SMP = 3, Sarjana = 4).

• Rasio : merupakan variabel numerik kontinu yang memilki nilai nol absolut dan rasio bermakna. Contoh : pendapatan (dalam juta rupiah). Untuk variabel rasio dapat digunakan langsung tanpa perlakuan khusus.

11.1 Simulasi Data

untuk mencoba memodelkan regresi logistik dilakukan simulasi data dengan 500 observasi. Data akan berisikan 500 observasi dengan variabel gender, pendidikan, income, dan status keberhasilan

n <- 500 

gender <- sample(c("Male","Female"), n, replace = TRUE)
education <- sample(c("HighSchool", "Bachelor", "Master", "PhD"), n, replace = TRUE, prob = c (0.4,0.3,0.2,0.1))
income <- rnorm(n, mean = 50, sd = 15)

logit_p <- -2 + 0.5*(gender == "Female") + 0.8*as.numeric(factor(education, ordered=TRUE)) + 0.03*income 
p <- 1/(1+exp(-logit_p))

set.seed(123)
success <- rbinom(n,1,p)

library(tibble)
sim_data <- tibble(success, gender, education, income)

head(sim_data)

## # A tibble: 6 × 4
##   success gender education  income
##     <int> <chr>  <chr>       <dbl>
## 1       1 Female Bachelor     26.6
## 2       1 Female HighSchool   54.9
## 3       1 Male   Master       47.6
## 4       0 Female HighSchool   63.2
## 5       0 Male   Bachelor     61.3
## 6       1 Female Bachelor     54.5

11.2 Eksplorasi Data

sim_data %>% 
  dplyr::group_by(success) %>% 
  dplyr::summarise(
    Jumlah = dplyr::n(), 
    Rata2_Income = mean(income)
  )

## # A tibble: 2 × 3
##   success Jumlah Rata2_Income
##     <int>  <int>        <dbl>
## 1       0    112         43.4
## 2       1    388         50.9

Dari hasil eksplorasi data, distribusi jumlah sukses dan tidak sukses adalah seperti yang ada pada tabel. Pada kolom ke 3 adalah rata-rata pendapatan dari masing-masing grup.

11.3 Perlakuan Variabel Ordinal

11.3.1 Diperlakukan sebagai Nominal (Dummy)

Dari variabel ordinal yang ada (tingkat pendidikan), variabel tersebut diberlakukan sebagai variabel dummy yang kemudian dibuat model regresi logistiknya.

sim_data_nominal <- sim_data %>%
mutate(
  education = factor(education, levels = c("HighSchool", "Bachelor", "Master", "PhD"))
)
model_nominal <- glm(success ~ gender + education + income, data = sim_data_nominal, family = binomial)
summary(model_nominal) 

## 
## Call:
## glm(formula = success ~ gender + education + income, family = binomial, 
##     data = sim_data_nominal)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.057084   0.416643   0.137 0.891023    
## genderMale        -0.818238   0.234399  -3.491 0.000482 ***
## educationBachelor -0.789155   0.256701  -3.074 0.002111 ** 
## educationMaster    0.211401   0.335221   0.631 0.528282    
## educationPhD       1.712471   0.635252   2.696 0.007023 ** 
## income             0.037117   0.007902   4.697 2.64e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 531.92  on 499  degrees of freedom
## Residual deviance: 470.04  on 494  degrees of freedom
## AIC: 482.04
## 
## Number of Fisher Scoring iterations: 5

exp(coef(model_nominal))

##       (Intercept)        genderMale educationBachelor   educationMaster 
##         1.0587451         0.4412083         0.4542283         1.2354078 
##      educationPhD            income 
##         5.5426384         1.0378145

Model menggunakan education sebagai variabel dummy dengan baseline “High School”. Setiap koefisien membandingkan kategori terhadap baseline

Interpretasi Koefisien :

Intercept (0.328351)

• Merupakan log-odds dasar untuk individu berjenis kelamin Female, berpendidikan High School, dan memiliki income = 0

• Nilai odds rationya adalah 1.39 yang menunjukkan peluang sukses lebih tinggi dibandingkan baseline. Namun hal ini tidak signifikan karena p-value lebih besar dibandingkan 0.05

Gender Male (-0.852196)

• Individu male memiliki log-odds sukses yang lebih rendah sebesar 0.85 dibandingkan dengan Female (baseline).

• Dengan p-value lebih kecil dibandingkan 5% menunjukkan bahwa gender ini signifikan dibandingkan dengan Female.

• Nilai odds rationya adalah 0.43 yang menunjukkan adanya peluang 43% seorang dengan gender Male untuk sukses dibandingkan gender Female.

Education Bachelor (-0.856816)

• Individu dengan pendidikan Bachelor memiliki log-odds sukses yang lebih rendah sebesar 0.86 dibandingkan dengan HighSchool (baseline)

• Hal ini signifikan dan menunjukkan bahwa perbedaan tingkat pendidikan menjadi hal yang signifikan dibandingkan dengan baseline

• Nilai odds rationya adalah 0.42 yang menunjukkan bahwa hanya 42% seorang dengan pendidikan Bachelor sukses dibandingkan dengan pendidikan HighSchool

Education Master (1.047507)

• Individu dengan pendidikan Master memiliki log-odds sukses lebih tinggi sebesar 1.05 dibandingkan dengan HighSchool dan hal ini signifikan

• Nilai odds rationya sebesar 2.85 yang menunjukkan bahwa peluang seseorang dengan pendidikan master untuk sukses 2 kali lipat lebih tinggi dibandingkan education highschool

Education PhD (1.169522)

• Individu dengan pendidikan PhD memiliki log-odds subses lebih tinggi sebesar 1.17 dibandingkan dengan HighSchool tetapi hal ini tidak signifikan

• Nilai odds rationya sebesar 3.22 yang menunjukkan bahwa peluang individu dengan pendidikan PhD untuk sukses adalah 3 kali lipat dibandingkan dengan education high school tetapi hal ini tidak signifikan.

Income (0.030674)

• Setiap kenaikan 1 unit income maka akan menaikkan log-oddsnya sebesar 0.031 dan hal ini signifikan.

• Nilai odds rationya sebesar 1.031 yang artinya setiap kenaikan 1 unit income maka akan meningkatkan 3.1% peluang sukses.

Goodness of Fit

• Dari hasil nilai null deviance dan residual deviance ada penurunan nilai yang menunjukkan bahwa model ini memiliki informasi yang cukup baik.

• Nilai AIC sebesar 468.09 digunakan untuk membandingkan model dengan model lain. Semakin kecil nilai AIC maka semakin baik pula model tersebut.

Signifikansi Model

• Variabel income dan pendidikan (Bachelor dan Master) menjadi prediktor penting dalam peluang sukses.

• Variabel gender juga berpotensi berpengaruh tetapi memerlukan variabel tambahan.

Kesimpulan

• Variabel income dan pendidikan yang semakin tinggi menjadi prediktor penting dalam menentukan peluang sukses.

• Variabel gender berpengaruh tetapi memerlukan informasi tambahan dari variabel lain.

• Model cukup baik dibandingkan model null, tetapi masih dapat ditingkatkan dengan penambahan variabel.

11.3.2 Diperlakukan sebagai Rasio (Numeric Rank)

Variabel pendidikan diubah menjadi rank dengan nilai 1 pada HighSchool, 2 pada Bachelor, 3 pada Master, dan 4 pada PhD.

sim_data_numeric <- sim_data %>%
mutate(
education_numeric = case_when(
education == "HighSchool" ~ 1,
education == "Bachelor" ~ 2,
education == "Master" ~ 3,
education == "PhD" ~ 4
)
)
model_numeric <- glm(success ~ gender + education_numeric + income, data = sim_data_numeric, family = binomial)
summary(model_numeric)

## 
## Call:
## glm(formula = success ~ gender + education_numeric + income, 
##     family = binomial, data = sim_data_numeric)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.456882   0.447097  -1.022 0.306835    
## genderMale        -0.792548   0.228495  -3.469 0.000523 ***
## education_numeric  0.257599   0.115044   2.239 0.025147 *  
## income             0.034871   0.007546   4.621 3.82e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 531.92  on 499  degrees of freedom
## Residual deviance: 493.76  on 496  degrees of freedom
## AIC: 501.76
## 
## Number of Fisher Scoring iterations: 4

exp(coef(model_numeric))

##       (Intercept)        genderMale education_numeric            income 
##         0.6332552         0.4526899         1.2938205         1.0354859

Interpretasi Koefisien :

Intercept (-0.333852)

• Merupakan log-odds dasar untuk gender perempuan, pendidikan high school (education_numeric = 1) dan income = 0

• Memiliki nilai odds ratio 0.72 yang menunjukkan bahwa peluang dasar lebih kecil dari baseline tetapi tidak signifikan karena p-valuenya lebih besar dibandingkan taraf.

Gender Male (-0.820981)

• Individu male memiliki log-odds lebih rendah 0.82 dibandingkan dengan female

• Dengan nilai odds ratio 0.44 menunjukkan bahwa laki-laki memiliki peluang sukses hanya 44% dibandingkan dengan perempuan dan hal ini signifikan

Education_Numeric (0.335171)

• Setiap kenaikan tingkat pendidikan akan meningkatkan nilai log-oddsnya sebesar 0.34. Hal ini signifikan menunjukkan bahwa semakin tinggi pendidikan semakin tinggi juga peluang suksesnya.

• Nilai odds ratio sebesar 1.40 menunjukkan bahwa setiap kenaikan 1 tingkat pendidikan akan meningkatkan 40% peluang sukses.

Income (0.029346)

• Setiap kenaikan 1 satuan pendapatan (juta rupiah) akan meingkatkan nilai log-odds sebesar 0.03. Hal ini signifikan menunjukkan bahwa semakin tinggi pendapatan seseorang semakin tinggi juga peluang suksesnya.

• Nilai odds rationya sebesar 1.03 yang menunjukkan setiap kenaikan 1 satuan pendapatan akan meningkatkan 3% peluang sukses.

Goodness of Fit

• Dari hasil nilai null deviance dan residual deviance ada penurunan nilai yang menunjukkan bahwa model ini memiliki informasi yang cukup baik.

• Nilai AIC sebesar 495.09 digunakan untuk membandingkan model, semakin kecil nilai AIC maka semakin baik pula model tersebut.

Signifikansi Model

• Variabel income dan education_numeric menjadi variabel prediktor penting dalam penentuan peluang sukses

• Variabel gender menjadi salah satu variabel penting tetapi tidak bisa berdiri sendiri dan membutuhkan variabel-variabel lain

Kesimpulan

• Semakin tinggi tingkat pendidikan maka akan semakin tinggi peluang sukses

• Semakin tinggi income maka semakin tinggi pula peluang suksesnya

• Gender perempuan memiliki peluang sukses yang lebih tinggi dibandingkan gender laki-laki

Perbandingan Model

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

## $AIC_Nominal
## [1] 482.0377
## 
## $AIC_Numeric
## [1] 501.7551

nullmod <- glm(success ~ 1, data = sim_data, family = binomial)
r2_nominal <- 1 - (logLik(model_nominal)/logLik(nullmod))
r2_numeric <- 1 - (logLik(model_numeric)/logLik(nullmod))
list(
McFadden_R2_Nominal = r2_nominal,
McFadden_R2_Numeric = r2_numeric
)

## $McFadden_R2_Nominal
## 'log Lik.' 0.1163446 (df=6)
## 
## $McFadden_R2_Numeric
## 'log Lik.' 0.0717566 (df=4)

Dengan membandingkan nilai AIC dan McFadden R-Squared didapatkan bahwa model nominal lebih baik dibandingkan dengan model numeric karena nilai AIC yang lebih kecil dan nilai R-squared yang lebih besar sehingga dapat menjelaskan model dengan lebih baik.

11.3.3 Visualisasi Model

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

sim_data_nominal %>%
  ggplot(aes(x = income, y = predicted, color = education)) +
  geom_point(alpha = 0.6) +
  labs(
    title = "Prediksi Probabilitas Sukses Berdasarkan Income dan Pendidikan",
    x = "Income",
    y = "Probabilitas Prediksi",
    color = "Pendidikan"
  ) +
  theme_minimal()

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

sim_data_numeric %>%
  ggplot(aes(x = income, y = predicted, color = as.factor(education_numeric))) +
  geom_point(alpha = 0.6) +
  labs(
    title = "Prediksi Probabilitas (Ordinal sebagai Numeric)",
    x = "Income",
    y = "Probabilitas Prediksi",
    color = "Education Level (Numeric)"
  ) +
  theme_minimal()

11.3.4 Ringkasan

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

# Ringkasan model nominal
tidy(model_nominal) %>%
  kable(format = "latex", booktabs = TRUE,
        caption = "Ringkasan Koefisien Model dengan Education sebagai Variabel Nominal") %>%
  kable_styling(latex_options = c("hold_position", "striped"))

# Ringkasan model numeric
tidy(model_numeric) %>%
  kable(format = "latex", booktabs = TRUE,
        caption = "Ringkasan Koefisien Model dengan Education sebagai Variabel Numeric") %>%
  kable_styling(latex_options = c("hold_position", "striped"))

Dari hasil pemodelan,

• Gender perempuan memiliki peluang sukses yang lebih tinggi dibandingkan dnegan gender laki-laki

• Semakin tinggi pendidikan, semakin besar juga peluang untuk sukses. JIka variabel pendidikan dianggap sebagai dummy, variabel HighSchool menjadi baseline, jika dianggap sebagai numeric maka setiap kenaikan 1 tingkat akan meningkatkan peluang sukses.

• Income, semakin tinggi pendapatan individu maka semakin tinggi juga peluang suksesnya.

12 PEMILIHAN MODEL REGRESI LOGISTIK DAN EVALUASI

12.1 Membangun Model Regresi Logistik : Pendekatan Confirmatory dan Exploratory

Dalam analisis regresi logistik, pemilihan model menjadi penting untuk mendapatkan model yang baik dalam memprediksi probabilitas kejadian suatu peristiwa (respon biner). Pendekatan yang digunakan untuk membangun model adalah pendekatan Confirmatory dan Exploratory.

12.1.1 Pendekatan Confirmatory

Pendekatan ini digunakan ketika peneliti telah memiliki teori atau hipotesis yang jelas mengenai efek atau hubungan antara variabel prediktor dan respon

Ciri-ciri :

• Model dibangun berdasarkan teori atau hasil penelitian sebelumnya

• Tujuan utama membangun model dengan pendekatan ini tidak hanya mencari model terbaik tetapi juga menguji apakah efek terkait signifikan.

• Peneliti biasa menyusun model penuh terlebih dahulu kemudian dilakukan pengujian dengan penambahan atau pengurangan efek memberikan peningkatan model secara signifikan atau tidak.

• Uji signifikansi dilakukan dengan membandingkan model dengan efek tertentu dan model tanpa efek tersebut, misalnya dengan Likelihood Ratio Test.

Contoh penggunaan :

Misal, teori menyatakan bahwa faktor x1 dan x2 mempengaruhi probabilitaas seseorang membeli produk. Maka model logistik langsung dibangun dengan x1 dan x2 lalu diuji apakah kontribusi x2 benar-benar signifikan.

12.1.2 Pendekatan Exploratory

Pendekatan ini digunakan ketika peneliti belum memiliki teori yang pasti atau ingin mengeksplorasi hubungan potensial antar variabel.

Ciri-ciri :

• Model dibangun secara bertahap dengan tujuan menemukan kombinasi prediktor terbaik

• Pemilihan variabel dilakukan berdasarkan kriteria statistik, seperti AIC, deviance, atau log-llikelihood.

12.1.3 Tujuan

Pemilihan model terbaik dilakukan melalui :

• Forward Selection : mulai dari model kosong, satu per satu variabel dimasukkan jika signifikan

• Backward Elimination : mulai dari model penuh, variabel yang tidak signifikan dikeluarkan

• Stepwise Selection : penggabungan dari kedua metode, variabel dapat dimasukkan atau dikeluarkan secara dinamis.

12.2 Metode Stepwise : Forward, Backward, dan Kedua Arah

Dilakukan simulasi data dengan jumlah 200 data

library(knitr)
library(dplyr)
library(ggplot2)
library(MASS)

## 
## Attaching package: 'MASS'

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

library(caret)

## Loading required package: lattice

## 
## Attaching package: 'caret'

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

library(pROC)

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

## 
## Attaching package: 'pROC'

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

library(DescTools)
set.seed(123)

n <- 200
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)
lin_pred <- -0.5 + 1.2 * x1 - 0.8 * x2 + 0.5 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
df <- data.frame(y = as.factor(y), x1, x2, x3)
head(df)

##   y          x1 x2         x3
## 1 0 -0.56047565  1 -0.7152422
## 2 0 -0.23017749  0 -0.7526890
## 3 1  1.55870831  1 -0.9385387
## 4 1  0.07050839  1 -1.0525133
## 5 1  0.12928774  0 -0.4371595
## 6 1  1.71506499  0  0.3311792

Data simulasi ini memiliki 3 variabel prediktor yaitu x1 (ratio), x2 variabel nominal yang dijadikan variabel dummy, dan x3 (ratio) Dari data simulasi, akan dibuat model regresi logistiknya.

#Model Full
model_full <- glm(y ~ x1 + x2 + x3, data = df, family = binomial)
summary(model_full)

## 
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.7148     0.2470  -2.894  0.00381 ** 
## x1            1.4029     0.2315   6.061 1.35e-09 ***
## x2           -0.2507     0.3463  -0.724  0.46903    
## x3            0.3567     0.1704   2.094  0.03630 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 257.72  on 199  degrees of freedom
## Residual deviance: 202.67  on 196  degrees of freedom
## AIC: 210.67
## 
## Number of Fisher Scoring iterations: 4

Setelah didapatkan modelnya akan dilakukan pemilihan modelnya dengan stepwise dan membandingkan nilai AIC nya.

null_model <- glm(y ~ 1, data = df, family = binomial)
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)
AIC(model_full, step_forward, step_backward, step_both)

##               df      AIC
## model_full     4 210.6739
## step_forward   3 209.1998
## step_backward  3 209.1998
## step_both      3 209.1998

Dari hasil stepwise didapatkan bahwa nilai AIC setalah dilakukan forward, backward, atau keduanya memiliki AIC yang lebih kecil dibandingkan dnegan model full sehingga dilakukan evaluasi model.

12.3 Prinsip Parsimony

Model yang kompleks sering memiliki AIC dan deviance yang lebih kecil, tetapi model sederhana lebih mudah diinterpretasikan, Jika penurunan AIC tidak signifikan, pilih model yang lebih sederhana. Model yang sedeerhana atau parsimony mencegah adanya overfitting.

12.3.1 AIC

AIC adalah ukuran untuk menilai model berdasarkan kombinasi antara goodness of fit (melalui kompleksitas (melalui jumlah parameter k. Semakin kecil nilai AIC, maka semakin baik model tersebut karena AIC menghukum model yang terlalu kompleks meskipun memiliki likelihood tinggi. Rumus AIC adalah sebagai berikut :

\[ \text{AIC} = -2(log L - k ) = -2logL + 2k \]

12.3.2 Deviance

Deviance mengukur seberapa jauh model saat ini dibandingkan dengan model sempurna (model saturasi). Nilai deviance yang kecil menunjukkan model memberikan prediksi yang mendekati nilai aktual. Rumus Deviance adalah sebagai berikut :

\[ D = -2[logL(\text{model}) - logL(\text{model saturasi}) \]

12.3.3 Likelihood Ratio

Statistik Likelihood Ratio digunakan untuk menguji apakah penambahan variabel dalam model secara signifikan meningkatkan kecocokan model. Jika \(G^{2}\) besar dan p-value kecil, maka model kompleks lebih baik dari model sederhana secara statistik. Rumus likelihood ratio adalah sebagai berikut :

\[ G^{2} = -2(logL_{0} - logL_{1}) \]

12.4 Evaluasi Model : ROC dan AUC

ROC curve atau Receiever Operating Characteristic Curve adalah grafik yang menggambarkan performa suatu model di berbagai threshold. Kurva ini menunjukkan trade-off antara True Positive Rate dan False Positive Rate pada berbagai threshold klasifikasi.

1. Definisi

• Sumbu Y :

  \[ \text{Sensitivity} = \text{True Postitive Rate} = \frac{TP}{(TP+FN)} \]

• Sumbu X :

  \[ \text{1-specificity} = \text{False Positive Rate} = \frac{FP}{FP+TN} \]

• Garis diagonal (dari kiri bawah ke kanan atas) menunjukkan performa acak (random guess).

• Kurva yang mendekati pojok kiri ata menunjukkan performa klasifikasi yang lebih baik.

2. Cut-off dan Pergerakan Kurva

Saat cut-off menurun, model mengklasifikasikan lebih banyak pengamatan sebagai positif. Sensitivitas naik dan spesifisitas turun. Saat cut-off meningkat, model menjadi lebih konservatif dimana sensitivitas turun dan spesifisitas naik.

3. Kurva ROC Ideal

Kurva ROC yang ideal memiliki bentuk :

• Naik tajam secara vertikal hingga mencapai sensitivitas = 1

• Bergerak secara horizontal menuju 1-specificity = 1

• Area Under the Curve (AUC) mendekati 1

Untuk mendapatkan kurva yang ideal, AUC atau Area Under Curve yang merupakan luas di bawa kurva ROC. Nilai AUC berada di antara 0 dan 1. Semakin dekat nilai AUC dengan 1 maka model akan semakin sempurna dan semakin kecil nilai AUC nya (<0.5) berarti model tersebut buruk atau tidak lebih baik dari tebak acak.

4. Kegunaan Kurva ROC

• Membandingkan performa beberapa model klasifikasi

• Memilih threshold (cut-off) optimal berdasarkan kebutuhan aplikasi (misalnya : lebih penting menghindari false negative atau false positive)

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.7964

5. Pemilihan Threshold Optimal

Pemiliihan threshold terbaik dapat dengan mengevaluasi sensitivitas dan spesifisitas pada berbagai cut-off

thresholds <- seq(0.1, 0.9, by = 0.05)
results <- data.frame(Threshold = thresholds)

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

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

print(results)

##    Threshold Sensitivity Specificity
## 1       0.10  0.98550725   0.2137405
## 2       0.15  0.92753623   0.3358779
## 3       0.20  0.84057971   0.5114504
## 4       0.25  0.81159420   0.5877863
## 5       0.30  0.75362319   0.6870229
## 6       0.35  0.69565217   0.7633588
## 7       0.40  0.68115942   0.8091603
## 8       0.45  0.62318841   0.8396947
## 9       0.50  0.53623188   0.8854962
## 10      0.55  0.47826087   0.9083969
## 11      0.60  0.42028986   0.9465649
## 12      0.65  0.33333333   0.9465649
## 13      0.70  0.30434783   0.9694656
## 14      0.75  0.18840580   0.9847328
## 15      0.80  0.13043478   0.9847328
## 16      0.85  0.05797101   0.9847328
## 17      0.90  0.02898551   1.0000000

Cut-off yang optimal dipilih berdasarkan nilai maksimum dari sensitivity + specificity atau mempertimbangkan trade-off sesuai tujuan aplikasi.

12.5 Precision-Recall Curve (PR Curve)

Kurva Precision-Recall (PR) adalah alat evaluasi performa model klasifikasi seperti ROC tetapi PR Curve ini sangat berguna saat bekerja dengan data yang tidak seimbang (class imbalance)

1. Definisi

• Precision (Presisi) merupakan proposi prediksi positif yang benar-benar positif

\[ \text{Precision} = \frac{TP}{TP+FP} \]

• Recall (Sensitivitas) merupakan proporsi kasus positif yang berhasil diprediksi positif

\[ \text{Recall} = \frac{TP}{TP+FN} \]

2. Interpretasi

• PR Curve menunjukkan bagaimana presisis berubah saat recall meningkat.

• Idealnya, kita ingin nilai presisi dan recall keduanya tinggi, tetapi biasanya ada trade-off

• MOdel dengan performa baik memiliki PR Curve yang melengkung ke pojok kanan atas

3. Area Under PR Curve

• Luas kurva (AUPRC) mendekati 1 menandakan bahwa model sangat baik

• BAseline AUPRC = prevalensi kelas positif dalam data

4. Perbandingan ROC dan PRC

\[ \begin{array}{ccc} \hline \textbf{Aspek} & \textbf{ROC Curve} & \textbf{Precision-Recall Curve}\\ \hline \textbf{Fokus} & \text{Semua Kelas} & \text{Kelas Positif Saja} \\ \textbf{Kuat di} & \text{Data Seimbang} & \text{Data Tidak Seimbang}\\ \textbf{Sumbu Y} & \text{Sensitivitas (Recall)} & \text{Precision} \\ \textbf{Sumbu X} & \text{1-Spesifisitas} & \text{Recall} \\ \hline \end{array} \]

library(PRROC)

## Loading required package: rlang

set.seed(123)
x1 <- rnorm(300)
x2 <- rbinom(300, 1, 0.5)
x3 <- rnorm(300)
lin_pred <- -1 + 1.2 * x1 - 0.6 * x2 + 0.8 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(300, 1, p)
data <- data.frame(y = y, x1, x2, x3)
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
prob <- predict(model, type = "response")
pr <- pr.curve(scores.class0 = prob[data$y == 1],
scores.class1 = prob[data$y == 0],
curve = TRUE)
plot(pr)

12.6 Pseudo R-Squared

Dalam model regresi linear, digunakan R-squared untuk melihat seberapa besar variasi dalam data yang bisa dijelaskan oleh model, tetapi dalam regresi logistik hal ini tidak dapat digunakan karena model tidak memodelkan nilai kontinu dan varians residualnya tidak terdefinisi dengan cara yang sama seperti regresi linear.

Untuk regresi logisitik digunakan psuedo R-Squared dengan McFadden’s R-squared, Cox & Snell, dan Nagelkerke (adjusted Cox-Snell)

Rumusnya adalah sebagai berikut :

\[ R^2_{\text{Cox and Snell}} = 1 - \left( \frac{L_0}{L_M} \right)^{2/n} \]

\[ R^2_{\text{McFadden}} = 1 - \frac{\log L_M}{\log L_0} \]

Dengan keterangan :

• \(L_0\) : log-likelihood dari model null (tanpa prediktor)
• \(L_ M\): log-likelihood dari model penuh
• \(n\) : jumlah observasi

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

##   CoxSnell Nagelkerke   McFadden 
##  0.2385981  0.3294000  0.2115439

Dengan hasil nilai Psuedo R-Squared ada beberapa hal yang ditunjukkan, dengan nilai McFadden’s R-squared sebesar 0.21 model sudah termasuk dalam kategori model yang baik karena interpretasi dari McFadden ini lebih longgar dan ketika nilainya berada di antara 0.2 sampai 0.4 menandakan model sudah baik. Selain itu dengan nilai Nagelkerke R-Squared 0.33 menandakan bahwa model mampu menjelaskan sepertiga variasi dalam data. Dari semua informasi ini, variabel variabel dalam model memiliki kontribusi yang berarti dalam menjelaskan kemungkinan sukses.

12.7 Tabel Klasifikasi dan Evaluasi

Tabel klasifikasi atau Confusion Matrix adalah alat untuk mengevaluasi model. Tabel ini akan membandingkan prediksi dari model dengan hasil sebenarnya.

pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
conf_matrix <- confusionMatrix(factor(pred_class), df$y, positive = "1")
conf_matrix

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 116  32
##          1  15  37
##                                        
##                Accuracy : 0.765        
##                  95% CI : (0.7, 0.8219)
##     No Information Rate : 0.655        
##     P-Value [Acc > NIR] : 0.0005028    
##                                        
##                   Kappa : 0.4478       
##                                        
##  Mcnemar's Test P-Value : 0.0196041    
##                                        
##             Sensitivity : 0.5362       
##             Specificity : 0.8855       
##          Pos Pred Value : 0.7115       
##          Neg Pred Value : 0.7838       
##              Prevalence : 0.3450       
##          Detection Rate : 0.1850       
##    Detection Prevalence : 0.2600       
##       Balanced Accuracy : 0.7109       
##                                        
##        'Positive' Class : 1            
## 

conf_matrix$byClass[c("Sensitivity", "Specificity")]

## Sensitivity Specificity 
##   0.5362319   0.8854962

Dari hasil confusion matrix, menunjukkan bahwa model ini sudah baik karena sudah menunjukkan 76,5 accuracy dan model sudah cukup akurat.

12.7.1 Sensitivitas dan Spesifisitas

• Sensitivitas adalah kemampuan model mendeteksi kelas positif secara benar (True Positive Rate). Rumus untuk sensitivity adalah

  \[ \text{Sensitivity} = \frac{TP}{TP + FN} \]

• Spesifisitas adalah kemampuan model mendeteksi kelas negatif secara benar (True Negative Rate)

  \[ \text{Specificity}=\frac{TN}{TN+FP} \]

13 DISTRIBUSI MULTINOMIAL

Distribusi multinomial adalah perluasan dari distribusi ninomial dengan lebih dari dua kategori. Jika \(X_{1},X_{2}, ... , X_{k}\) menyatakan banyaknya kejadian dalam masing-masing dari k kategori, maka :

\[ P(X_{1} = x_1, ... , X_k = x_k) = \frac{n!}{x_{1}!x_{2}!...x_{k}!}p_{1}^{x_1}p_{2}^{x_2}...p_{k}^{x_k} \]

dengan \(\sum_{i=1}^{n} x_i = n \text{dan} \sum_{i=1}^{n} p_i = 1\)

13.1 Studi Kasus

Sebuah survei dilakukan terhadap 10 orang yang diminta memilih satu dari tiga jenis buah favorit : Apel (A), Jeruk (J), dan Pisang (P).

Hasil survei :

• Apel : 4 orang

• Jeruk : 3 orang

• Pisang : 3 orang

Probabilitas teoretik preferensi :

• \(p_A = 0.3\)

• \(p_J = 0.4\)

• \(p_P = 0.3\)

Pertanyaan : Berapa peluang bahwa dalam 10 orang akan ada 4 yang memilih apel, 3 memilih jeruk, dan 3 memilih pisang

Perhitungan di R

n <- 10
x <- c(4, 3, 3)
p <- c(0.3, 0.4, 0.3)
# Hitung komponen-koefisien
faktorial_total <- factorial(n)
faktorial_x <- prod(factorial(x))
koefisien <- faktorial_total / faktorial_x
# Hitung peluang
peluang <- koefisien*prod(p^x)
peluang

## [1] 0.05878656

Probabilitas bahwa 4 orang memilih apel, 3 jeruk, dan 3 pisang adalah 0.058787. Distribusi multinomial digunakan untuk menghitung peluang dalam percobaan dengan beberapa kategori hasil. Rumus dasarnya merupakan generalisasi dari binomial untuk lebih dari dua kategori

13.2 Multinomial Logistic Regression

Model ini digunakan untuk memodelkan hubungan antara satu variabel respon kategorik (lebih dari 2 kategori) dan satu atau lebih variabel prediktor.

Mislakan Y memiliki K kategori dan kita pilih referensi (baseline) kateghori K, maka model logit untuk kategori j adalah :

\[ \log\left(\frac{P(Y=j}{P(Y=K)}\right) = \beta_{j0} + \beta_{j1}x_1 + \dots + \beta_{jp}x_p \]

untuk j = 1,2,…, K-1

13.2.1 Baseline-category Logit Model

Baseline-category logit model adalah model regresi logistik untuk variabel respon kategorik dengaan lebih dari dua kategori (nominal). Model ini menggunakan satu kategori sebagai acuan (baseline) dan membandingkan kategori lainnya terhadap baseline tersebut dalam bentuk logit .

\[ \log \left (\frac{\pi_j}{\pi_c}\right), j = 1,\dots,c-1 \]

dengan :

• \(\pi_j\) adalah probabilitas espon berada di kategori j

• \(\pi_c\) adalah probabilitas respon berada di kategori acuan (baseline)

Terdapat sebanyak (c-1) fungsi logit.

Jika terdapat satu prediktor x, maka bentuk umum model logit-nya adalah :

\[ \log \left (\frac{\pi_j}{\pi_c}\right) = \alpha_j + \beta_jx, j = 1,\dots,c-1 \]

Contoh Kasus : 3 kategori respon.

Misalkan respon Y memiliki tiga kategori: \(Y \in \{1,2,3\}\) dan kita gunakan kategori ke-3 sebagai baseline. Maka :

\[ \log \left(\frac{\pi_1}{\pi_3}\right) = \alpha_1 + \beta_1x \] \[ \log \left(\frac{\pi_2}{\pi_3}\right) = \alpha_2 + \beta_2x \]

Terdapat dua model logit satu untuk membandingkan kategori 1 dengan 3 dan untuk membandingkan kategori 2 dan 3

Relasi Antar Kategori : Jika ingin menghitung logit antara kategori 1 dan 2 :

\[ \log \left(\frac{\pi_1}{\pi_3}\right) = \log\left(\frac{\frac{\pi_1}{\pi_3}}{\frac{\pi_2}{\pi_3}}\right) = \log\left(\frac{\pi_1}{\pi_3}\right) - \log\left(\frac{\pi_2}{\pi_3}\right) \]

\[ = (\alpha_1 + \beta_1x) - (\alpha_2+\beta_2x) = (\alpha_1 - \alpha_2) + (\beta_1 - \beta_2)x \] Model baseline-category logit :

• Digunakan untuk respon dengan kategori lebih dari 2

• menghasilkan (c-1) fungsi logit terhadap satu baseline

• Logit antara kategori selain baseline dapat dihitung dengan selisih dau logit terhadap baseline

13.3 Estimasi Parameter

Estimasi parameter dilakukan dengan menggunakna metode maximum likelihood dengan algoritma iteratif seperti Newton-Raphson.

Log-likelihood :

\[ \ell(\theta) = \sum_{i=1}^{n} \sum_{j=1}^{K} y_{ij} \log(\pi_j) \quad \text{dengan} \quad \pi_j = P(Y_i = j \mid x_i), \quad y_{ij} = \begin{cases} 1 & \text{jika } Y_i = j \\ 0 & \text{lainnya} \end{cases} \]

13.4 Contoh Kasus

13.4.1 Kasus 1

Sebuah perusahaan teknologi ingin memahami faktor-faktro yang mempengaruhi preferensi karyawan terhadap jenis perngkat kerja yang diberikan : Laptop, Desktop, atau Tablet. Perusahaan melakukan survei terhadap 150 karyawan dan mengumpulkan data berikut :

• Device : jenis perngkat yang dipilih (Laptop, Desktop, Tablet)

• Age : Usia karyawan

• Department : divisi tempat bekerja (IT, Marketing, HR)

• Experience : tahun pengalman kerja

Tujuannya adalah : mengetahui bagaimana usia, divisi, dan pengalaman kerja mempengaruhi pilihan perangkat kerja

13.4.1.1 Simulasi Data

set.seed(123)
n <- 150
Department <- sample(c("IT", "Marketing", "HR"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 35, sd = 7))
Experience <- round(pmax(rnorm(n, mean = 7, sd = 3), 0))
# Simulasikan Device berdasarkan probabilitas berbeda per Department
Device <- sapply(Department, function(dep) {
if (dep == "IT") {
sample(c("Laptop", "Desktop", "Tablet"), size = 1, prob = c(0.6, 0.2, 0.2))
} else if (dep == "Marketing") {
sample(c("Laptop", "Desktop", "Tablet"), size = 1, prob = c(0.3, 0.3, 0.4))
} else {
sample(c("Laptop", "Desktop", "Tablet"), size = 1, prob = c(0.4, 0.4, 0.2))
}
})
df <- data.frame(Device = factor(Device), Age, Department = factor(Department), Experience)
df$Device <- relevel(df$Device, ref = "Laptop") # baseline
head(df)

##    Device Age Department Experience
## 1 Desktop  36         HR          8
## 2 Desktop  44         HR         10
## 3  Tablet  32         HR          3
## 4  Laptop  37  Marketing          9
## 5  Laptop  34         HR         14
## 6  Tablet  38  Marketing          2

13.4.1.2 Estimasi Model

library(nnet)
model_mnlogit <- multinom(Device ~ Age + Department + Experience, data = df)

## # weights:  18 (10 variable)
## initial  value 164.791843 
## iter  10 value 153.109490
## final  value 153.057472 
## converged

summary(model_mnlogit)

## Call:
## multinom(formula = Device ~ Age + Department + Experience, data = df)
## 
## Coefficients:
##         (Intercept)         Age DepartmentIT DepartmentMarketing  Experience
## Desktop  -0.1947009  0.00396078  -1.04132876           0.2153637 -0.03133456
## Tablet   -0.2612636 -0.02768129   0.03737652           1.2451673  0.04419583
## 
## Std. Errors:
##         (Intercept)        Age DepartmentIT DepartmentMarketing Experience
## Desktop    1.139007 0.02802812    0.5359687           0.4742033 0.06782814
## Tablet     1.168694 0.02813761    0.5419639           0.4997768 0.06972093
## 
## Residual Deviance: 306.1149 
## AIC: 326.1149

13.4.1.3 Nilai p-value dan Interpretasi

z <- summary(model_mnlogit)$coefficients / summary(model_mnlogit)$standard.errors
pval <- 2 * (1 - pnorm(abs(z)))
round(pval, 4)

##         (Intercept)    Age DepartmentIT DepartmentMarketing Experience
## Desktop      0.8643 0.8876        0.052              0.6497     0.6441
## Tablet       0.8231 0.3252        0.945              0.0127     0.5261

• Koefisien untuk kategori Desktop dan Tablet dibandingkan dengan baseline Laptop

• Nilai p-value yang kecil menunjukkan bahwa variabel tersebut signifikan memengaruhi preferensi karyawan memilih perangkat

13.4.1.4 Prediksi dan Validasi

df$Predicted <- predict(model_mnlogit)
table(Predicted = df$Predicted, Actual = df$Device)

##          Actual
## Predicted Laptop Desktop Tablet
##   Laptop      50      26     21
##   Desktop      3       0      1
##   Tablet      13      15     21

13.4.1.5 Kesimpulan

Dari model regresi logistik yang terbentuk didapatkan bahwa :

• Terdapat 2 model logit yang membandingkan kategori Desktop dan Tablet terhadap Laptop

• Tidak ada faktor yang signifikan membedakan preferensi karyawan memilih dekstop dan laptop.

• karyawan dari departemen marketing secaara signifikan lebih besar kemungkinan memilih tablet dibandingkan dengan laptop.

13.4.2 Kasus 2

Digunakan dataset iris untuk memodelkan jenis spesies bunga berdasarkan lebar dan panjang kelopak

13.4.2.1 Dataset

data(iris)
iris <- iris %>% mutate(Species = relevel(Species, ref = "setosa"))
model <- multinom(Species ~ Petal.Length + Petal.Width, data = iris)

## # weights:  12 (6 variable)
## initial  value 164.791843 
## iter  10 value 12.657828
## iter  20 value 10.374056
## iter  30 value 10.330881
## iter  40 value 10.306926
## iter  50 value 10.300057
## iter  60 value 10.296452
## iter  70 value 10.294046
## iter  80 value 10.292029
## iter  90 value 10.291154
## iter 100 value 10.289505
## final  value 10.289505 
## stopped after 100 iterations

summary(model)

## Call:
## multinom(formula = Species ~ Petal.Length + Petal.Width, data = iris)
## 
## Coefficients:
##            (Intercept) Petal.Length Petal.Width
## versicolor   -22.79944      6.92122    7.878496
## virginica    -67.82521     12.64721   18.261016
## 
## Std. Errors:
##            (Intercept) Petal.Length Petal.Width
## versicolor     44.3859     37.58715    81.00888
## virginica      46.3939     37.65702    81.09482
## 
## Residual Deviance: 20.57901 
## AIC: 32.57901

13.4.2.2 Interpetasi Koefisien

z <- summary(model)$coefficients / summary(model)$standard.errors
p_values <- 2 * (1 - pnorm(abs(z)))
round(p_values, 4)

##            (Intercept) Petal.Length Petal.Width
## versicolor      0.6075       0.8539      0.9225
## virginica       0.1438       0.7370      0.8218

Hasil dari nilai p-value tidak menunjukkan bahwa adanya variabel prediktor yang berpengaruh signifikan terhadap log-odds dibandingkan dengan baseline (setosa)

13.4.2.3 Prediksi dan Visualisasi

iris$predicted <- predict(model, newdata = iris)
table(Predicted = iris$predicted, Actual = iris$Species)

##             Actual
## Predicted    setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         47         3
##   virginica       0          3        47

ggplot(iris, aes(x = Petal.Length, y = Petal.Width, color = predicted)) +
geom_point(size = 2) +
labs(title = "Multinomial Logistic Regression Predictions",
x = "Petal Length", y = "Petal Width") +
theme_minimal()

13.5 Kesimpulan

Model multinomial logistic regression efektif digunakna untuk klasifikasi kategori nominal lebih dari 2. Model ini memberikan estimasi log-odds terhadap baseline dan prediksi kategori berdaasarkan kombinasi prediktor

14 REGRESI LOGISTIK ORDINAL

Regresi logistik ordinal digunakan ketika variabel respon Y bersifat ordinal (memiliki urutan), misalnya tingkat kepuasan : rendah, sedang, tinggi

Perbedaan dengan model regresi logistik lainnya adalah :

• Regresi logistik biner : hanya 2 kategori

• Regresi logistik mulultinomial : kategori lebih dari 2 tetapi tidak berurutan.

14.1 Konsep Cumulative Logit Model

Model yang digunakan ada;aj cumulative logit model dengan asumsi proportional odds:

\[ \log \left( \frac{P(Y \leq j)}{P(Y > j)} \right) = \alpha_j + \beta x \]

dengan :

• \(\alpha_j\) : intercept khusus untuk kategori ke-j

• \(\beta\) : koefisien regresi (sama untuk semua kategori kumulatif)

Sama seperti regresi logistik multinomial, jika terdapat c kategori, maka terdapat c-1 model logit kumulatif

14.2 Interpretasi Koefisien

Koefisien \(\beta\) menjelaskan efek \(x\) terhadap kemungkinan berada pada kategori yang lebih rendah atau sama.

• Jika \(\beta > 0\) : semakin besar \(x\), semakin besar peluang berada di kategori rendah
• Jika \(\beta < 0\) : semakin besar \(x\), semakin besar peluang berada di kategori tinggi

Perhitungan odds ratio : \(OR = e^\beta\)

14.3 Contoh Kasus

Misal terdapat data simulasi tingkat kepuassan pelanggan (1 : Tidak Puas, 2: cukup, 3: puas) terhadap kecepatan layanan :

set.seed(123)
n <- 200 
speed <- round(runif(n, 1, 10))
satisfaction <- cut(5 + 0.5*speed + rnorm(n),
breaks = c(-Inf, 5.5, 7.5, Inf),
labels = c("Tidak Puas", "Cukup", "Puas"),
ordered_result = TRUE)
df <- data.frame(satisfaction, speed)
head(df)

##   satisfaction speed
## 1        Cukup     4
## 2         Puas     8
## 3        Cukup     5
## 4         Puas     9
## 5         Puas     9
## 6   Tidak Puas     1

14.3.1 Estimasi Model Ordinal

library(MASS)
model_ord <- polr(satisfaction~speed, data = df, Hess = TRUE)
summary (model_ord)

## Call:
## polr(formula = satisfaction ~ speed, data = df, Hess = TRUE)
## 
## Coefficients:
##        Value Std. Error t value
## speed 0.9096     0.1094   8.315
## 
## Intercepts:
##                  Value  Std. Error t value
## Tidak Puas|Cukup 1.3015 0.4377     2.9738 
## Cukup|Puas       4.4734 0.5718     7.8232 
## 
## Residual Deviance: 237.2312 
## AIC: 243.2312

Nantinya akan ada 2 model regresi logistik ordinal yang membedakan adalah interceptnya.

14.3.2 Nilai P-Value

ctable <- coef(summary(model_ord))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = round(p, 4)))

##                      Value Std. Error  t value p value
## speed            0.9095585  0.1093925 8.314630  0.0000
## Tidak Puas|Cukup 1.3015075  0.4376597 2.973789  0.0029
## Cukup|Puas       4.4733938  0.5718127 7.823180  0.0000

Dari hasil nilai p-value baik intercept maupun slope memiliki pengaruh signifikan terhadap log-odds kepuasan pelanggan

14.3.3 Prediksi Probabilitas

newdata <- data.frame(speed = 5:9)
predict(model_ord, newdata = newdata, type = "probs")

##    Tidak Puas     Cukup      Puas
## 1 0.037460604 0.4439482 0.5185912
## 2 0.015430723 0.2566765 0.7278928
## 3 0.006271788 0.1245723 0.8691559
## 4 0.002535158 0.0546231 0.9428417
## 5 0.001022461 0.0228089 0.9761686

14.4 Goodness-of-Fit dan Proportional odds

Model cumulative logit mengasumsikan efek prediktor sama untuk setiap cutoff. Jika tidak, pertimbangkan model non-proportional odds seperti generalized ordinal model.

14.5 Alternatif Model Ordinal

Selain cumulative logit, model untuk respon ordinal lainnya adalah :

• Adjacent-category logit

• Continuation-ratio (sequential) logit

Kedua model digunakan ketika asumsi proportional odds tidak terpenuhi.

14.6 Asumsi Paralelisme dalam Regresi Logistik Ordinal

Model regresi logistik ordinal yang paling umum adalah cumulative logit model yang harus memenuhi asumsi proportional odds. Asumsi proportional odds dikenal juga dengan asumsi paralelisme (parallel lines assumption)

Asumsi paralelisme menyatakan bahwa koefisien regresi \((\beta)\) sama untuk setiap kategori kumulatif dari variabel respon.

Dalam asumsi paralelisme, kurve logit kumulatif dari tiap kategori terhadap prediktor akan memiliki kemiringan yang sama (paralel), hanya berbeda posisi (intercept)

Jika asumsi ini tidak terpenuhi :

• Efek prediktor berbeda untuk setiap batas kategori

• Model cumulative logit tidak valid

• Perlu gunakan model alternatif (Generalized Ordinal Logistic Regression dan Partial Proporyional Odds Model)

Pengujian asumsi paralelisme dapat menggunakan Likelihood Ratio Test antara model proportional dan non-proportional atau Brant Test. Contoh pengujiannya adalah sebagai berikut

library(brant)
brant(model_ord)

## -------------------------------------------- 
## Test for X2  df  probability 
## -------------------------------------------- 
## Omnibus      1.13    1   0.29
## speed        1.13    1   0.29
## -------------------------------------------- 
## 
## H0: Parallel Regression Assumption holds

Hasil dari pengujian asumsi paralelisme dari simulasi data sebelumnya adalah model tersebut memenuhi asumsi paralelisme sehingga model cumulative logit tervalidasi.

14.7 Kesimpulan

• Regresi ordinal efektif untuk respon berurutan

• Model cumulative logit menginterpretasikan efek dalam bentuk log-odds kumulatif

• Asumsi paralelisme harus terpenuhi untuk model cumulative logit

• Jika asumsi paralelisme tidak terpenuhi harus menggunakan model alternatif.

15 LOG LINEAR MODEL

Ketika membangun model statistik yang dapat mengendalikan efek dari banyak variabel dan interaksinya secara simultan, model log-linear akan sangat berguna. Model log-linear adalah bentuk khusus dari Generalized Linear Model (GLM) yang digunakan pada frekuensi sel dalam tabel kontingensi dan mengasumsikan distribusi Poisson. Model log linear tidak seperti model regresi logistik karena model log-lienar tidak menetapkan variabel mana yang dependen dan independen, setiap variabel diperlakukan secara simetris. Model ini cocok untuk memahami struktur asosiasi atau independensi antar variabel dan bukan prediksi

Struktur model log-linear ditentukan berdasarkan efek utama dari masing-masing variabel serta interaksi di antara variabel-variabel tersebut. Misalnya, dalam tabel kontingensi tiga arah, model log-linear dapat membedakan apakah interaksi dua variabel cukup menjelaskan data atau diperlukan interaksi tiga arah untuk menangkap struktur asosianya. Penyesuaian model ini menggunakan metode likelihood ratio test untuk membandingkan model yang sederhana dengan yang lebih kompleks.

Sehingga dalam analisis data kategorik dapat dianalisis dengan beberapa pendekatan yaitu tabel kontingensi, model log-linear, dan model regresi logistik. Setiap pendekatan memiliki kelebihan dan kekuatan tergantung tujuan analisis dan struktur data. Tabel kontingensi bersifat deskriptif, model log-linear bersifat eksploratif terhadap hubungan simetris, dan regresi logistik bersifat prediktif terhadap outcome kategorik. Singkatnya, tabel kontingensi menyajikan frekuensi gabungan dari dua atau lebih variabel kategorik, model loglinear digunakan untuk memodelkan struktur asosiasi di dalam tabel kontingensi tanpa menganggap ada variabel dependen, dan model regresi logistik digunakan untuk memodelkan probabilitas dari kategori variabel dependen berdasarkan variabel independen.

15.1 Model Analisis Kategorik

Pada bab sebelumnya, tabel kontingensi sudah dipelajari yaitu tabel yang menyajikan jumlah frekuensi dari kombinasi kategori antar variabel. contoh tabel 2x2 adalah sebagai berikut

table_data <- matrix(c(30, 20, 50, 70), nrow=2,
dimnames = list(Obat = c("Timolol", "Placebo"),
Serangan = c("Ya", "Tidak")))
table_data

##          Serangan
## Obat      Ya Tidak
##   Timolol 30    50
##   Placebo 20    70

dengan catatan tabel kontingensi tersebut bersifat deskriptif dan tidak melibatkan pemodelan probabilitas

Model loglinear memodelkan logaritma dari ekspektasi frekuensi sel dalam tabel kontingensi yang ada. Bentuk umumnya adalah sebagai berikut

\[ \log(\mu_{ij}) = \mu + \lambda_i^A +\lambda_j^B + \lambda_{ij}^{AB} \]

Model loglinear ini :

• cocok untuk menyelidiki asosiasi dan independensi antar variabel,

• tidak membedakan variabel respon dan penjelas,

• umumnya digunakan pada tabel muti-dimensi (2 arah, 3 arah, dst)

Contoh model loglinear dengan tabel kontingensi contoh sebelumnya

library(MASS)
loglm(~Obat*Serangan, data = table_data)

## Call:
## loglm(formula = ~Obat * Serangan, data = table_data)
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1

Selain itu analisis kategorik dapat menggunakan model regresi logistik biner, dengan model sebagai berikut :

\[ \log \left(\frac{p}{1-p}\right) = \beta_0 + \beta_1x \]

model regresi logistik biner :

• digunakan jika ada variabel dependen kategorik (biasanya bersifat biner)

• bertujuan untuk memprediksi probabilitas

• umumnya digunakan dalam studi observasional atau eksperimental

Contoh penggunaan regresi logistik biner adalah sebagai berikut

data_glm <- data.frame(
Serangan = c(1, 0, 1, 0),
Obat = factor(c("Timolol", "Timolol", "Placebo", "Placebo")),
Frek = c(30, 20, 50, 70)
)
model_logit <- glm(Serangan ~ Obat, weights = Frek, family = binomial, data = data_glm)
summary(model_logit)

## 
## Call:
## glm(formula = Serangan ~ Obat, family = binomial, data = data_glm, 
##     weights = Frek)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.3365     0.1852  -1.817   0.0692 .
## ObatTimolol   0.7419     0.3430   2.163   0.0305 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 235.08  on 3  degrees of freedom
## Residual deviance: 230.31  on 2  degrees of freedom
## AIC: 234.31
## 
## Number of Fisher Scoring iterations: 4

15.2 Model Saturated

Model saturated atau model penuh menyertakan seluruh efek utama dan interaksi. Model saturated contohnya adalah sebagai berikut (menggunakan tabel kontingensi 2x2) :

library(MASS)
data <- matrix(c(35, 65, 45, 55), nrow=2, byrow=TRUE)
dimnames(data) <- list(Obat = c("Timolol", "Placebo"), Serangan = c("Ya", "Tidak"))
ftable(data)

##         Serangan Ya Tidak
## Obat                     
## Timolol          35    65
## Placebo          45    55

model_saturated <- loglm(~ Obat * Serangan, data = data)
summary(model_saturated)

## Formula:
## ~Obat * Serangan
## attr(,"variables")
## list(Obat, Serangan)
## attr(,"factors")
##          Obat Serangan Obat:Serangan
## Obat        1        0             1
## Serangan    0        1             1
## attr(,"term.labels")
## [1] "Obat"          "Serangan"      "Obat:Serangan"
## attr(,"order")
## [1] 1 1 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1

Model saturatednya adalah sebagai berikut :

\[ \log(\mu_{ij}) = \lambda + \lambda_{i}^{Obat} + \lambda_{j}^{Serangan} + \lambda_{ij}^{Obat,Serangan} \]

dan hasil dari pengujian signifikansi model dengan likelihood ratio maupun pearson menunjukkan bahwa model tersebut cocok untuk menjelaskan semua variasi pada data.

15.3 Model Independent

Tidak seperti model saturated, model independen mengasumsikan bahwa tidak ada interaksi antara variabelnya, sehingga model umumnya adalah sebagai berikut :

\[ \log(\mu_{ij}) = \lambda + \lambda_{i}^{A} + \lambda_{j}^{B} \]

Model ini menguji hipotesis bahwa variabel X dan Y saling independen

model_indep <- loglm(~ Obat + Serangan, data = data)
summary(model_indep)

## Formula:
## ~Obat + Serangan
## attr(,"variables")
## list(Obat, Serangan)
## attr(,"factors")
##          Obat Serangan
## Obat        1        0
## Serangan    0        1
## attr(,"term.labels")
## [1] "Obat"     "Serangan"
## attr(,"order")
## [1] 1 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
## 
## Statistics:
##                       X^2 df  P(> X^2)
## Likelihood Ratio 2.087576  1 0.1485015
## Pearson          2.083333  1 0.1489147

dengan \(H_0\) bahwa variabel bersifat saling independen, dengan p-value yang dihasilkan dari pengujian dengan likelihood ratio maupun pearson menunjukkan bahwa variabel variabel bersifat independen sehingga model independen ini cocok digunakan.

15.4 Odds Ratio dan Interpretasi

Odds ratio untuk tabel 2x2 :

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

Interpretasi nilai Odds Ratio :

• OR = 1 : tidak ada asosiasi

• OR > 1 : asosiasi positif

• OR < 1 : asosiasi negatiff

15.5 Estimasi Parameter

Dalam model saturated :

• estimasi dilakukan dengan pembatasan seperti sum-to-zero

• estimasi parameter dilakukan dengan iterative proportional fitting (IPF)

logOR <- log((data[1,1] * data[2,2]) / (data[1,2] * data[2,1]))
logOR

## [1] -0.4183685

Dari hasil perhitungan nilai odds rationya, variabel obat dan serangan memiliki asosiasi negatif.

15.6 Perbandingan Model

Dalam pembuatan model, pasti harus dipilih model yang terbaik. Perbandingan antar model dilakukan dengan mneggunakan statistik deviance \((G^2)\) atau likelihood ratio test. Contohnya adalah sebagai berikut

anova(model_indep, model_saturated)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Obat + Serangan 
## Model 2:
##  ~Obat * Serangan 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   2.087576  1                                    
## Model 2   0.000000  0   2.087576         1         0.1485
## Saturated 0.000000  0   0.000000         0         1.0000

dari hasil pengujian didapatkan nilai p-value lebih besar dari 0.05 yang menunjukkan bahwa model tanpa interaksi sudah cukup baik untuk menjelaskan data. Karena model yang diinginkan adalah model yang parsimony, maka model yang dipilih adalah model independen

15.7 Contoh Kasus

Studi dari Agresti (2019) membahas hubungan antara kebahagiaan dan kepercayaan terhadap kehidupan akhirat. Datanya adalah sebagai berikut :

data_survey <- matrix(c(32,190,
113,611,
51,326),
nrow = 3, byrow = TRUE,
dimnames = list(Kebahagiaan = c("Tidak", "Cukup", "Sangat"),
Surga = c("Tidak Percaya", "Percaya")))
ftable(data_survey)

##             Surga Tidak Percaya Percaya
## Kebahagiaan                            
## Tidak                        32     190
## Cukup                       113     611
## Sangat                       51     326

modelsaturated <- loglm(~ Kebahagiaan*Surga, data = data_survey)
modelindep <- loglm(~ Kebahagiaan + Surga, data = data_survey)
anova (modelsaturated, modelindep)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Kebahagiaan + Surga 
## Model 2:
##  ~Kebahagiaan * Surga 
## 
##            Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   0.8911136  2                                    
## Model 2   0.0000000  0  0.8911136         2        0.64047
## Saturated 0.0000000  0  0.0000000         0        1.00000

Dari hasil pengujian modle saturated dan model independen didapatkan bahwa model independen sudah cukup untuk menjelaskan variabel. maka model log linearnya adalah sebagai berikut :

\[ \log(\mu_{ij}) = \lambda + \lambda_{i}^{Kebahagiaan} + \lambda_{j}^{Kepercayaan} \]

16 MODEL LOG LINEAR DUA ARAH

16.1 Model log-linear pada tabel kontingensi

model log linear adalah model yang digunakan untuk menganalisis hubungan antara dua atau lebih variabel kategorik yang disajikan dalam tabel kontingensi. Model ini mengasumsikan bahwa logaritma dari nilai ekspektasi frekuensi sel \[ (\mu_{ij})\ \] dapat dinyatakan sebagai penjumlahan efek variabel dan (bila perlu) interaksinya. Model log-linear untuk tabel kontingensi 2x2 adalah :

\[ \log(\mu_{ij}) = \lambda + \lambda_i^A + \lambda_j^B + \lambda_{ij}^{AB} \]

16.2 Perbedaan Model Log-Linear dan Model Regresi Logistik

• Model log linear digunakan untuk memeodelkan frekuensi pada tabel kontingensi dan menguji asosiasi antar variabel kategorik, tanpa menganggap ada variabel respon dan prediktor

• Model regresi logistik digunakan untuk memodelkan probabilitas kejadian suatu outcome (biner) berdasarkan satu atau lebih prediktor (bisa kategorik atau kontinu)

16.3 Estimasi Parameter Log Linear 2 Arah

Sistem peramaan model log linear :

\[ \begin{aligned} \log(\mu_{11}) &= \lambda + \lambda_1^A + \lambda_1^B + \lambda_{11}^{AB} \\ \log(\mu_{12}) &= \lambda + \lambda_1^A + \lambda_2^B + \lambda_{12}^{AB} \\ \log(\mu_{21}) &= \lambda + \lambda_2^A + \lambda_1^B + \lambda_{21}^{AB} \\ \log(\mu_{22}) &= \lambda + \lambda_2^A + \lambda_2^B + \lambda_{22}^{AB} \end{aligned} \]

Constraint Sum-to-Zero

\[ \begin{aligned} \lambda_1^A + \lambda_2^A &= 0 \\ \lambda_1^B + \lambda_2^B &= 0 \\ \lambda_{11}^{AB} + \lambda_{12}^{AB} + \lambda_{21}^{AB} + \lambda_{22}^{AB} &= 0 \end{aligned} \]

Rumus Estimasi Parameter dengan Sum-to-Zero Constraint

\[ \lambda_1^A = \frac{1}{2} \left[ (\log \mu_{11} + \log \mu_{12}) - (\log \mu_{21} + \log \mu_{22}) \right] \]

\[ \lambda_1^B = \frac{1}{2} \left[ (\log \mu_{11} + \log \mu_{21}) - (\log \mu_{12} + \log \mu_{22}) \right] \]

\[ \lambda_{12}^{AB} = \frac{1}{4} \left[ \log \mu_{12} - \log \mu_{11} - \log \mu_{22} + \log \mu_{21} \right] \]

16.4 Contoh

16.4.1 Analisis Data Tabel Kontingensi 2x2

Diberikan Data

Merokok	Sakit	Sehat
Ya	30	20
Tidak	10	40

16.4.2 Model Log-Linear pada Tabel 2x2

\[ \log(\mu_{ij}) = \lambda + \lambda_i^A + \lambda_j^B + \lambda_{ij}^{AB} \]

dengan constraint sum-to-zero :

\[ \sum_i \lambda_i^A = 0,\quad \sum_j \lambda_j^B = 0,\quad \sum_{i,j} \lambda_{ij}^{AB} = 0 \]

16.4.3 Estimasi Parameter Model

Misalkan :

• A1 = Merokok (Ya), A2 = Tidak

• B1 = Sakit, B2 = Sehat

Observasi :

• \(n_{11} = 30, n_{12} = 20\)

• $n_{21} = 10, n_{22} = 40 $

Model yang akan terbentuk

\[ \log(\mu_{11}) = \lambda + \lambda_1^A + \lambda_1^B + \lambda_{11}^{AB} \]\[ \log(\mu_{12}) = \lambda + \lambda_1^A + \lambda_2^B + \lambda_{12}^{AB} \]\[ \log(\mu_{21}) = \lambda + \lambda_2^A + \lambda_1^B + \lambda_{21}^{AB} \]\[ \log(\mu_{22}) = \lambda + \lambda_2^A + \lambda_2^B + \lambda_{22}^{AB} \]

dengan constraint sum-to-zero :

\(\lambda_1^A + \lambda_2^A =0\)

\(\lambda_1^B + \lambda_2^B =0\)

\(\lambda_{11}^{AB} + \lambda_{12}^{AB} + \lambda_{21}^{AB} + \lambda_{22}^{AB} =0\)

Perhitungan :

1. Rata-rata log frekuensi sel

\[ \frac{1}{4} \sum_{i=1}^2 \sum_{j=1}^2 \log(n_{ij}) = \frac{1}{4} \left[\log(30) + \log(20) + \log(10) + \log(40)\right] \]

{r} n <- matrix(c(30,20,10,40), nrow = 2, byrow = TRUE) mean_log_n <- mean(log(n)) ; mean_log_n}

\[ \lambda = 3.0971 \]

2. Efek utama A (Merokok)

\[ \lambda_1^A = \frac{1}{2}[(\log(30) + \log(20))] - [\log(10) + \log(40))] \]

n <- matrix(c(30,20,10,40),nrow = 2, byrow = TRUE)
log_n <- log(n)
str(log_n)

##  num [1:2, 1:2] 3.4 2.3 3 3.69

alpha_1 <- 0.5 * (sum(log_n[1, ]) - sum(log_n[2, ]))
alpha_2 <- -alpha_1

alpha_1

## [1] 0.2027326

alpha_2

## [1] -0.2027326

\[ \lambda_1^A = 0.2027 \]\[ \lambda_2^A = -0.2027 \]

3. Efek utama B (Status)

\[ \lambda_1^B = \frac{1}{2}[(\log(30) + \log(10))] - [\log(20) + \log(40))] \]

alpha_1 <- 0.5 * (sum(log_n[,1]) - sum(log_n[,2])) 
alpha_2 <- -alpha_1  
alpha_1 

## [1] -0.4904146

alpha_2

## [1] 0.4904146

\[ \lambda_1^B = -0.4904 \]\[ \lambda_2^B = 0.4904 \]

4. Efek interaksi AB

\[ \lambda_{11}^{AB} = \frac{1}{4}[\log(30) - \log(10) - \log(20) + \log(40)] \]

alpha_11 <- 0.25*(log(30)-log(10)-log(20)+log(40)) 
alpha_11

## [1] 0.4479399

\[ \lambda_{11}^{AB} = 0.4479 \]\[ \lambda_{12}^{AB} = -\lambda_{11}^{AB} = -0.4479 \]\[ \lambda_{21}^{AB} = 0.4479 \]\[ \lambda_{12}^{AB} = -0.4479 \]

16.4.4 Hitung Odds Ratio dan Interval Kepercayaan

\[OR = \frac{n_{11}n_{22}}{n_{12}n_{21}} = \frac{30\times40}{20\times10}\]

or <- (30*40)/(20*10) ; or  

## [1] 6

log_or <- log(or) ; log_or

## [1] 1.791759

\[ OR = 6 \]\[ \log(OR) = 1.7198 \]

Standard Error (SE)

\[ SE = \sqrt{\frac{1}{n_{11}}+\frac{1}{n_{12}}+\frac{1}{n_{21}}+\frac{1}{n_{22}}} = 0.4564 \]

se <- sqrt(1/30 + 1/20 + 1/10 + 1/40) ; se

## [1] 0.4564355

Confidence Interval untuk log(OR) :

\[ \log(OR) \pm 1.96 \times SE \]

lower <- log_or - 1.96*se
upper <- log_or + 1.96*se

lower_backtf <- exp(lower)
upper_backtf <- exp(upper)
cat(lower_backtf,upper_backtf)

## 2.452593 14.67834

\[ CI : 2.45 - 14.68 \]

16.4.5 Fitting Model Log-Linear

# Data 2x2
tabel <- matrix(c(30, 20, 10, 40), nrow = 2, byrow = TRUE)
colnames(tabel) <- c("Sakit", "Sehat")
rownames(tabel) <- c("Ya", "Tidak")
tabel

##       Sakit Sehat
## Ya       30    20
## Tidak    10    40

data <- as.data.frame(as.table(tabel))
colnames(data) <- c("Merokok", "Status", "Freq")
data

##   Merokok Status Freq
## 1      Ya  Sakit   30
## 2   Tidak  Sakit   10
## 3      Ya  Sehat   20
## 4   Tidak  Sehat   40

Model tanpa interaksi (independen)

fit_no_inter <- loglm(Freq ~ Merokok + Status, family = poisson, data = data) 
summary(fit_no_inter)

## Formula:
## Freq ~ Merokok + Status
## attr(,"variables")
## list(Freq, Merokok, Status)
## attr(,"factors")
##         Merokok Status
## Freq          0      0
## Merokok       1      0
## Status        0      1
## attr(,"term.labels")
## [1] "Merokok" "Status" 
## attr(,"order")
## [1] 1 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(Freq, Merokok, Status)
## attr(,"dataClasses")
##      Freq   Merokok    Status 
## "numeric"  "factor"  "factor" 
## 
## Statistics:
##                       X^2 df     P(> X^2)
## Likelihood Ratio 17.26092  1 3.258188e-05
## Pearson          16.66667  1 4.455709e-05

Model dengan interaksi (saturated)

fit_inter <- loglm(Freq ~ Merokok * Status, family = poisson, data = data) 
summary(fit_inter)

## Formula:
## Freq ~ Merokok * Status
## attr(,"variables")
## list(Freq, Merokok, Status)
## attr(,"factors")
##         Merokok Status Merokok:Status
## Freq          0      0              0
## Merokok       1      0              1
## Status        0      1              1
## attr(,"term.labels")
## [1] "Merokok"        "Status"         "Merokok:Status"
## attr(,"order")
## [1] 1 1 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(Freq, Merokok, Status)
## attr(,"dataClasses")
##      Freq   Merokok    Status 
## "numeric"  "factor"  "factor" 
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1

Dari hasil pengujian model independen dan model saturated (interaksi) model menunjukkan adanya interaksi karena variabel tidak independen sehingga model yang digunakan adalah model saturated.

16.4.6 Interpretasi Parameter

• intercept atau parameter utama menunjukkan rata-rata log-frekuensi sel.

• Efek “merokok” dan “status” menunjukkan perbedaan log frekuensi antar kategori

• interaksi signifikan menunjukkan adanya asosisasi antara merokok dan status.

16.5 Model Log Linear untuk Tabel 2x3

Bentuk umum model og-linear untuk tabel 2x3 (dengan sum-to-zero constraint) :

\[ \log(\mu_{ij}) = \lambda + \lambda_i^A + \lambda_j^B +\lambda_{ij}^{AB} \]

dengan :

• \(\lambda_{ij}\) : ekspektasi frekuensi pada bari ske-i dan komlom ke-j

• A : kategori dalam baris

• B : kategori dalam kolom

Misalkan :

Suatu survei dilakukan untuk mengetahui hubungan antara jenis kelamin (perempuan/laki-laki) dan kategori BMI (kurus/normal/gemuk) :

\[ \begin{array}{c|ccc} \text{} & \text{Kurus} & \text{Normal} & \text{Gemuk} \\ \hline \text{Laki-laki} & 12 & 20 & 8 \\ \text{Perempuan} & 18 & 24 & 10 \\ \end{array} \]

Maka bentuk model loglinearnya adalah sebagai berikut :

\[ \log(\mu_{ij}) = \lambda + \lambda_1^A + \lambda_j^B +\lambda_{ij}^{AB} \]

dengan

A : jenis kelamin (i = 1 : Laki-laki, i = 2 : Perempuan)

B : Kategori BMI :j = 1 : Kurus, j = 2 : Normal, j = 3 : Gemuk)

16.6 Fitting Model Log-Linear di R

#Tabel Kotingensi
tabel2x3 <- matrix(c(12, 20, 8, 18, 24, 10), nrow = 2, byrow = TRUE)
colnames(tabel2x3) <- c("Kurus", "Normal", "Gemuk")
rownames(tabel2x3) <- c("Laki-laki", "Perempuan")
tabel2x3

##           Kurus Normal Gemuk
## Laki-laki    12     20     8
## Perempuan    18     24    10

#Tabel untuk glm 
data2x3 <- as.data.frame(as.table(tabel2x3))
colnames(data2x3) <- c("JenisKelamin", "BMI", "Freq")
data2x3

##   JenisKelamin    BMI Freq
## 1    Laki-laki  Kurus   12
## 2    Perempuan  Kurus   18
## 3    Laki-laki Normal   20
## 4    Perempuan Normal   24
## 5    Laki-laki  Gemuk    8
## 6    Perempuan  Gemuk   10

#Model log-linear tanpa interaksi (independen)
fit_no_inter <- glm(Freq ~ JenisKelamin + BMI, family = poisson, data = data2x3)
summary(fit_no_inter)

## 
## Call:
## glm(formula = Freq ~ JenisKelamin + BMI, family = poisson, data = data2x3)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             2.5683     0.2179  11.789   <2e-16 ***
## JenisKelaminPerempuan   0.2624     0.2103   1.248   0.2122    
## BMINormal               0.3830     0.2368   1.618   0.1058    
## BMIGemuk               -0.5108     0.2981  -1.713   0.0866 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 13.06443  on 5  degrees of freedom
## Residual deviance:  0.22527  on 2  degrees of freedom
## AIC: 35.26
## 
## Number of Fisher Scoring iterations: 3

#Model loglinear dengan interaksi (saturated model)
fit_inter <- glm(Freq ~ JenisKelamin * BMI, family = poisson, data = data2x3)
summary(fit_inter)

## 
## Call:
## glm(formula = Freq ~ JenisKelamin * BMI, family = poisson, data = data2x3)
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                       2.4849     0.2887   8.608   <2e-16 ***
## JenisKelaminPerempuan             0.4055     0.3727   1.088    0.277    
## BMINormal                         0.5108     0.3651   1.399    0.162    
## BMIGemuk                         -0.4055     0.4564  -0.888    0.374    
## JenisKelaminPerempuan:BMINormal  -0.2231     0.4802  -0.465    0.642    
## JenisKelaminPerempuan:BMIGemuk   -0.1823     0.6032  -0.302    0.762    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  1.3064e+01  on 5  degrees of freedom
## Residual deviance: -9.0719e-30  on 0  degrees of freedom
## AIC: 39.034
## 
## Number of Fisher Scoring iterations: 3

#Membandingkan model 
anova(fit_no_inter, fit_inter)

## Analysis of Deviance Table
## 
## Model 1: Freq ~ JenisKelamin + BMI
## Model 2: Freq ~ JenisKelamin * BMI
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1         2    0.22526                     
## 2         0    0.00000  2  0.22526   0.8935

Interpretasi :

• Dari hasil pengujian model yang lebih baik adalah model tanpa interaksi karena interaksi yang ada pada model tidak signifikan ditunjukkan pada analysis of deviance tabel

• Jenis Kelamin dan BMI independen

• Model tanpa interkasi menggunakan jenis kelaimin laki-laki dan bmi kurus sebagai base-line.

17 MODEL LOG LINEAR TIGA ARAH

Penambahan dari model loglinear dua arah, model loglinear tiga arah memungkinkan adanya interaksi yang lebih banyak dibandingkan yang dua arah. Model ini melibatkan tuga variabel yang mungkin diinteraksikan secara bersamaan

17.1 Model Log-Linear untuk Tabel Tiga Arah

Dalam tabel tiga arah yang melibatkan tiga variabel kategorik (X,Y dan X) dapat dibentuk berbagai bentuk model tergantung interaksi yang ingin dimasukkan. Bentuknya adalah sebagai berikut :

1. Model Saturated

\[ \log(\mu_{ijk}) = \lambda + \lambda_{i}^{X} + \lambda_{j}^{Y} + \lambda_{k}^{Z} + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} + \lambda_{ijk}^{XYZ} \]

Model ini memuat semua kemungkinana interaksi, termasuk interaksi tiga arahnya

2. Model Homogen

\[ \log(\mu_{ijk}) = \lambda + \lambda_{i}^{X} + \lambda_{j}^{Y} + \lambda_{k}^{Z} + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \]

Model ini memuat interaksi dua arahnya tetapi tidak interaksi tiga arahnya

3. Model Conditional

Model conditional hanya memasukkan interakasi dari variabel payang diinginkan

• Conditional pada X :

\[ \log(\mu_{ijk}) = \lambda + \lambda_{i}^{X} + \lambda_{j}^{Y} + \lambda_{k}^{Z} + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} \]

• Conditional pada Y

\[ \log(\mu_{ijk}) = \lambda + \lambda_{i}^{X} + \lambda_{j}^{Y} + \lambda_{k}^{Z} + \lambda_{ij}^{XY} + \lambda_{jk}^{YZ} \]

• Conditional pada Z

  \[ \log(\mu_{ijk}) = \lambda + \lambda_{i}^{X} + \lambda_{j}^{Y} + \lambda_{k}^{Z} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \]

4. Model Joint Independence

Model ini hanya memasukkan interkasi variabel yang independensi

• Independensi antara X dan Y

\[ \log(\mu_{ijk}) = \lambda + \lambda_{i}^{X} + \lambda_{j}^{Y} + \lambda_{k}^{Z} + \lambda_{ij}^{XY} \]

• Independensi antara X dan Z

\[ \log(\mu_{ijk}) = \lambda + \lambda_{i}^{X} + \lambda_{j}^{Y} + \lambda_{k}^{Z} + \lambda_{ik}^{XZ} \]

• Independensi antara Y dan Z

\[ \log(\mu_{ijk}) = \lambda + \lambda_{i}^{X} + \lambda_{j}^{Y} + \lambda_{k}^{Z} + \lambda_{jk}^{YZ} \]

5. Model Tanpa Interaksi

Model ini hanya memasukkan efek utama tanpa interaksi antar variabelnya

\[ \log(\mu_{ijk}) = \lambda + \lambda_{i}^{X} + \lambda_{j}^{Y} + \lambda_{k}^{Z} \]

17.2 Pengujian Interaksi dalam Model Log-Linear Tiga Arah

Pengujian ini dilakukan untuk mengetahui apakah ada interaksi antar variabel atau tidak. Pengujiann ini tidak dapat dilakukan langsung tetapi harus bertahap dari yang interaksi tertinggi ke yang lebih rendah. Tahapannya adalah sebagai berikut :

1. Pengujian interaksi tiga arah
   • Bandingkan model saturated dengan model homogenous
2. Pengujian Interaksi dua arah

   • Bandingkan model homogenous dengan model conditional

   • Bandingkan model conditiosnal dengan model joint independence

   • Bandingkan model joint independence dengan model tanpa interakasi

Setiap pengujian dilakukan untuk menilai kecocokan model dan menentukan struktur interaksi mana yang paling sesuai dengan data

17.3 Studi Kasus

Tabel berikut menyajikan data dari survei General Sosial Survey (GSS) tahun 1994 mengenai jenis kelamin responden, tingkat fundamentalisme, dan sikap terhadap hukuman mati untuk kasus pembunuhan. Susun dan interpretasikan model log-linear paling sederhana (paling parsimonious) untuk data ini. Jelaskan proses yang anda lakukan dalam menentukan model terbaik serta asosiasi apa saja yang teridentifikasi. Tunjukkan juga bagaimana nilai yang diprediksi dari model menggambarkan asosiasi tersebut.

\[ \begin{array}{c|c|c|c|c} \text{Fundamentalisme} & \text{Jenis Kelamin} & \text{Mendukung} & \text{Menolak} & \text{Total} \\ \hline \text{Fundamentalist} & \text{Laki-laki} & 128 & 32 & 160 \\ \text{Fundamentalist} & \text{Perempuan} & 123 & 73 & 196 \\ \text{Fundamentalist} & \text{Total} & 251 & 105 & 356 \\ \hline \text{Moderate} & \text{Laki-laki} & 182 & 56 & 238 \\ \text{Moderate} & \text{Perempuan} & 168 & 105 & 273 \\ \text{Moderate} & \text{Total} & 350 & 161 & 511 \\ \hline \text{Liberal} & \text{Laki-laki} & 119 & 49 & 168 \\ \text{Liberal} & \text{Perempuan} & 111 & 70 & 181 \\ \text{Liberal} & \text{Total} & 230 & 119 & 349 \\ \end{array} \]

17.3.1 Model Saturated

Input Data :

# Input data sesuai tabel praktikum
z.fund <- factor(rep(c("1fund", "2mod", "3lib"), each = 4))
x.sex <- factor(rep(c("1M", "2F"), each = 2, times = 3))
y.fav <- factor(rep(c("1fav", "2opp"), times = 6))
counts <- c(128, 32, 123, 73, 182, 56, 168, 105, 119, 49, 111, 70)
data <- data.frame(
Fundamentalisme = z.fund,
Jenis_Kelamin = x.sex,
Sikap = y.fav,
Frekuensi = counts
)
data

##    Fundamentalisme Jenis_Kelamin Sikap Frekuensi
## 1            1fund            1M  1fav       128
## 2            1fund            1M  2opp        32
## 3            1fund            2F  1fav       123
## 4            1fund            2F  2opp        73
## 5             2mod            1M  1fav       182
## 6             2mod            1M  2opp        56
## 7             2mod            2F  1fav       168
## 8             2mod            2F  2opp       105
## 9             3lib            1M  1fav       119
## 10            3lib            1M  2opp        49
## 11            3lib            2F  1fav       111
## 12            3lib            2F  2opp        70

#Membentuk Tabel Kontingensi 3 Arah 
table3d <- xtabs(Frekuensi ~ Fundamentalisme + Jenis_Kelamin + Sikap,data = data)
ftable(table3d)

##                               Sikap 1fav 2opp
## Fundamentalisme Jenis_Kelamin                
## 1fund           1M                   128   32
##                 2F                   123   73
## 2mod            1M                   182   56
##                 2F                   168  105
## 3lib            1M                   119   49
##                 2F                   111   70

17.3.1.1 Pembentukan Model

Model saturated adalh model log-linear saturated memasukkan semua interaksi ketiga variabel yang terlibat.

x.sex <- relevel(x.sex, ref = "2F")
y.fav <- relevel(y.fav, ref = "2opp")
z.fund <- relevel(z.fund, ref = "3lib")
model_saturated <- glm(counts ~ x.sex + y.fav + z.fund +
                         x.sex*y.fav + x.sex*z.fund + y.fav*z.fund +
                         x.sex*y.fav*z.fund,
                       family = poisson(link = "log"))
summary(model_saturated)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     x.sex * z.fund + y.fav * z.fund + x.sex * y.fav * z.fund, 
##     family = poisson(link = "log"))
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    4.248495   0.119523  35.545  < 2e-16 ***
## x.sex1M                       -0.356675   0.186263  -1.915  0.05551 .  
## y.fav1fav                      0.461035   0.152626   3.021  0.00252 ** 
## z.fund1fund                    0.041964   0.167285   0.251  0.80193    
## z.fund2mod                     0.405465   0.154303   2.628  0.00860 ** 
## x.sex1M:y.fav1fav              0.426268   0.228268   1.867  0.06185 .  
## x.sex1M:z.fund1fund           -0.468049   0.282210  -1.659  0.09721 .  
## x.sex1M:z.fund2mod            -0.271934   0.249148  -1.091  0.27507    
## y.fav1fav:z.fund1fund          0.060690   0.212423   0.286  0.77511    
## y.fav1fav:z.fund2mod           0.008969   0.196903   0.046  0.96367    
## x.sex1M:y.fav1fav:z.fund1fund  0.438301   0.336151   1.304  0.19227    
## x.sex1M:y.fav1fav:z.fund2mod   0.282383   0.301553   0.936  0.34905    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.4536e+02  on 11  degrees of freedom
## Residual deviance: 5.9952e-15  on  0  degrees of freedom
## AIC: 100.14
## 
## Number of Fisher Scoring iterations: 3

exp(model_saturated$coefficients)

##                   (Intercept)                       x.sex1M 
##                    70.0000000                     0.7000000 
##                     y.fav1fav                   z.fund1fund 
##                     1.5857143                     1.0428571 
##                    z.fund2mod             x.sex1M:y.fav1fav 
##                     1.5000000                     1.5315315 
##           x.sex1M:z.fund1fund            x.sex1M:z.fund2mod 
##                     0.6262231                     0.7619048 
##         y.fav1fav:z.fund1fund          y.fav1fav:z.fund2mod 
##                     1.0625694                     1.0090090 
## x.sex1M:y.fav1fav:z.fund1fund  x.sex1M:y.fav1fav:z.fund2mod 
##                     1.5500717                     1.3262868

Model ini memodelkan hubungan antara jenis kelamin, sikap terhadap hukuman, dan tingkat fundamentalisme terhadap frekuensi responden. Model ini memperhatikan interakksi antar dua faktor dan tiga faktor dengan jenis kelamin perempuan, sikap menolak, dan fundamentalisme liberal sebagai baseline.

17.3.1.2 Interpretasi Koefisien

• Intercept : rata-rata log jumlah kasus untuk baseline adalah 4.25

• x.sex1M : Laki-laki memiliki expected count 0.7 kali dari perempuan, tetapi koefisien tidak signifikan

• y.fav1fav : orang yang mendukung hukuman mati memiliki expected count 1.59 kali lipat lebih tinggi dibandingkan yang menolak dan hal ini signifikan

• z.fund1fund : kelompok fundamentalist tidak berbeda signifikan dari yang liberal

• z.fund2mod : kelompok moderate memiliki expected count 1.5 kali lebih besar dibandingkan liberal dan hal ini signifikan

• Untuk interaksi dua arah dan tiga arah tidak signifikan sehingga tidak ada bukti kuat adaanya efek interaksi antar variabel

17.3.1.3 Goodness-of-Fit

• Residual deviance : nilainya sangat kecil dan mendekati 0 sehingga model ini fit terhadap data

• AIC = 100.14 yang nantinya dapat digunakan untuk membandingkan model dan memilih model yang lebih sederhana

17.3.1.4 Kesimpulan

• Model saturated sangat fit dengan data, tetapi masih ada beberapa parameter dan interaksi yang signifikan

• Efek yang paling signifikan :

  • Sikap mendukung hukuman mati

  • Kelompok Moderate

• Model masih perlu dipertimbangkan dengan model yang lebih sederhana mengingat interaksinya tidak signifikan

17.3.2 Model Homogenous

Model log-linear homogenous memasukkan semua efek utama dan interaksi dua arahnya.

17.3.2.1 Pembentukan Model

model_homogenous <- glm(counts ~ x.sex + y.fav + z.fund +
                          x.sex*y.fav + x.sex*z.fund + y.fav*z.fund,
                        family = poisson(link = "log"))
summary(model_homogenous)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     x.sex * z.fund + y.fav * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            4.31096    0.10522  40.972  < 2e-16 ***
## x.sex1M               -0.51575    0.13814  -3.733 0.000189 ***
## y.fav1fav              0.35707    0.12658   2.821 0.004788 ** 
## z.fund1fund           -0.06762    0.14452  -0.468 0.639854    
## z.fund2mod             0.33196    0.13142   2.526 0.011540 *  
## x.sex1M:y.fav1fav      0.66406    0.12728   5.217 1.81e-07 ***
## x.sex1M:z.fund1fund   -0.16201    0.15300  -1.059 0.289649    
## x.sex1M:z.fund2mod    -0.08146    0.14079  -0.579 0.562887    
## y.fav1fav:z.fund1fund  0.23873    0.16402   1.455 0.145551    
## y.fav1fav:z.fund2mod   0.13081    0.14951   0.875 0.381614    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.361  on 11  degrees of freedom
## Residual deviance:   1.798  on  2  degrees of freedom
## AIC: 97.934
## 
## Number of Fisher Scoring iterations: 3

exp(model_homogenous$coefficients)

##           (Intercept)               x.sex1M             y.fav1fav 
##            74.5122374             0.5970531             1.4291312 
##           z.fund1fund            z.fund2mod     x.sex1M:y.fav1fav 
##             0.9346132             1.3936992             1.9426624 
##   x.sex1M:z.fund1fund    x.sex1M:z.fund2mod y.fav1fav:z.fund1fund 
##             0.8504296             0.9217741             1.2696313 
##  y.fav1fav:z.fund2mod 
##             1.1397492

17.3.2.2 Membandingkan Model Saturated dengan Model Homogenous

Pengujian ini menggunakan residual deviance dari kedua model

Hipotesis :

\(H_0\) : Tidak ada interaksi tiga arah (model homogenous sudah cukup)

\(H_1\) : Ada interaksi tiga arah (model saturated dibutuhkan)

Statistik uji : selisih deviance

deviance <- model_homogenous$deviance - model_saturated$deviance
deviance

## [1] 1.797977

Daerah kritis :

df <- model_homogenous$df.residual - model_saturated$df.residual
df

## [1] 2

chi.tabel <- qchisq(1 - 0.05, df = df)
chi.tabel

## [1] 5.991465

Keputusan

keputusan <- ifelse(deviance <= chi.tabel, "Terima H0", "Tolak H0")
keputusan

## [1] "Terima H0"

Kesimpulan :

Dengan taraf signifikansi 5%, tidak ada interaksi tiga arah sehingga model homogenous sudah cukup.

17.3.3 Model Conditional pada X

Model log-linear ini memasukkan efek utama dan interaksi dua arah antara X dan Y, dan X dan Z.

17.3.3.1 Pembentukan Model

model_conditional_X <- glm(counts ~ x.sex + y.fav + z.fund +
                             x.sex*y.fav + x.sex*z.fund,
                           family = poisson(link = "log"))
summary(model_conditional_X)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     x.sex * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          4.23495    0.08955  47.293  < 2e-16 ***
## x.sex1M             -0.52960    0.13966  -3.792 0.000149 ***
## y.fav1fav            0.48302    0.08075   5.982 2.20e-09 ***
## z.fund1fund          0.07962    0.10309   0.772 0.439916    
## z.fund2mod           0.41097    0.09585   4.288 1.81e-05 ***
## x.sex1M:y.fav1fav    0.65845    0.12708   5.181 2.20e-07 ***
## x.sex1M:z.fund1fund -0.12841    0.15109  -0.850 0.395405    
## x.sex1M:z.fund2mod  -0.06267    0.13908  -0.451 0.652274    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   3.9303  on  4  degrees of freedom
## AIC: 96.067
## 
## Number of Fisher Scoring iterations: 4

17.3.3.2 Membandingkan Model Conditional pada X dengan Model Homogenous

Hipotesis :

\(H_0\) : Tidak ada interaksi antara Y dan Z (model conditional pada X sudah cukup)

\(H_1\) : Ada interaksi antara Y dan Z (model homogenous dibutuhkan)

Statistik uji : selisih deviance

deviance <- model_conditional_X$deviance - model_homogenous$deviance 
deviance

## [1] 2.132302

Daerah Kritis

df <- model_conditional_X$df.residual - model_homogenous$df.residual
df

## [1] 2

chi.tabel <- qchisq(1 - 0.05, df = df)
chi.tabel

## [1] 5.991465

Keputusan

keputusan <- ifelse(deviance <= chi.tabel, "Terima H0", "Tolak H0")
keputusan

## [1] "Terima H0"

Kesimpulan :

Dengan taraf signifikansi 5%, tidak ada interaksi antara Y dan Z sehingga model conditional pada X sudah cukup.

17.3.4 Model Conditional pada Y

Model log-linear ini memasukkan efek utama dan interaksi dua arah antara X dan Y, dan Y dan Z.

17.3.4.1 Pembentukan Model

model_conditional_Y <- glm(counts ~ x.sex + y.fav + z.fund +
x.sex*y.fav + y.fav*z.fund,
family = poisson(link = "log"))
summary(model_conditional_Y)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     y.fav * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            4.33931    0.09919  43.748  < 2e-16 ***
## x.sex1M               -0.59345    0.10645  -5.575 2.48e-08 ***
## y.fav1fav              0.37259    0.12438   2.996  0.00274 ** 
## z.fund1fund           -0.12516    0.13389  -0.935  0.34989    
## z.fund2mod             0.30228    0.12089   2.500  0.01240 *  
## x.sex1M:y.fav1fav      0.65845    0.12708   5.181 2.20e-07 ***
## y.fav1fav:z.fund1fund  0.21254    0.16205   1.312  0.18966    
## y.fav1fav:z.fund2mod   0.11757    0.14771   0.796  0.42606    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   2.9203  on  4  degrees of freedom
## AIC: 95.057
## 
## Number of Fisher Scoring iterations: 4

17.3.4.2 Membandingkan Model Conditional pada Y dengan Model Homogenous

Hipotesis :

\(H_0\) : Tidak ada interaksi antara X dan Z (model conditional pada Y sudah cukup)

\(H_1\) : Ada interaksi antara X dan Z (model homogenous dibutuhkan)

Statistik uji : selisih deviance

deviance <- model_conditional_Y$deviance - model_homogenous$deviance 
deviance

## [1] 1.122315

Daerah Kritis

df <- model_conditional_Y$df.residual - model_homogenous$df.residual
df

## [1] 2

chi.tabel <- qchisq(1 - 0.05, df = df)
chi.tabel

## [1] 5.991465

Keputusan

keputusan <- ifelse(deviance <= chi.tabel, "Terima H0", "Tolak H0")
keputusan

## [1] "Terima H0"

Kesimpulan :

Dengan taraf signifikansi 5%, tidak ada interaksi antara X dan Z sehingga model conditional pada Y sudah cukup.

17.3.5 Model Conditional pada Z

Model log-linear ini memasukkan efek utama dan interaksi dua arah antara X dan Z, dan Y dan Z.

17.3.5.1 Pembentukan Model

model_conditional_Z <- glm(counts ~ x.sex + y.fav + z.fund +
                             x.sex*z.fund + y.fav*z.fund,
                           family = poisson(link = "log"))
summary(model_conditional_Z)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * z.fund + 
##     y.fav * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            4.12255    0.10518  39.195  < 2e-16 ***
## x.sex1M               -0.07453    0.10713  -0.696    0.487    
## y.fav1fav              0.65896    0.11292   5.836 5.36e-09 ***
## z.fund1fund           -0.06540    0.15126  -0.432    0.665    
## z.fund2mod             0.33196    0.13777   2.410    0.016 *  
## x.sex1M:z.fund1fund   -0.12841    0.15109  -0.850    0.395    
## x.sex1M:z.fund2mod    -0.06267    0.13908  -0.451    0.652    
## y.fav1fav:z.fund1fund  0.21254    0.16205   1.312    0.190    
## y.fav1fav:z.fund2mod   0.11757    0.14771   0.796    0.426    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.361  on 11  degrees of freedom
## Residual deviance:  29.729  on  3  degrees of freedom
## AIC: 123.87
## 
## Number of Fisher Scoring iterations: 4

17.3.5.2 Membandingkan Model Conditional pada Y dengan Model Homogenous

Hipotesis :

\(H_0\) : Tidak ada interaksi antara X dan Z (model conditional pada Y sudah cukup)

\(H_1\) : Ada interaksi antara X dan Z (model homogenous dibutuhkan)

Statistik uji : selisih deviance

deviance <- model_conditional_Z$deviance - model_homogenous$deviance 
deviance

## [1] 27.93095

Daerah Kritis

df <- model_conditional_Z$df.residual - model_homogenous$df.residual
df

## [1] 1

chi.tabel <- qchisq(1 - 0.05, df = df)
chi.tabel

## [1] 3.841459

Keputusan

keputusan <- ifelse(deviance <= chi.tabel, "Terima H0", "Tolak H0")
keputusan

## [1] "Tolak H0"

Kesimpulan :

Dengan taraf signifikansi 5%, ada interaksi antara X dan Y sehingga model dengan interaksi antara X dan Y dibutuhkan

17.3.6 Ringkasan

\[ \begin{array}{l|l|c} \textbf{Model} & \textbf{Parameter} & \textbf{AIC} \\ \hline \text{Saturated} & \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} + \lambda_{ijk}^{XYZ} & 100.14 \\ \text{Homogenous} & \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} & 97.934 \\ \text{Conditional on X} & \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} & 96.067 \\ \text{Conditional on Y} & \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{jk}^{YZ} & 95.057 \\ \text{Conditional on Z} & \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} & 123.87 \\ \end{array} \] Ringkasan Pengujian Interaksi \[ \begin{array}{l|l|c} \textbf{Interaksi} & \textbf{Pengujian} & \textbf{Keterangan} \\ \hline \text{XYZ} & \text{Saturated vs Homogenous}& \text{tidak ada interaksi} \\ \text{YZ} & \text{Conditional pada X vs Homogenous}& \text{tidak ada interaksi} \\ \text{XZ} & \text{Conditional pada Y vs Homogenous}& \text{tidak ada interaksi} \\ \text{XY} & \text{Conditional pada Z vs Homogenous}& \text{ ada interaksi} \\ \end{array} \] Sehingga model terbaik dengan hanya asosaisi yang nyata adalah :

\[ \log(\mu_{ijk}) = \lambda + \lambda_{i}^{X} + \lambda_{j}^{Y} + \lambda_{k}^{Z} + \lambda_{ij}^{XY} \]

17.3.7 Model terbaik

bestmodel <- glm(counts ~ x.sex + y.fav + z.fund +
x.sex*y.fav,
family = poisson(link = "log"))
summary(bestmodel)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav, 
##     family = poisson(link = "log"))
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        4.26518    0.07794  54.721  < 2e-16 ***
## x.sex1M           -0.59345    0.10645  -5.575 2.48e-08 ***
## y.fav1fav          0.48302    0.08075   5.982 2.20e-09 ***
## z.fund1fund        0.01986    0.07533   0.264    0.792    
## z.fund2mod         0.38130    0.06944   5.491 4.00e-08 ***
## x.sex1M:y.fav1fav  0.65845    0.12708   5.181 2.20e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   4.6532  on  6  degrees of freedom
## AIC: 92.79
## 
## Number of Fisher Scoring iterations: 4

Interpretasi Koefisien

data.frame(
  koef = bestmodel$coefficients,
  exp_koef = exp(bestmodel$coefficients)
)

##                          koef   exp_koef
## (Intercept)        4.26517861 71.1776316
## x.sex1M           -0.59344782  0.5524194
## y.fav1fav          0.48302334  1.6209677
## z.fund1fund        0.01985881  1.0200573
## z.fund2mod         0.38129767  1.4641834
## x.sex1M:y.fav1fav  0.65845265  1.9318008

• Peluang seseorang berjenis kelamin laki-laki (tanpa memperhatikan fundamentalisme dan pendapat mengenai hukuman mati) adalah 0.55 kali dibandingkan perempuan

• Peluang seseorang mendukung hukuman mati (tanpa memperhatikan jenis kelamin atau fundamentalisme) adalah 1.62 kali dibandingkan yang menolak

• Peluang seseorang fundamentalist adalah 1.02 kali dibandingkan liberal dan peluang seseorang moderate adalah 1.46 kali dibandingkan liberal.

• Peluang jika hukuman mati jika dia laki-laki adalah 1.932 kali dibandingkan peluang yang sama jika dia perempuan

Nilai dugaan model terbaik

data.frame(
  Fund = z.fund,
  sex = x.sex,
  favor = y.fav,
  counts = counts,
  fitted = bestmodel$fitted.values
)

##     Fund sex favor counts    fitted
## 1  1fund  1M  1fav    128 125.59539
## 2  1fund  1M  2opp     32  40.10855
## 3  1fund  2F  1fav    123 117.69079
## 4  1fund  2F  2opp     73  72.60526
## 5   2mod  1M  1fav    182 180.27878
## 6   2mod  1M  2opp     56  57.57155
## 7   2mod  2F  1fav    168 168.93257
## 8   2mod  2F  2opp    105 104.21711
## 9   3lib  1M  1fav    119 123.12582
## 10  3lib  1M  2opp     49  39.31990
## 11  3lib  2F  1fav    111 115.37664
## 12  3lib  2F  2opp     70  71.17763

18 KASUS RIL I

Peneliti ingin meneliti apakah jumlah jam belajar dan hasil ujian tengah semester di Universitas Sumatera Utara program studi farmasi mempengaruhi status kelulusan mahasiswa. Data ini merupakan data yang diambil dari jurnal “Analisis Regresi Logistik Biner untuk Memprediksi Probabilitas Kelulusan Ujian Akhir Semester Mahasiswa yang Mengambil Mata Kuliah Matematika Farmasi” yang ditulis oleh Siti Fatimah Sihotang dan dipublikasikan dalam Journal of Mathematics Education and Science (MES), Vol. 8, No. 2, April 2023 yang merupakan hasil survey langsung pada mahasiswa di USU. Akan dilakukan analisis regresi logisitk biner untuk melihat apakah betul jam belajar dan hasil UTS mempengaruhi status kelulusan mahasiswa.

\[ \begin{array}{|c|c|c|c|} \hline \textbf{Mahasiswa} & \textbf{Kelulusan} & \textbf{Jam} & \textbf{UTS} \\ \hline 1 & 0 & 1 & 5 \\ 2 & 0 & 1 & 5 \\ 3 & 0 & 1 & 5 \\ 4 & 0 & 1 & 5 \\ 5 & 1 & 1 & 5 \\ 6 & 0 & 1 & 5 \\ 7 & 0 & 1 & 5 \\ 8 & 0 & 1 & 6 \\ 9 & 0 & 1 & 7 \\ 10 & 1 & 1 & 6 \\ 11 & 0 & 2 & 5 \\ 12 & 0 & 2 & 5 \\ 13 & 0 & 2 & 6 \\ 14 & 1 & 2 & 6 \\ 15 & 1 & 2 & 7 \\ 16 & 1 & 2 & 7 \\ 17 & 1 & 2 & 7 \\ 18 & 1 & 2 & 6 \\ 19 & 1 & 2 & 5 \\ 20 & 1 & 2 & 5 \\ 21 & 0 & 3 & 5 \\ 22 & 0 & 3 & 6 \\ 23 & 1 & 3 & 5 \\ 24 & 1 & 3 & 5 \\ 25 & 1 & 3 & 6 \\ 26 & 1 & 3 & 6 \\ 27 & 1 & 3 & 7 \\ 28 & 1 & 3 & 7 \\ 29 & 1 & 3 & 8 \\ 30 & 1 & 3 & 8 \\ 31 & 0 & 4 & 6 \\ 32 & 1 & 4 & 5 \\ 33 & 1 & 4 & 6 \\ 34 & 1 & 4 & 7 \\ 35 & 1 & 4 & 7 \\ 36 & 1 & 4 & 8 \\ 37 & 1 & 4 & 8 \\ 38 & 1 & 4 & 8 \\ 39 & 1 & 4 & 6 \\ 40 & 1 & 4 & 8 \\ \hline \end{array} \]

data_mahasiswa <- data.frame(
  Kelulusan = c(0,0,0,0,1,0,0,0,0,1,
                0,0,0,1,1,1,1,1,1,1,
                0,0,1,1,1,1,1,1,1,1,
                0,1,1,1,1,1,1,1,1,1), 
  Jam = as.factor(rep(1:4, each = 10)), 
  UTS = c(5,5,5,5,5,5,5,6,7,6,
          5,5,6,6,7,7,7,6,5,5,
          5,6,5,5,6,6,7,7,8,8,
          6,5,6,7,7,8,8,8,6,8)
)
n <- 40
attach(data_mahasiswa)

18.1 Membentuk Model Regresi Logistik

model_full <- glm(Kelulusan~Jam+UTS, data= data_mahasiswa, family=binomial)
summary(model_full)

## 
## Call:
## glm(formula = Kelulusan ~ Jam + UTS, family = binomial, data = data_mahasiswa)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -6.7280     3.0749  -2.188   0.0287 *
## Jam2          1.9969     1.1052   1.807   0.0708 .
## Jam3          2.3317     1.1825   1.972   0.0486 *
## Jam4          2.6467     1.4058   1.883   0.0597 .
## UTS           0.9677     0.5266   1.838   0.0661 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 51.796  on 39  degrees of freedom
## Residual deviance: 34.720  on 35  degrees of freedom
## AIC: 44.72
## 
## Number of Fisher Scoring iterations: 5

exp(model_full$coefficients)

##  (Intercept)         Jam2         Jam3         Jam4          UTS 
##  0.001196894  7.366385812 10.295886879 14.106959425  2.631971272

Interpretasi koefisien :

• Intercept (-6.728) : log-odds dasar untuk mahasiswa dengan jam belajar kategori 1 (4-6 jam) dan uts 0

• Jam 2 (1.9969) : mahasiswa dengan jam belajar kategori 2 (7-9 jam) memiliki log-odds lulus 1.9969 lebih tinggi dibandingkan mahasiswa dengan jam belajar kategori 1. Nilai odds rationya adalah 7.366 yang menunjukkan bahwa mahasiswa dengan jam belajar kategori 2 memiliki peluang 7 kali lebih tinggi dibandingkan jam belajar kategori 1, tetapi hal ini tidak signifikan (mendekati signifikan)

• Jam 3 (2.3317) : mahasiswa dengan jam belajar kategori 3 (10-12 jam) memiliki log-odds lulus 2.3317 lebih tinggi dibandingkan mahasiswa dengan jam belajar kategori 1. Mahasiswa dengan jam belajar 10-12 jam memiliki peluang lulus 10 kali dibandingkan jam belajar 4-6 jam

• Jam 4 (2.6467) : mahasiswa dengan jam belajar kategori 4 memiliki log-odds lulus 2.6467 lebih tinggi dibandingkan mahasiswa dengan jam belajar kategori 1. Namun hal ini tidak signifikan (mendekati signifikan )

• UTS (0.9677) : setiap kenaikan 1 satuan nilai UTS akan meningkatkan 0.9677 log-oddsnya.

• Namun dari faktor-faktor yang ada , hanya kelompok jam belajar 3 yang signifikan sehingga dapat pertimbangkan untuk mengurangi faktornya.

model_jam <- glm(Kelulusan~Jam, data= data_mahasiswa, family=binomial)
summary(model_jam)

## 
## Call:
## glm(formula = Kelulusan ~ Jam, family = binomial, data = data_mahasiswa)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -1.3863     0.7906  -1.754  0.07951 . 
## Jam2          2.2336     1.0494   2.128  0.03330 * 
## Jam3          2.7726     1.1180   2.480  0.01314 * 
## Jam4          3.5835     1.3175   2.720  0.00653 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 51.796  on 39  degrees of freedom
## Residual deviance: 38.735  on 36  degrees of freedom
## AIC: 46.735
## 
## Number of Fisher Scoring iterations: 4

18.2 Pemilihan Model

list (
  AIC_ModelFull = AIC (model_full), 
  AIC_ModelJam = AIC(model_jam), 
  Deviance_Full = deviance(model_full), 
  Deviance_Jam = deviance(model_jam)
)

## $AIC_ModelFull
## [1] 44.71992
## 
## $AIC_ModelJam
## [1] 46.73504
## 
## $Deviance_Full
## [1] 34.71992
## 
## $Deviance_Jam
## [1] 38.73504

anova(model_full,model_jam, test = "LRT")

## Analysis of Deviance Table
## 
## Model 1: Kelulusan ~ Jam + UTS
## Model 2: Kelulusan ~ Jam
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1        35     34.720                       
## 2        36     38.735 -1  -4.0151  0.04509 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dari hasil nilai AIC, didapatkan bahwa model full memiliki nilai AIC yang lebih kecil sehingga sebenarnya model full lebih baik dibandingkan model yang hanya mengambil jam sebagai faktor. Selain itu deviance juga mendukung pernyataan sebelumnya. Dengan likelihood test, model lebih baik menggunakan uts sebagai salah satu faktor untuk menentukan tingkat kelulusan.

18.3 Goodness-of-Fit-Model

model_null <- glm(Kelulusan ~ 1, data = data_mahasiswa, family = binomial)

logL0 <- logLik(model_null)
logLM <- logLik(model_full)
L0 <- exp(logL0)
LM <- exp(logLM)
n <- nobs(model_full)
cox_snell <- 1 - (L0 / LM)^(2 / n)
mcfadden <- 1 - (as.numeric(logLM) / as.numeric(logL0))
r2 <- data.frame(
R2_Cox_Snell = cox_snell,
R2_McFadden = mcfadden
)
r2

##   R2_Cox_Snell R2_McFadden
## 1    0.3474681    0.329676

Model ini sudah cukup menjelaskan data dengan baik karena nilai R-squared McFadden jika sudah lebih besar dari 0.2 sudah cukup menjelaskan bahwa model memiliki kecocokan yang baik dan juga R-squared Cox & Snell juga sama

18.4 Confusion Matrix

library(caret)
pred_prob <- predict(model_full, type = "response")
pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
conf_matrix <- confusionMatrix(factor(pred_class), factor(data_mahasiswa$Kelulusan), positive = "1")
conf_matrix

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0  7  2
##          1  7 24
##                                           
##                Accuracy : 0.775           
##                  95% CI : (0.6155, 0.8916)
##     No Information Rate : 0.65            
##     P-Value [Acc > NIR] : 0.06444         
##                                           
##                   Kappa : 0.4611          
##                                           
##  Mcnemar's Test P-Value : 0.18242         
##                                           
##             Sensitivity : 0.9231          
##             Specificity : 0.5000          
##          Pos Pred Value : 0.7742          
##          Neg Pred Value : 0.7778          
##              Prevalence : 0.6500          
##          Detection Rate : 0.6000          
##    Detection Prevalence : 0.7750          
##       Balanced Accuracy : 0.7115          
##                                           
##        'Positive' Class : 1               
## 

Dari confusion matrix didapatkan bahwa model ini dapat menjelaskan secara akurat 71.15%. Model juga lebih baik dalam menjelaskan peluang lulus dibandingkan peluang tidak lulus.

18.5 Kesimpulan

Hasil dari pengujian dan pembentukan model didapatkan bahwa faktor-faktor yang mempengaruhi peluang seseorang untuk lulus adalah jam belajar dan nilai UTS. Model nya adalah sebagai berikut :

\[ \log\left(\frac{P(Kelulusan =1)}{1-P(Kelulusan=1)}\right) = -6.7280 + 1.9969. Jam_2 + 2.3317.Jam_3 + 2.6467.Jam_4 + 0.9677. UTS \]

Model ini sudah cukup menjelaskan data dengan baik dan dapat digunakan untuk mmengetahui peluang lulusnya.

19 KASUS RIL II

Akan dibantuk model log-linear dari data jumlah kecelakaan lalu lintas di kota Makassar dengan tiga jenis variabel yaitu : jenis kendaraan (1 : kendaraan roda 2 dan 2: kendaraan roda 4) , umur, dan pendidikan terakhir. Datanya adalah sebagai berikut

\[ \begin{array}{cc|cccc|c} \textbf{Jenis} & \textbf{Umur} & \textbf{SD} & \textbf{SMP} & \textbf{SMA} & \textbf{PT} & \textbf{Total} \\ \hline 1 & \text{Anak-Anak} & 1 & 2 & 0 & 0 & 3 \\ 1 & \text{Remaja} & 1 & 5 & 10 & 0 & 16 \\ 1 & \text{Dewasa} & 7 & 17 & 77 & 25 & 126 \\ 2 & \text{Anak-Anak} & 0 & 0 & 0 & 0 & 0 \\ 2 & \text{Remaja} & 0 & 0 & 2 & 0 & 2 \\ 2 & \text{Dewasa} & 0 & 0 & 28 & 7 & 35 \\ \hline \text{}&{\textbf{Total}} & 9 & 24 & 117 & 32 & 182 \\ \end{array} \]

# Buat data frame
data_loglin <- data.frame(
  Jenis = c(
    rep(1, 12),  # Jenis 1: 3 Umur × 4 Pendidikan
    rep(2, 12)   # Jenis 2: 3 Umur × 4 Pendidikan
  ),
  Umur = rep(rep(c("Anak", "Remaja", "Dewasa"), each = 4), 2),
  Pendidikan = rep(c("SD", "SMP", "SMA", "PT"), 6),
  Count = c(
    # Jenis 1
    1, 2, 0, 0,     # Anak-Anak
    1, 5, 10, 0,    # Remaja
    7, 17, 77, 25,  # Dewasa
    # Jenis 2
    0, 0, 0, 0,     # Anak-Anak
    0, 0, 2, 0,     # Remaja
    0, 0, 28, 7     # Dewasa
  )
)

19.1 Membentuk Model

19.1.1 Model Saturated

Model saturated adalah model yang menyertakan semua faktor utama, interaksi dau arah, dan interaksi tiga arah

model_saturated <- glm(Count~Jenis*Umur*Pendidikan, data = data_loglin, family=poisson())
summary(model_saturated)

## 
## Call:
## glm(formula = Count ~ Jenis * Umur * Pendidikan, family = poisson(), 
##     data = data_loglin)
## 
## Coefficients:
##                                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)                    -2.530e+01  4.234e+05       0        1
## Jenis                           6.840e-09  2.678e+05       0        1
## UmurDewasa                      2.979e+01  4.234e+05       0        1
## UmurRemaja                      1.382e-08  5.987e+05       0        1
## PendidikanSD                    5.061e+01  4.638e+05       0        1
## PendidikanSMA                   5.204e-08  5.987e+05       0        1
## PendidikanSMP                   5.199e+01  4.638e+05       0        1
## Jenis:UmurDewasa               -1.273e+00  2.678e+05       0        1
## Jenis:UmurRemaja               -6.651e-09  3.787e+05       0        1
## Jenis:PendidikanSD             -2.530e+01  3.279e+05       0        1
## Jenis:PendidikanSMA            -2.506e-08  3.787e+05       0        1
## Jenis:PendidikanSMP            -2.600e+01  3.279e+05       0        1
## UmurDewasa:PendidikanSD        -2.590e+01  5.009e+05       0        1
## UmurRemaja:PendidikanSD        -1.389e-08  6.559e+05       0        1
## UmurDewasa:PendidikanSMA        8.636e-01  5.987e+05       0        1
## UmurRemaja:PendidikanSMA        2.921e+01  7.333e+05       0        1
## UmurDewasa:PendidikanSMP       -2.551e+01  5.009e+05       0        1
## UmurRemaja:PendidikanSMP        1.833e+00  6.559e+05       0        1
## Jenis:UmurDewasa:PendidikanSD  -6.729e-01  3.787e+05       0        1
## Jenis:UmurRemaja:PendidikanSD   6.723e-09  4.638e+05       0        1
## Jenis:UmurDewasa:PendidikanSMA  2.614e-01  3.787e+05       0        1
## Jenis:UmurRemaja:PendidikanSMA -1.609e+00  4.638e+05       0        1
## Jenis:UmurDewasa:PendidikanSMP -8.671e-01  3.787e+05       0        1
## Jenis:UmurRemaja:PendidikanSMP -9.163e-01  4.638e+05       0        1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4.9755e+02  on 23  degrees of freedom
## Residual deviance: 2.4629e-10  on  0  degrees of freedom
## AIC: 93.584
## 
## Number of Fisher Scoring iterations: 23

19.1.2 Model Homogenous

Model homogenous adalah model yang menyertakan semua faktor utama dan interaksi dua arah,.

model_homogenous <- glm(Count ~ Jenis + Umur + Pendidikan +
                          Jenis*Umur + Jenis*Pendidikan + Umur*Pendidikan,
                       data = data_loglin, family = poisson())
summary(model_homogenous)

## 
## Call:
## glm(formula = Count ~ Jenis + Umur + Pendidikan + Jenis * Umur + 
##     Jenis * Pendidikan + Umur * Pendidikan, family = poisson(), 
##     data = data_loglin)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)              -2.502e+01  2.184e+05   0.000    1.000
## Jenis                    -1.867e-01  1.189e+05   0.000    1.000
## UmurDewasa                2.952e+01  2.184e+05   0.000    1.000
## UmurRemaja                2.461e+00  2.464e+05   0.000    1.000
## PendidikanSD              5.013e+01  1.839e+05   0.000    1.000
## PendidikanSMA            -3.909e-01  1.893e+05   0.000    1.000
## PendidikanSMP             5.185e+01  1.844e+05   0.000    1.000
## Jenis:UmurDewasa         -1.086e+00  1.189e+05   0.000    1.000
## Jenis:UmurRemaja         -1.684e+00  1.189e+05   0.000    1.000
## Jenis:PendidikanSD       -2.492e+01  9.721e+04   0.000    1.000
## Jenis:PendidikanSMA       2.614e-01  4.812e-01   0.543    0.587
## Jenis:PendidikanSMP      -2.595e+01  1.013e+05   0.000    1.000
## UmurDewasa:PendidikanSD  -2.648e+01  1.441e+05   0.000    1.000
## UmurRemaja:PendidikanSD  -7.771e-01  1.838e+05   0.000    1.000
## UmurDewasa:PendidikanSMA  1.254e+00  1.893e+05   0.000    1.000
## UmurRemaja:PendidikanSMA  2.687e+01  2.211e+05   0.000    1.000
## UmurDewasa:PendidikanSMP -2.629e+01  1.441e+05   0.000    1.000
## UmurRemaja:PendidikanSMP  1.392e-01  1.838e+05   0.000    1.000
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4.9755e+02  on 23  degrees of freedom
## Residual deviance: 3.0449e-10  on  6  degrees of freedom
## AIC: 81.584
## 
## Number of Fisher Scoring iterations: 23

exp(model_homogenous$coefficients)

##              (Intercept)                    Jenis               UmurDewasa 
##             1.356213e-11             8.297009e-01             6.583459e+12 
##               UmurRemaja             PendidikanSD            PendidikanSMA 
##             1.171888e+01             5.922488e+21             6.764737e-01 
##            PendidikanSMP         Jenis:UmurDewasa         Jenis:UmurRemaja 
##             3.301813e+22             3.374710e-01             1.856090e-01 
##       Jenis:PendidikanSD      Jenis:PendidikanSMA      Jenis:PendidikanSMP 
##             1.500536e-11             1.298701e+00             5.383046e-12 
##  UmurDewasa:PendidikanSD  UmurRemaja:PendidikanSD UmurDewasa:PendidikanSMA 
##             3.150703e-12             4.597427e-01             3.505827e+00 
## UmurRemaja:PendidikanSMA UmurDewasa:PendidikanSMP UmurRemaja:PendidikanSMP 
##             4.650560e+11             3.825854e-12             1.149357e+00

19.1.3 Model Conditional

Model conditional pada X

model_conditional_X <- glm(Count ~ Jenis + Umur + Pendidikan +
                          Jenis*Umur + Jenis*Pendidikan,
                       data = data_loglin, family = poisson())

## Warning: glm.fit: fitted rates numerically 0 occurred

summary(model_conditional_X)

## 
## Call:
## glm(formula = Count ~ Jenis + Umur + Pendidikan + Jenis * Umur + 
##     Jenis * Pendidikan, family = poisson(), data = data_loglin)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)  
## (Intercept)           18.8628  6362.2761   0.003   0.9976  
## Jenis                -19.5220  6362.2761  -0.003   0.9976  
## UmurDewasa           -14.5963  6362.2761  -0.002   0.9982  
## UmurRemaja           -15.8615  6362.2762  -0.002   0.9980  
## PendidikanSD          18.8486  7887.3432   0.002   0.9981  
## PendidikanSMA          1.0388     0.6182   1.680   0.0929 .
## PendidikanSMP         20.8103  7887.3433   0.003   0.9979  
## Jenis:UmurDewasa      18.3340  6362.2760   0.003   0.9977  
## Jenis:UmurRemaja      17.5355  6362.2761   0.003   0.9978  
## Jenis:PendidikanSD   -19.8703  7887.3432  -0.003   0.9980  
## Jenis:PendidikanSMA    0.2083     0.4772   0.436   0.6625  
## Jenis:PendidikanSMP  -20.8511  7887.3433  -0.003   0.9979  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 497.552  on 23  degrees of freedom
## Residual deviance:  18.332  on 12  degrees of freedom
## AIC: 87.917
## 
## Number of Fisher Scoring iterations: 19

exp(model_conditional_X$coefficients)

##         (Intercept)               Jenis          UmurDewasa          UmurRemaja 
##        1.556000e+08        3.324174e-09        4.580482e-07        1.292552e-07 
##        PendidikanSD       PendidikanSMA       PendidikanSMP    Jenis:UmurDewasa 
##        1.534119e+08        2.825760e+00        1.090929e+09        9.169340e+07 
##    Jenis:UmurRemaja  Jenis:PendidikanSD Jenis:PendidikanSMA Jenis:PendidikanSMP 
##        4.126203e+07        2.346624e-09        1.231527e+00        8.799840e-10

Model Conditional Y

model_conditional_Y <- glm(Count ~ Jenis + Umur + Pendidikan +
                          Jenis*Umur + Umur*Pendidikan,
                       data = data_loglin, family = poisson())

## Warning: glm.fit: fitted rates numerically 0 occurred

summary(model_conditional_Y)

## 
## Call:
## glm(formula = Count ~ Jenis + Umur + Pendidikan + Jenis * Umur + 
##     Umur * Pendidikan, family = poisson(), data = data_loglin)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)               3.540e-01  9.716e+03   0.000    1.000
## Jenis                    -1.933e+01  5.513e+03  -0.004    0.997
## UmurDewasa                4.148e+00  9.716e+03   0.000    1.000
## UmurRemaja               -1.764e+01  1.335e+04  -0.001    0.999
## PendidikanSD              1.897e+01  8.000e+03   0.002    0.998
## PendidikanSMA            -2.231e-07  1.131e+04   0.000    1.000
## PendidikanSMP             1.967e+01  8.000e+03   0.002    0.998
## Jenis:UmurDewasa          1.805e+01  5.513e+03   0.003    0.997
## Jenis:UmurRemaja          1.725e+01  5.513e+03   0.003    0.998
## UmurDewasa:PendidikanSD  -2.049e+01  8.000e+03  -0.003    0.998
## UmurRemaja:PendidikanSD   2.683e-01  1.215e+04   0.000    1.000
## UmurDewasa:PendidikanSMA  1.188e+00  1.131e+04   0.000    1.000
## UmurRemaja:PendidikanSMA  2.173e+01  1.455e+04   0.001    0.999
## UmurDewasa:PendidikanSMP -2.030e+01  8.000e+03  -0.003    0.998
## UmurRemaja:PendidikanSMP  1.185e+00  1.215e+04   0.000    1.000
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 497.552  on 23  degrees of freedom
## Residual deviance:  14.937  on  9  degrees of freedom
## AIC: 90.521
## 
## Number of Fisher Scoring iterations: 18

exp(model_conditional_Y$coefficients)

##              (Intercept)                    Jenis               UmurDewasa 
##             1.424784e+00             4.034343e-09             6.327734e+01 
##               UmurRemaja             PendidikanSD            PendidikanSMA 
##             2.193684e-08             1.739716e+08             9.999998e-01 
##            PendidikanSMP         Jenis:UmurDewasa         Jenis:UmurRemaja 
##             3.479431e+08             6.885329e+07             3.098398e+07 
##  UmurDewasa:PendidikanSD  UmurRemaja:PendidikanSD UmurDewasa:PendidikanSMA 
##             1.257389e-09             1.307785e+00             3.281251e+00 
## UmurRemaja:PendidikanSMA UmurDewasa:PendidikanSMP UmurRemaja:PendidikanSMP 
##             2.730209e+09             1.526830e-09             3.269461e+00

Model Conditional pada Z

model_conditional_Z <- glm(Count ~ Jenis + Umur + Pendidikan +
                           Jenis*Pendidikan + Umur*Pendidikan,
                       data = data_loglin, family = poisson())
summary(model_conditional_Z)

## 
## Call:
## glm(formula = Count ~ Jenis + Umur + Pendidikan + Jenis * Pendidikan + 
##     Umur * Pendidikan, family = poisson(), data = data_loglin)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)   
## (Intercept)              -2.140e+01  4.493e+04   0.000  0.99962   
## Jenis                    -1.273e+00  4.276e-01  -2.977  0.00291 **
## UmurDewasa                2.589e+01  4.493e+04   0.001  0.99954   
## UmurRemaja               -5.375e-07  6.353e+04   0.000  1.00000   
## PendidikanSD              4.538e+01  5.554e+04   0.001  0.99935   
## PendidikanSMA            -3.083e-01  6.435e+04   0.000  1.00000   
## PendidikanSMP             4.712e+01  5.617e+04   0.001  0.99933   
## Jenis:PendidikanSD       -2.271e+01  3.266e+04  -0.001  0.99945   
## Jenis:PendidikanSMA       2.083e-01  4.772e-01   0.436  0.66251   
## Jenis:PendidikanSMP      -2.376e+01  3.371e+04  -0.001  0.99944   
## UmurDewasa:PendidikanSD  -2.395e+01  4.493e+04  -0.001  0.99957   
## UmurRemaja:PendidikanSD   5.375e-07  6.353e+04   0.000  1.00000   
## UmurDewasa:PendidikanSMA  1.239e+00  6.435e+04   0.000  0.99998   
## UmurRemaja:PendidikanSMA  2.496e+01  7.848e+04   0.000  0.99975   
## UmurDewasa:PendidikanSMP -2.375e+01  4.493e+04  -0.001  0.99958   
## UmurRemaja:PendidikanSMP  9.163e-01  6.353e+04   0.000  0.99999   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 497.55244  on 23  degrees of freedom
## Residual deviance:   0.61319  on  8  degrees of freedom
## AIC: 78.198
## 
## Number of Fisher Scoring iterations: 21

exp(model_conditional_Z$coefficients)

##              (Intercept)                    Jenis               UmurDewasa 
##             5.085648e-10             2.800000e-01             1.755641e+11 
##               UmurRemaja             PendidikanSD            PendidikanSMA 
##             9.999995e-01             5.131225e+19             7.347145e-01 
##            PendidikanSMP       Jenis:PendidikanSD      Jenis:PendidikanSMA 
##             2.915103e+20             1.368594e-10             1.231527e+00 
##      Jenis:PendidikanSMP  UmurDewasa:PendidikanSD  UmurRemaja:PendidikanSD 
##             4.818056e-11             3.987148e-11             1.000001e+00 
## UmurDewasa:PendidikanSMA UmurRemaja:PendidikanSMA UmurDewasa:PendidikanSMP 
##             3.451597e+00             6.925449e+10             4.841537e-11 
## UmurRemaja:PendidikanSMP 
##             2.500001e+00

19.1.4 Model Tanpa Interaksi

model_tanpa_interaksi <- glm(Count ~ Jenis + Umur + Pendidikan,
                       data = data_loglin, family = poisson())
summary(model_tanpa_interaksi)

## 
## Call:
## glm(formula = Count ~ Jenis + Umur + Pendidikan, family = poisson(), 
##     data = data_loglin)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     0.4989     0.6389   0.781 0.434894    
## Jenis          -1.3658     0.1842  -7.416 1.21e-13 ***
## UmurDewasa      3.9828     0.5827   6.835 8.19e-12 ***
## UmurRemaja      1.7918     0.6236   2.873 0.004063 ** 
## PendidikanSD   -1.2685     0.3773  -3.362 0.000774 ***
## PendidikanSMA   1.2964     0.1995   6.499 8.10e-11 ***
## PendidikanSMP  -0.2877     0.2700  -1.065 0.286710    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 497.552  on 23  degrees of freedom
## Residual deviance:  37.946  on 17  degrees of freedom
## AIC: 97.53
## 
## Number of Fisher Scoring iterations: 6

exp(model_tanpa_interaksi$coefficients)

##   (Intercept)         Jenis    UmurDewasa    UmurRemaja  PendidikanSD 
##     1.6468830     0.2551724    53.6666666     6.0000000     0.2812500 
## PendidikanSMA PendidikanSMP 
##     3.6562500     0.7500000

19.2 Pengujian Interaksi

19.2.1 Homogenous vs Saturated

deviance <- model_homogenous$deviance - model_saturated$deviance
deviance

## [1] 5.820315e-11

df <- model_homogenous$df.residual - model_saturated$df.residual
df

## [1] 6

chi.tabel <- qchisq(1 - 0.05, df = df)
chi.tabel

## [1] 12.59159

keputusan <- ifelse(deviance <= chi.tabel, "Terima H0", "Tolak H0")
keputusan

## [1] "Terima H0"

19.2.2 Homogenous vs Conditional X

deviance <- model_conditional_X$deviance - model_homogenous$deviance
deviance

## [1] 18.33246

df <- model_conditional_X$df.residual - model_homogenous$df.residual
df

## [1] 6

chi.tabel <- qchisq(1 - 0.05, df = df)
chi.tabel

## [1] 12.59159

keputusan <- ifelse(deviance <= chi.tabel, "Terima H0", "Tolak H0")
keputusan

## [1] "Tolak H0"

19.2.3 Homogenous vs Conditional on Y

deviance <- model_conditional_Y$deviance - model_homogenous$deviance
deviance

## [1] 14.93654

df <- model_conditional_Y$df.residual - model_homogenous$df.residual
df

## [1] 3

chi.tabel <- qchisq(1 - 0.05, df = df)
chi.tabel

## [1] 7.814728

keputusan <- ifelse(deviance <= chi.tabel, "Terima H0", "Tolak H0")
keputusan

## [1] "Tolak H0"

19.2.4 Homogenous vs Conditional Z

deviance <- model_conditional_Z$deviance - model_homogenous$deviance
deviance

## [1] 0.6131902

df <- model_conditional_Z$df.residual - model_homogenous$df.residual
df

## [1] 2

chi.tabel <- qchisq(1 - 0.05, df = df)
chi.tabel

## [1] 5.991465

keputusan <- ifelse(deviance <= chi.tabel, "Terima H0", "Tolak H0")
keputusan

## [1] "Terima H0"

19.2.5 Model Conditional Z vs Model Tanpa Interaksi

deviance <- model_tanpa_interaksi$deviance - model_conditional_Z$deviance
deviance

## [1] 37.33255

df <- model_tanpa_interaksi$df.residual - model_conditional_Z$df.residual
df

## [1] 9

chi.tabel <- qchisq(1 - 0.05, df = df)
chi.tabel

## [1] 16.91898

keputusan <- ifelse(deviance <= chi.tabel, "Terima H0", "Tolak H0")
keputusan

## [1] "Tolak H0"

Dari hasil pengujian didapatkan bahwa terdapat interkasi antara variabel X dan Y sehingga model log-linearnya adalah sebagai berikut :

\[ \log (\mu_{ijk}) = \lambda + \lambda_{i}^{X} + \lambda_{j}^{Y} + \lambda_{k}^{Z} + \lambda_{ij}^{XY} \]

19.3 Model Terbaik

Model terbaiknya adalah model dengan faktor utama dan interaksi antara X dan Y saja.

model_terbaik <- glm(Count~Jenis+Umur+Pendidikan+
                       Jenis*Umur, data= data_loglin, 
                     family = poisson())
summary(model_terbaik)

## 
## Call:
## glm(formula = Count ~ Jenis + Umur + Pendidikan + Jenis * Umur, 
##     family = poisson(), data = data_loglin)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        17.0183  2391.4543   0.007 0.994322    
## Jenis             -17.6579  2391.4541  -0.007 0.994109    
## UmurDewasa        -12.6393  2391.4543  -0.005 0.995783    
## UmurRemaja        -13.9045  2391.4545  -0.006 0.995361    
## PendidikanSD       -1.2685     0.3773  -3.362 0.000774 ***
## PendidikanSMA       1.2964     0.1995   6.499  8.1e-11 ***
## PendidikanSMP      -0.2877     0.2700  -1.065 0.286710    
## Jenis:UmurDewasa   16.3770  2391.4541   0.007 0.994536    
## Jenis:UmurRemaja   15.5785  2391.4542   0.007 0.994802    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 497.552  on 23  degrees of freedom
## Residual deviance:  35.301  on 15  degrees of freedom
## AIC: 98.885
## 
## Number of Fisher Scoring iterations: 16

exp(model_terbaik$coefficients)

##      (Intercept)            Jenis       UmurDewasa       UmurRemaja 
##     2.460026e+07     2.144175e-08     3.241993e-06     9.148480e-07 
##     PendidikanSD    PendidikanSMA    PendidikanSMP Jenis:UmurDewasa 
##     2.812500e-01     3.656250e+00     7.500000e-01     1.295500e+07 
## Jenis:UmurRemaja 
##     5.829748e+06

data.frame(
  Jenis = data_loglin$Jenis, 
  Umur=data_loglin$Umur,
  Pendidikan=data_loglin$Pendidikan,
  Count=data_loglin$Count,
  fitted = model_terbaik$fitted.values
)

##    Jenis   Umur Pendidikan Count       fitted
## 1      1   Anak         SD     1 1.483516e-01
## 2      1   Anak        SMP     2 3.956044e-01
## 3      1   Anak        SMA     0 1.928571e+00
## 4      1   Anak         PT     0 5.274725e-01
## 5      1 Remaja         SD     1 7.912088e-01
## 6      1 Remaja        SMP     5 2.109890e+00
## 7      1 Remaja        SMA    10 1.028571e+01
## 8      1 Remaja         PT     0 2.813187e+00
## 9      1 Dewasa         SD     7 6.230769e+00
## 10     1 Dewasa        SMP    17 1.661538e+01
## 11     1 Dewasa        SMA    77 8.100000e+01
## 12     1 Dewasa         PT    25 2.215385e+01
## 13     2   Anak         SD     0 3.180919e-09
## 14     2   Anak        SMP     0 8.482450e-09
## 15     2   Anak        SMA     0 4.135195e-08
## 16     2   Anak         PT     0 1.130993e-08
## 17     2 Remaja         SD     0 9.890110e-02
## 18     2 Remaja        SMP     0 2.637363e-01
## 19     2 Remaja        SMA     2 1.285714e+00
## 20     2 Remaja         PT     0 3.516484e-01
## 21     2 Dewasa         SD     0 1.730769e+00
## 22     2 Dewasa        SMP     0 4.615385e+00
## 23     2 Dewasa        SMA    28 2.250000e+01
## 24     2 Dewasa         PT     7 6.153846e+00

20 REFERENSI

• Fisher, R. A. (1935). The Design of Experiments. Oliver and Boyd

• Agresti, A. (2013). Categorical Data Analysis. John Wiley & Sons.

• Everitt, B. S. (1992). The analysis of contingency tables (2nd ed.). Chapman & Hall.

• Mehta, C. R., & Patel, N. R. (1983). A Network Algorithm for Performing Fisher’s Exact Test in r × c Contingency Tables. Journal of the American Statistical Association, 78(382), 427–434.

• Howell, D. C. (2012). Statistical Methods for Psychology (8th ed.). Cengage Learning.

• Nely. (2019). Korelasi Hasil Belajar Siswa Tentang Materi Shalat Fardhu Terhadap Pengamalan Ibadah Siswa Sehari-Hari di SMP Negeri 6 Pelaihari (Skripsi, UIN Antasari Banjarmasin). Tersedia di repositori resmi UIN Antasari Banjarmasin.

• Nugraha, J. (2013). Pengantar analisis data kategorik: Metode dan aplikasi menggunakan program R. Deepublish.

• Maryana. (2013). Model Log Linier yang Terbaik untuk Analisis Data Kualitatif pada Tabel Kontingensi Tiga Arah. Malikussaleh Industrial Engineering Journal, 2(2), 32–37. ISSN 2302-934X.

• Sihotang, S. F. (2023). Analisis regresi logistik biner untuk memprediksi probabilitas kelulusan ujian akhir semester mahasiswa yang mengambil mata kuliah matematika farmasi. Journal of Mathematics Education and Science (MES), 8(2), April. https://jurnal.uisu.ac.id/index.php/mesuisu

• Nurman, T. A. (2011). Analisis data kategori dengan log linier menggunakan prinsip hirarki (studi kasus jumlah kecelakaan lalu lintas di Kota Makassar tahun 2011). UIN Alauddin Makassar.