E-BOOK Analisis Data Kategorik 1

KATA PENGANTAR

Puji syukur saya panjatkan atas terselesaikannya eBook yang berjudul “Analisis Data Kategorik” ini. Buku ini saya susun sebagai panduan praktis dan komprehensif untuk membantu pembaca memahami serta menguasai teknik-teknik analisis data kategorik, sebuah jenis data kualitatif yang sangat umum dijumpai dalam berbagai bidang ilmu dan penelitian.

Data kategorik, yang merepresentasikan kelompok atau kategori seperti jenis kelamin, pekerjaan, atau preferensi, memiliki karakteristik yang berbeda dengan data numerik. Oleh karena itu, analisisnya memerlukan pendekatan khusus agar dapat mengungkap pola distribusi, hubungan antar variabel, serta membangun model prediksi yang akurat. Dalam buku ini, saya memaparkan berbagai metode analisis data kategorik seperti tabel kontingensi, uji chi-square, regresi logistik, dan model log-linear, yang telah banyak digunakan dan dikembangkan dalam literatur statistik.

Saya berharap buku ini dapat menjadi sumber belajar yang bermanfaat bagi mahasiswa, peneliti, dan praktisi yang ingin memperdalam kemampuan analisis data non-numerik. Semoga materi yang disajikan dapat membantu meningkatkan kualitas penelitian dan pengambilan keputusan berbasis data.

Akhir kata, saya mengucapkan terima kasih kepada semua pihak yang telah mendukung tersusunnya buku ini. Kritik dan saran yang membangun sangat saya harapkan demi kesempurnaan buku ini di masa mendatang.

Selamat membaca

Jatinangor, 17 Juni 2025

Hadeelina Shakira Zalianty

1 Pendahuluan

Data kategorik adalah jenis data kualitatif yang merepresentasikan kelompok atau kategori, seperti jenis kelamin, pekerjaan, atau preferensi. Berbeda dari data numerik, data ini tidak memiliki nilai angka yang dapat diolah secara langsung dengan metode statistik biasa. Oleh karena itu, diperlukan pendekatan khusus dalam analisisnya. Analisis data kategorik berguna untuk memahami pola distribusi kategori, mengevaluasi hubungan antar variabel, serta membuat model prediksi berbasis klasifikasi. Beberapa metode umum yang digunakan adalah tabel kontingensi, uji chi-square, regresi logistik, dan model log-linear. Seperti dijelaskan oleh Agresti (2013), analisis data kategorik memerlukan teknik tersendiri karena karakteristiknya yang tidak memenuhi asumsi data numerik. Oleh sebab itu, penguasaan metode ini penting bagi peneliti yang bekerja dengan data non-numerik.

1.1 Tujuan Analisis Data Kategori

1.1.1 Memahami Distribusi Kategori dalam Data

Analisis data kategorik bertujuan untuk mengidentifikasi danmemahami bagaimana data tersebar di antara berbagai kategori. Hal ini membantu peneliti mengetahui proporsi atau frekuensi masing-masing kategori dalam dataset.

1.1.2 Menganalisis Hubungan Antar Variabel Kategorik

Tujuan lainnya adalah mengevaluasi apakah terdapat hubungan atau asosiasi antara dua atau lebih variabel kategorik. Misalnya, menentukan apakah terdapat keterkaitan antara tingkat pendidikan dan status pekerjaan.

1.1.3 Membangun Model Prediktif untuk Variabel Kategorik

Analisis data kategorik digunakan untuk membangun model yang dapat memprediksi kategori suatu variabel berdasarkan variabel lain. Contohnya, menggunakan regresi logistik untuk memprediksi apakah seorang pelanggan akan membeli produk berdasarkan karakteristik demografisnya.

1.1.4 Mengidentifikasi Pola dan Tren dalam Data Kategorik Dengan

menganalisis data kategorik, peneliti dapat mengungkap pola atau tren tertentu yang mungkin tidak terlihat secara langsung, seperti preferensi konsumen terhadap produk tertentu berdasarkan kelompok usia.

1.1.5 Menguji Hipotesis Terkait Variabel Kategorik

Analisis data kategorik memungkinkan peneliti untuk menguji hipotesis statistik mengenai distribusi atau hubungan antar kategori, seperti menggunakan uji chi-square untuk menguji independensi antara dua variabel kategorik.

1.2 Definisi dan Ruang Lingkup Analisis Data Kategorik Analisis data

kategorik adalah cabang statistik yang memfokuskan pada teknik untuk menganalisis data yang diklasifikasikan ke dalam kategori diskrit, baik dalam bentuk label tanpa urutan maupun label dengan urutan tertentu.

1.2.1 Jenis Data Kategorik

• Nominal

Data nominal terdiri atas kategori yang tidak memiliki urutan logis. Masing-masing kategori berdiri sendiri dan tidak menunjukkan tingkatan. Contoh: Warna mobil (merah, putih, hitam) Jenis kelamin (pria, wanita) Agama (Islam, Kristen, Hindu, Buddha)

• Ordinal

Data ordinal memiliki urutan atau tingkatan, namun perbedaan antara satu kategori dengan lainnya tidak dapat diukur secara kuantitatif. Contoh: Tingkat pendidikan (SD, SMP, SMA, Perguruan Tinggi) Skor kepuasan (sangat tidak puas, tidak puas, netral, puas, sangat puas)

1.2.2 Data Biner dan Multikategori

• Data Biner

Data biner merupakan data kategorik dengan dua kemungkinan kategori. Data ini sangat umum dalam penelitian, terutama ketika hasil yang diukur bersifat ya/tidak, benar/salah, atau hadir/tidak hadir. Contoh: Apakah seseorang membeli produk? (ya/tidak) Status merokok (perokok/bukan perokok) Hasil tes (positif/negatif)

• Data Multikategori

Data multikategori memiliki lebih dari dua kategori. Data ini bisa bersifat nominal maupun ordinal, tergantung pada apakah kategorinya memiliki urutan atau tidak. Contoh: Jenis pekerjaan (guru, dokter, pengacara, pengusaha) Tingkat stres (rendah, sedang, tinggi) Preferensi makanan (nasi, roti, mie, pasta)

1.3 Perbedaan dengan Data Kuantitatif

1.3.2 Sifat Data

o Data Kategorik: Data yang dikelompokkan ke dalam kategori atau label tanpa nilai numerik intrinsik.

Contohnya termasuk jenis kelamin (pria/wanita), status pernikahan (lajang/menikah), atau warna mata (biru/coklat). Data ini bersifat kualitatif dan tidak dapat diukur secara numerik.

o Data Kuantitatif: Data yang berupa angka dan mencerminkan besaran atau jumlah tertentu.

Contohnya termasuk tinggi badan (170 cm), berat badan (65 kg), atau jumlah anak (2). Data ini bersifat numerik dan dapat diukur serta dianalisis secara matematis.

1.3.2 Skala Pengukuran

o Data Kategorik: Biasanya diukur pada skala nominal atau ordinal. Skala nominal digunakan untuk kategori tanpa urutan, seperti jenis kelamin atau warna favorit. Skala ordinal digunakan untuk kategori dengan urutan tertentu, seperti tingkat pendidikan (SD, SMP, SMA) atau tingkat kepuasan (puas, netral, tidak puas).

o Data Kuantitatif: Diukur pada skala interval atau rasio. Skala interval memiliki jarak yang sama antara nilai, tetapi tanpa nol absolut, seperti suhu dalam Celsius. Skala rasio memiliki nol absolut, memungkinkan perbandingan rasio, seperti berat badan atau pendapatan.

1.3.3 Metode Analisis

o Data Kategorik: Analisis data ini melibatkan teknik seperti tabel kontingensi, uji chi-square, regresi logistik, dan model log-linear. Metode ini dirancang untuk menangani data yang tidak memenuhi asumsi distribusi normal dan varians homogen.

o Data Kuantitatif: Dianalisis menggunakan teknik seperti analisis regresi linier, uji t, ANOVA, dan korelasi Pearson. Metode ini mengasumsikan data berdistribusi normal dengan varians homogen.

1.3.4 Penyajian Data

o Data Kategorik: Sering disajikan dalam bentuk tabel frekuensi, diagram batang, atau diagram lingkaran untuk menunjukkan distribusi kategori.

oData Kuantitatif: Dapat disajikan menggunakan histogram, diagram garis, atau scatter plot untuk menggambarkan distribusi dan hubungan antar variabel numerik.

1.3.5Tujuan Pengukuran

o Data Kategorik: Bertujuan untuk mengklasifikasikan individu atau objek ke dalam kelompok atau kategori tertentu berdasarkan atribut atau karakteristik.

o Data Kuantitatif: Bertujuan untuk mengukur besaran atau jumlah suatu atribut, memungkinkan perhitungan statistik seperti rata-rata, median, dan standar deviasi.

1.4 Manfaat Analisis Data Kategorik di Berbagai Bidang

1.4.1 Kesehatan

Dalam bidang kesehatan, analisis data kategorik digunakan untuk memahami hubungan antara faktor risiko dan status kesehatan. Misalnya, menganalisis hubungan antara jenis kelamin dan prevalensi penyakit tertentu, atau antara status perokok dan kejadian penyakit jantung.

1.4.2 Pendidikan

Di sektor pendidikan, analisis ini membantu mengevaluasi faktor-faktor yang mempengaruhi prestasi siswa, seperti hubungan antara tingkat pendidikan orang tua dan prestasi akademik anak, atau antara metode pengajaran dan tingkat kelulusan siswa.

1.4.3 Pemasaran dan Bisnis

Perusahaan menggunakan analisis data kategorik untuk memahami preferensi konsumen, segmentasi pasar, dan efektivitas kampanye pemasaran. Misalnya, menganalisis hubungan antara jenis kelamin dan preferensi produk tertentu, atau antara usia dan loyalitas merek.

1.4.4 Sosial dan Ekonomi

Dalam penelitian sosial dan ekonomi, analisis data kategorik membantu mengidentifikasi pola dalam distribusi pendapatan, status pekerjaan, atau akses terhadap layanan publik. Misalnya, mengevaluasi hubungan antara status pekerjaan dan tingkat kepuasan hidup, atau antara lokasi geografis dan akses terhadap pendidikan.

1.4.5 Psikologi

Dalam psikologi, analisis ini digunakan untuk mempelajari hubungan antara variabel psikologis, seperti antara jenis kelamin dan preferensi terhadap jenis terapi tertentu, atau antara usia dan kecenderungan terhadap perilaku tertentu.

2 Metode dalam Analisis Data Kategorik

2.1 Distribusi Frekuensi dan Proporsi

Menghitung dan menganalisis distribusi frekuensi atau proporsi dari kategori dalam satu variabel untuk memahami pola dasar data.

2.2 Tabulasi Silang (Contingency Tables)

Menyusun tabel yang menunjukkan distribusi frekuensi dari dua atau lebih variabel kategorik untuk mengeksplorasi hubungan antar variabel tersebut.

2.3 Uji Chi-Square

Menggunakan uji statistik untuk menguji independensi atau kesesuaian antara variabel kategorik, seperti uji chi-square untuk tabel kontingensi.

2.4 Regresi Logistik Membangun model prediktif ketika

Variabel respons bersifat biner atau multinomial, memungkinkan analisis pengaruh variabel prediktor terhadap probabilitas kategori tertentu.

2.5 Analisis Asosiasi dan Ukuran Keterkaitan Mengukur kekuatan dan

Arah hubungan antara dua variabel kategorik. Ukuran seperti Odds Ratio, Risk Ratio, dan Koefisien Kontingensi digunakan untuk menilai derajat asosiasi.

2.6 Analisis Correspondence (Correspondence Analysis)

Analisis correspondence adalah teknik multivariat yang digunakan untuk menganalisis hubungan antara dua atau lebih variabel kategorik dengan menampilkan hubungan tersebut dalam bentuk visual.

2.7 Decision Tree

Decision tree adalah metode yang digunakan untuk memodelkan keputusan dan prediksi dengan cara membagi data ke dalam berbagai kategori berdasarkan serangkaian pertanyaan yang diteruskan secara berurutan. Setiap “cabang” pada pohon keputusan mewakili keputusan berdasarkan suatu atribut, dan “daun” mewakili hasil atau prediksi.

2.8 Random Forest

Random forest adalah metode ensemble yang menggabungkan beberapa pohon keputusan untuk meningkatkan akurasi prediksi. Ini dilakukan dengan membuat banyak pohon keputusan dengan menggunakan subset acak dari data dan fitur yang ada, kemudian menggabungkan hasilnya.

3 Distribusi Probabilitas dalam Data Kategorik

3.1 Distribusi Bernoulli

Distribusi Bernoulli adalah distribusi probabilitas diskrit untuk percobaan yang hanya memiliki dua hasil, yaitu sukses (1) dan gagal (0). Rumus: \[ P(X = 1) = p,\quad P(X = 0) = 1 - p \] Contoh: Dalam melempar sebuah koin, peluang munculnya sisi “Gambar” (misal sukses) adalah 0,5. Distribusi ini mengikuti distribusi Bernoulli.

3.2 Distribusi Multinomial

Distribusi multinomial adalah generalisasi dari distribusi binomial untuk percobaan yang memiliki lebih dari dua kemungkinan hasil. Rumus: \[ P(X_1 = x_1, X_2 = x_2, ..., X_k = x_k) = \frac{n!}{x_1!x_2!...x_k!} \cdot p_1^{x_1} \cdot p_2^{x_2} \cdot ... \cdot p_k^{x_k} \] Contoh: Dalam survei terhadap 100 orang tentang minuman favorit mereka (kopi, teh, jus), hasilnya: 50 orang memilih kopi, 30 orang memilih teh, 20 orang memilih jus.

3.3 Distribusi Binomial

Distribusi binomial menggambarkan jumlah keberhasilan dalam sejumlah percobaan Bernoulli. Rumus: \[ P(X = k) = \binom{n}{k} p^k (1 - p)^{n - k}\] Contoh: Sebuah produk elektronik dicek kualitasnya. Peluang sebuah produk cacat adalah 0,1. Dalam sampel 10 produk, peluang terdapat tepat 2 produk cacat dihitung menggunakan distribusi binomial.

3.4 Distribusi Poisson

Distribusi Poisson digunakan untuk menghitung jumlah kejadian dalam interval waktu atau ruang tertentu. Rumus: \[ P(X = k) = \frac{\lambda^k e^{-\lambda}}{k!}\] Contoh: Jumlah keluhan pelanggan di sebuah toko dalam satu jam rata-rata adalah 3. Probabilitas bahwa dalam satu jam ada 5 keluhan dapat dihitung dengan distribusi Poisson.

4 Desain Sampling dalam Analisis Data Kategori

Secara umum, desain sampling dibagi menjadi dua kategori utama, yaitu sampling prospektif dan sampling retrospektif. Masing-masing pendekatan memiliki ciri khas tersendiri dalam karakteristik sampel dan metode pengumpulan data, baik dalam studi observasional maupun eksperimental seperti eksperimen, studi kohort, dan studi kasus-kontrol.

4.1 Sampling Prospektif

Sampling prospektif adalah teknik pengumpulan data di mana subjek diidentifikasi terlebih dahulu, lalu diikuti selama periode waktu tertentu untuk mengamati hasil yang terjadi di masa depan. Pendekatan ini umum digunakan dalam studi kohort dan eksperimental. Beberapa metode sampling dalam kategori ini antara lain:

4.1.1 Eksperimen

Pada penelitian eksperimental, subjek dipilih secara acak dan dialokasikan ke dalam kelompok perlakuan maupun kelompok kontrol.

• Simple Random Sampling (SRS):

Setiap individu dalam populasi memiliki kesempatan yang sama untuk dipilih secara acak.

• Stratified Sampling:

Populasi dibagi ke dalam subkelompok (strata) berdasarkan karakteristik tertentu, lalu sampel diambil dari tiap strata.

• Cluster Sampling:

Populasi dikelompokkan menjadi beberapa cluster, kemudian satu atau lebih cluster dipilih secara acak sebagai sampel.

4.1.2 Studi Kohort Dalam

penelitian kohort, kelompok individu diikuti dari waktu ke waktu untuk mengamati hasil yang diukur.

• Sensus Sampling: Seluruh anggota populasiyang memenuhi kriteria diikutsertakan dalam penelitian.

• Systematic Sampling: Subjek dipilih secara sistematis berdasarkan interval tertentu dari daftar populasi. • Matched Sampling: Individu dalam kelompok kohort dicocokkan dengan individu serupa dari kelompok lain berdasarkan karakteristik tertentu.

4.2 Sampling Retrospektif

Sampling retrospektif merupakan metode pengumpulan data dari peristiwa yang telah terjadi dimasa lalu. Teknik ini sering digunakan dalam studi observasional untuk meneliti hubungan antara faktor risiko dan hasil yang sudah terjadi.

4.2.1 Studi Kasus-Kontrol

Dalam studi ini, individu dengan kondisi tertentu (kasus) dibandingkan dengan individu tanpa kondisi tersebut (kontrol).

• Purposive Sampling:

Subjek dipilih secara sengaja berdasarkan karakteristik tertentu yang relevan dengan penelitian.

• Snowball Sampling:

Responden yang telah direkrut diminta merekomendasikan individu lain yang sesuai kriteria penelitian.

•Incidence Density Sampling:

Kasus dan kontrol dipilih dari populasi yang sama dengan memperhitungkan waktu kejadian kasus.

4.2.2 Studi Kohort Retrospektif

Dalam studi ini, data historis digunakan untuk mengelompokkan individu berdasarkan paparan sebelumnya, kemudian hasilnya diamati.

• Quota Sampling: Subjek dipilih untuk memenuhi kuota tertentu berdasarkan kategori yang telah ditetapkan.

• Case-Base Sampling: Sampel diambil dari basis data populasi umum yang sudah tersedia.

• Systematic Sampling: Pemilihan subjek dilakukan berdasarkan interval tertentu dalam catatan yang ada.

5. Tabel Kontingensi 2x2

Tabel kontingensi 2×2 adalah alat statistik yang digunakan untuk menganalisis hubungan antara dua variabel kategorik, masing-masing dengan dua kategori. Tabel ini menyajikan frekuensi pengamatan dalam empat sel yang merepresentasikan kombinasi kategori dari kedua variabel tersebut. Contoh aplikasinya misalkan dalam sebuah studi kesehatan, peneliti ingin mengetahui hubungan antara kebiasaan merokok (Ya/Tidak) dan kejadian penyakit jantung (Ya/Tidak). Data yang diperoleh disusun dalam tabel kontingensi 2×2 untuk dianalisis apakah terdapat asosiasi antara merokok dan penyakit jantung.

Tabel kontingensi 2×2 digunakan untuk:

• Mengidentifikasi hubungan atau asosiasi antara dua variabel kategorik.

• Menghitung ukuran asosiasi seperti Odds Ratio (OR) dan Risk Ratio (RR).

• Melakukan uji independensi seperti uji Chi-Square untuk menentukan apakah terdapat hubungan yang signifikan antara kedua variabel.

Contoh Tabel:

	Merokok (B1)	Tidak Merokok (B2)	Total
Sakit Jantung (A1)	30	10	40
Tidak Sakit (A2)	20	40	60
Total	50	50	100

5.1 Distribusi Peluang dalam Tabel Kontingensi 2 × 2

5.1.1 Peluang Bersama (Joint Probability)

Definisi: Peluang dua kejadian terjadi secara bersamaan.

Rumus: \[ P(A \cap B) = \frac{n(A \cap B)}{n}\]

5.1.2 Peluang Marjinal (Marginal Probability)

Definisi: Peluang terjadinya satu kejadian tanpa memperhatikan kejadian lainnya.

Rumus:

\[ P(A) = \frac{n(A)}{n} \]

5.1.3 Peluang Bersyarat (Conditional Probability)

Definisi:

Peluang suatu kejadian terjadi dengan syarat kejadian lain telah terjadi.

Rumus:

\[ P(A|B) = \frac{P(A \cap B)}{P(B)} = \frac{n(A \cap B)}{n(B)} \]

Contoh Kasus:

	Merokok (B1)	Tidak Merokok (B2)	Total
Sakit Jantung (A1)	30	10	40
Tidak Sakit (A2)	20	40	60
Total	50	50	100

Peluang Bersama (Joint Probability)

\[ P(\text{Sakit dan Merokok}) = \frac{30}{100} = 0.3\] Peluang Marjinal \[ P(\text{Sakit}) = \frac{40}{100} = 0.40\] \[ P(\text{Merokok}) = \frac{50}{100} = 0.50\]

Peluang Bersyarat \[ P(\text{Sakit} | \text{Merokok}) = \frac{30}{50} = 0.60\]

5.2 Ukuran Asosiasi

Ukuran asosiasi dalam tabel kontingensi 2×2 digunakan untuk mengukur kekuatan dan arah hubungan antara dua variabel kategorik. Berikut adalah beberapa ukuran asosiasi yang umum digunakan:

5.2.1 Odds Ratio (OR)

Odds Ratio digunakan untuk membandingkan “peluang” (odds) suatu kejadian antara dua kelompok. Ini bukan probabilitas, tapi rasio antara jumlah sukses dan gagal dalam tiap kelompok. OR banyak dipakai dalam penelitian kasus-kontrol.

🔹 Interpretasi:

OR = 1 → Tidak ada hubungan antara eksposur dan outcome.

OR > 1 → Peluang outcome lebih besar pada kelompok eksposur.

OR < 1 → Peluang outcome lebih kecil pada kelompok eksposur.

odds_ratio <- (30 * 40) / (20 * 10)
odds_ratio

## [1] 6

Interpretasi: Peluang perokok untuk sakit jantung 6 kali lebih besar dibanding yang tidak merokok.

5.2.2 Risk Ratio (RR) / Relative Risk

Risk Ratio menghitung perbandingan probabilitas outcome antara dua kelompok. Digunakan dalam penelitian kohort (prospektif).

🔹 Interpretasi:

RR = 1 → Risiko sama pada kedua kelompok.

RR > 1 → Risiko lebih tinggi di kelompok ekspos.

RR < 1 → Risiko lebih rendah di kelompok ekspos.

a <- 30  # Merokok dan Sakit
b <- 20  # Merokok dan Tidak Sakit
c <- 10  # Tidak Merokok dan Sakit
d <- 40  # Tidak Merokok dan Tidak Sakit
risk_merokok <- a / (a + b)        # 30 / 50
risk_tidak_merokok <- c / (c + d)  # 10 / 50
RR <- risk_merokok / risk_tidak_merokok
RR

## [1] 3

artinya perokok 3 kali lebih berisiko terkena sakit jantung dibanding non-perokok.

5.2.3. Risk Difference (RD)

Risk Difference (juga dikenal sebagai Absolute Risk Reduction) adalah selisih probabilitas outcome antar dua kelompok.

🔹 Interpretasi:

RD = 0 → Tidak ada perbedaan risiko.

RD > 0 → Risiko lebih tinggi di kelompok eksposur.

RD < 0 → Risiko lebih tinggi di kelompok kontrol.

rd <- (30 / 50) - (10 / 50)
rd

## [1] 0.4

Ada perbedaan risiko sebesar 40% antara kelompok merokok dan tidak merokok.

5.2.4. Chi-Square Test (χ²)

Uji Chi-Square digunakan untuk menguji apakah ada hubungan signifikan antara dua variabel kategorik.

🔹 Interpretasi:

p-value < 0.05 → Ada asosiasi signifikan.

p-value ≥ 0.05 → Tidak signifikan; tidak ada bukti hubungan.

5.2.5. Fisher’s Exact Test

Alternatif dari Chi-Square Test jika jumlah data kecil, terutama jika frekuensi harapan (expected cell frequency) < 5. Menggunakan distribusi eksak daripada aproksimasi.

🔹 Interpretasi:

p-value < 0.05 → asosiasi signifikan

Sangat cocok untuk data medis atau studi dengan sampel kecil

Contoh: Dalam uji klinis kecil, Fisher Test bisa menunjukkan bahwa obat A lebih efektif dari obat B secara signifikan.

6 Inferensi Tabel Kontingensi 2 arah

Inferensi pada tabel kontingensi dua arah bertujuan untuk menganalisis hubungan antara dua variabel kategorik. Analisis ini melibatkan estimasi parameter dan pengujian hipotesis untuk menentukan apakah terdapat asosiasi yang signifikan antara variabel-variabel tersebut.

6.1 Estimasi

Estimasi bertujuan untuk memperkirakan parameter populasi berdasarkan data sampel.

6.1.1 Estimasi Titik

Estimasi titik digunakan untuk menentukan satu nilai spesifik sebagai perkiraan terbaik dari parameter populasi \[\hat{p} = \frac{x}{n}\] \[\hat{p} = estimasi\quad titik\quad proporsi\] x = adalah jumlah individu dalam kategori tertentu n = adalah total jumlah individu dalam sampel.

6.1.2 Estimasi Interval

Estimasi interval bertujuan untuk memberikan rentang nilai yang diyakini mengandung parameter populasi dengan tingkat kepercayaan tertentu. \[ \hat{p} \pm Z_{\alpha/2} \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}\]

6.2 Uji Hipotesis

6.2.1 Uji Proporsi

Untuk uji proporsi dalam uji hipotesis tabel kontingensi 2x2, biasanya digunakan uji Chi-Square atau uji proporsi untuk membandingkan dua proporsi yang terpisah, dan untuk menentukan apakah terdapat hubungan yang signifikan antara dua variabel kategori.

Proporsi untuk masing-masing kelompok dihitung sebagai:

• \(p_{AX} = \frac{a}{a+b}\)
• \(p_{AY} = \frac{b}{a+b}\)
• \(p_{BX} = \frac{c}{c+d}\)
• \(p_{BY} = \frac{d}{c+d}\)

Estimasi proporsi gabungan:

• \(p_X = \frac{a+c}{N}\)
• \(p_Y = \frac{b+d}{N}\)

Statistik uji untuk proporsi dua sampel

Rumus statistik uji:

\[ Z = \frac{\hat{p}_1 - \hat{p}_2} {\sqrt{\hat{p}(1 - \hat{p}) \left( \frac{1}{n_1} + \frac{1}{n_2} \right)}} \] Hipotesis Nol : Tidak ada perbedaan proporsi antara dua kelompok • Hipotesis Alternatif : Terdapat perbedaan proporsi antara dua kelompok

6.2.1.1 Contoh

Sebuah studi dilakukan untuk membandingkan efektivitas dua metode pengajaran. Hasilnya adalah sebagai berikut:

• Kelompok A (metode tradisional): 60 dari 100 siswa lulus
• Kelompok B (metode baru): 75 dari 120 siswa lulus

Apakah terdapat perbedaan proporsi kelulusan secara signifikan antara kedua metode pada tingkat signifikansi \(\alpha = 0.05\)?

Tabel Kontingensi

\[ \begin{array}{|c|c|c|c|} \hline & \text{Lulus} & \text{Tidak Lulus} & \text{Total} \\ \hline \text{Metode A} & 60 & 40 & 100 \\ \text{Metode B} & 75 & 45 & 120 \\ \hline \end{array} \]

Estimasi Proporsi Sampel

Proporsi siswa yang lulus:

• Metode A: \[ \hat{p}_1 = \frac{60}{100} = 0.6 \]
• Metode B: \[ \hat{p}_2 = \frac{75}{120} = 0.625 \]

Proporsi Gabungan

Gabungan kelulusan dari dua kelompok: \[ \hat{p} = \frac{60 + 75}{100 + 120} = \frac{135}{220} \approx 0.6136 \] Statistik Uji Z

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

Substitusi nilai: \[ Z = \frac{0.6 - 0.625} {\sqrt{0.6136(1 - 0.6136) \left( \frac{1}{100} + \frac{1}{120} \right)}} \]

Hitung denominator: \[ \sqrt{0.6136 \times 0.3864 \left( \frac{1}{100} + \frac{1}{120} \right)} = \sqrt{0.2371 \times 0.01833} \approx \sqrt{0.00435} \approx 0.0659 \]

Sehingga: \[ Z = \frac{-0.025}{0.0659} \approx -0.3795 \]

Keputusan

Dengan \(\alpha = 0.05\) (dua sisi), nilai kritis \(Z_{0.025} = \pm1.96\).

Karena \(|Z| = 0.3795 < 1.96\), maka:

Gagal menolak H₀. Tidak ada cukup bukti untuk menyatakan bahwa kedua metode menghasilkan proporsi kelulusan yang berbeda secara signifikan.

6.2.2 Uji Asosiasi

Uji asosiasi dalam tabel kontingensi \(2 \times 2\) bertujuan untuk mengukur hubungan antara dua variabel kategori. Tiga ukuran utama dalam uji asosiasi adalah:

1. Risk Difference (RD) – Mengukur selisih risiko absolut antara dua kelompok.

2. Relative Risk (RR) – Mengukur perbandingan risiko antara dua kelompok.

3. Odds Ratio (OR) – Mengukur perbandingan odds antara dua kelompok.

---

Hipotesis Uji dalam Tabel Kontingensi \(2 \times 2\)

Untuk setiap uji asosiasi, hipotesis yang diuji adalah:

• Hipotesis Nol (H\(_0\)): Tidak ada asosiasi antara dua variabel.
• Hipotesis Alternatif (H\(_1\)): Terdapat asosiasi antara dua variabel.

---

6.2.2.1 Risk Difference (RD)

Risk Difference mengukur perbedaan absolut dalam probabilitas kejadian antara dua kelompok:

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

Standard Error untuk RD:

\[ SE(RD) = \sqrt{ \frac{\hat{p}_1(1-\hat{p}_1)}{n_{1\cdot}} + \frac{\hat{p}_2(1-\hat{p}_2)}{n_{2\cdot}} } \]

Statistik uji Z:

\[ Z_{RD} = \frac{RD}{SE(RD)} \]

---

6.2.2.2 Relative Risk (RR)

Relative Risk membandingkan kemungkinan kejadian antara dua kelompok:

\[ RR = \frac{\hat{p}_1}{\hat{p}_2} \]

Standard Error untuk log(RR):

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

Statistik uji Z:

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

---

6.2.2.3 Odds Ratio (OR)

Odds Ratio membandingkan peluang kejadian antara dua kelompok:

\[ OR = \frac{n_{11} \times n_{22}}{n_{12} \times 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}} } \]

Statistik uji Z:

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

6.2.3 Uji Independensi

Uji independensi digunakan untuk menentukan apakah ada hubungan statistik antara dua variabel kategorikal.

6.2.3.1 Uji Chi-Square

Uji Chi-Square digunakan untuk menguji apakah ada hubungan antara dua variabel kategorikal.

Rumus Chi-Square:

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

Dimana:

• \(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.

---

6.2.3.2 Partisi Chi-Square

Partisi Chi-Square digunakan untuk mengidentifikasi kategori mana dalam tabel kontingensi yang bertanggung jawab atas hubungan yang signifikan. Jika uji Chi-Square pada tabel kontingensi I × J signifikan, maka partisi Chi-Square memungkinkan kita untuk menguraikan efek hubungan dalam subkelompok yang lebih kecil.

Menurut Lancaster dan Irwin, statistik uji chi-square dapat dipecah menjadi komponen-komponen sebanyak derajat bebas dari tabel kontingensi, memungkinkan identifikasi kategori yang berkontribusi pada hubungan yang signifikan.

---

6.2.3.3 Contoh Kasus dan Perhitungan

Kasus: Peneliti ingin mengetahui apakah penggunaan media sosial berhubungan dengan tingkat stres pada remaja.

Data:

	Stres Tinggi	Stres Rendah	Total
Pengguna Intensif	30	20	50
Pengguna Ringan	10	40	50
Total	40	60	100

Langkah 1: Hitung nilai \(E\) (ekspektasi) untuk masing-masing sel:

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

Langkah 2: Hitung Chi-Square:

\[ \chi^2 = \frac{(30-20)^2}{20} + \frac{(20-30)^2}{30} + \frac{(10-20)^2}{20} + \frac{(40-30)^2}{30} \] \[ = \frac{100}{20} + \frac{100}{30} + \frac{100}{20} + \frac{100}{30} \] \[ = 5 + 3.33 + 5 + 3.33 \] \[ = 16.66\]

Langkah 3: Derajat bebas (df):

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

Langkah 4: Interpretasi: - Dengan \(\alpha = 0.05\), nilai kritis \(\chi^2\) dari tabel adalah 3.841. - Karena 16.66 > 3.841, maka tolak H₀. - Kesimpulan: Ada hubungan signifikan antara penggunaan media sosial dengan tingkat stres remaja.

6.2.3.4 Uji Likelihood Ratio (G²)

Uji Likelihood Ratio (\(G^2\)) 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 pada sel ke-\((i,j)\) dalam tabel kontingensi
  
• \(\hat{\mu}_{ij} = n \cdot p_{i+} \cdot p_{+j}\) = frekuensi ekspektasi, dihitung dari baris dan kolom marjinal

Statistik \(G^2\) mengikuti distribusi chi-square dengan derajat bebas: \[ df = (I - 1)(J - 1) \]

Tolak \(H_0\) jika: \[ G^2 \geq \chi^2_{(1 - \alpha), (I-1)(J-1)} \]

Contoh Kasus

Peneliti ingin mengetahui apakah terdapat hubungan antara jenis kelamin dan preferensi minuman. Data yang diperoleh sebagai berikut:

\[ \begin{array}{|c|c|c|c|} \hline & \text{Kopi} & \text{Teh} & \text{Total} \\ \hline \text{Laki-laki} & 40 & 30 & 70 \\ \text{Perempuan} & 20 & 50 & 70 \\ \hline \text{Total} & 60 & 80 & 140 \\ \hline \end{array} \]

Hitung Frekuensi Ekspektasi

Menggunakan rumus: \[ \hat{\mu}_{ij} = \frac{(\text{total baris } i) \times (\text{total kolom } j)}{n} \]

Ekspektasi masing-masing sel:

• \(\hat{\mu}_{11} = \frac{70 \times 60}{140} = 30\)
• \(\hat{\mu}_{12} = \frac{70 \times 80}{140} = 40\)
• \(\hat{\mu}_{21} = \frac{70 \times 60}{140} = 30\)
• \(\hat{\mu}_{22} = \frac{70 \times 80}{140} = 40\)

Hitung Statistik G² \[ G^2 = 2 \left[ 40 \ln\left(\frac{40}{30}\right) + 30 \ln\left(\frac{30}{40}\right) + 20 \ln\left(\frac{20}{30}\right) + 50 \ln\left(\frac{50}{40}\right) \right] \]

Hitung nilai logaritma:

• \(40 \ln(1.333) \approx 40 \times 0.2877 = 11.51\)
• \(30 \ln(0.75) \approx 30 \times (-0.2877) = -8.63\)
• \(20 \ln(0.6667) \approx 20 \times (-0.4055) = -8.11\)
• \(50 \ln(1.25) \approx 50 \times 0.2231 = 11.16\)

Sehingga:

\[ G^2 = 2 \times (11.51 - 8.63 - 8.11 + 11.16) = 2 \times 5.93 = 11.86 \]

Keputusan

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

Dengan \(\alpha = 0.05\), nilai kritis \(\chi^2_{0.95, 1} = 3.841\)

Karena \(G^2 = 11.86 > 3.841\), maka:

Tolak H₀. Terdapat hubungan yang signifikan antara jenis kelamin dan preferensi minuman.

6.2.3.5 Uji Exact Fisher

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

6.2.3.6 Uji Distribusi Hipergeomterik

digunakan untuk menghitung peluang terjadinya sejumlah keberhasilan dalam pengambilan sampel tanpa pengembalian dari populasi terbatas.

Rumus distribusi hipergeometrik:

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

Keterangan: - \(N\) = jumlah total populasi - \(K\) = jumlah item sukses dalam populasi - \(n\) = jumlah item yang diambil - \(x\) = jumlah sukses yang diinginkan

6.2.3.7 Contoh Kasus

Kasus:
Suatu pabrik memiliki 15 produk, terdiri dari 5 produk cacat dan 10 produk bagus.
Jika 4 produk diambil secara acak tanpa pengembalian, berapa peluang tepat 2 produk yang diambil adalah produk cacat?

Parameter:

• \(N = 15\) (total produk)
• \(K = 5\) (produk cacat)
• \(n = 4\) (produk diambil)
• \(x = 2\) (produk cacat yang diinginkan)

Perhitungan :

\[ P(X = 2) = \frac{{\binom{5}{2} \binom{10}{2}}}{\binom{15}{4}} \]

Hitung kombinasi: - \(\binom{5}{2} = 10\) - \(\binom{10}{2} = 45\) - \(\binom{15}{4} = 1365\)

Sehingga:

\[ P(X = 2) = \frac{10 \times 45}{1365} = \frac{450}{1365} \approx 0.3297 \]

Jadi, peluangnya sekitar 32.97%.

6.2.3.8 Uji Fisher’s Exact

Uji digunakan untuk menguji hubungan antara dua variabel kategorik dalam tabel kontingensi \(2 \times 2\), terutama jika ukuran sampel kecil sehingga asumsi chi-square tidak terpenuhi.

Probabilitas tabel \(2 \times 2\) tertentu dihitung menggunakan rumus:

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

dengan:

a,b,c,d = nilai dalam sel tabel

n = total keseluruhan

6.2.3.9 Contoh kasus

Sebuah studi ingin mengetahui apakah ada hubungan antara jenis kelamin (laki-laki/perempuan) dan hasil tes (lulus/gagal).

Data yang diperoleh:

\[ \begin{array}{|c|c|c|} \hline \textbf{} & \textbf{Lulus} & \textbf{Gagal} \\ \hline \textbf{Laki-laki} & 8 & 2 \\ \hline \textbf{Perempuan} & 1 & 5 \\ \hline \end{array} \] Sehingga: \[ a = 8, \quad b = 2, \quad c = 1, \quad d = 5, \quad n = 16\]

Gunakan rumus:

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

Substitusi:

\[ P = \frac{(8+2)!(1+5)!(8+1)!(2+5)!}{8!2!1!5!16!} \]

Hitung faktorial:

\[ P = \frac{10! \times 6! \times 9! \times 7!}{8! \times 2! \times 1! \times 5! \times 16!} \]

Maka hasil \(P\) diperoleh sekitar \(0.0107\).

Kesimpulan

Karena nilai \(P\) kecil (\(P < 0.05\)), maka terdapat hubungan yang signifikan antara jenis kelamin dan hasil tes.

6.3 Analisis Residual dalam Tabel Kontingensi

Residual dalam tabel kontingensi adalah selisih antara frekuensi observasi dengan ekspektasi, biasanya dibagi dengan deviasi standarnya agar bisa dibandingkan antar sel. Residual membantu mengidentifikasi pola atau hubungan spesifik dalam data yang mungkin tidak terlihat dari hasil uji global seperti \(\chi^2\).

6.3.1. Residual Standar (Pearson Residual)

\[ r_{ij} = \frac{n_{ij} - \hat{\mu}_{ij}}{\sqrt{\hat{\mu}_{ij}}} \]

• \(n_{ij}\): frekuensi observasi pada sel ke-\((i,j)\)
• \(\hat{\mu}_{ij}\): frekuensi harapan pada sel ke-\((i,j)\)

Residual ini menunjukkan berapa banyak standar deviasi nilai observasi dari nilai yang diharapkan.

6.3.2 Menggunakan Residual Untuk Deteksi Outlier

Pearson Residual:

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

Jika nilai \(|e_{ij}| > 2\), maka sel tersebut dianggap sebagai indikasi adanya outlier.

• Standardized Residual (Adjusted Residual):

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

Jika \(|r_{ij}| > 3\), maka sel tersebut dianggap sebagai outlier signifikan.

7 Tabel Kontingensi 3 Arah

Tabel kontingensi tiga arah adalah perpanjangan dari tabel kontingensi dua arah yang digunakan untuk menganalisis hubungan antara tiga variabel kategori secara simultan.

7.1 Tabel Parsial dan Marginal

Tabel yg meengabaikan nilai Z pada 3 arah X,Y,Z Struktur Tabel Kontingensi Tiga Arah:

\[\begin{array}{|c|c|c|c|c|c|c|} \hline & Y = 1 & Y = 2 & \cdots & Y = J & Y = . \\ \hline Z = 1\quad X=1 & n_{111} & n_{121} & \cdots & n_{1J1} & n_{1+1} \\ X=2 & n_{211} & n_{221} & \cdots & n_{2J1} & n_{2+1} \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ X=i & n_{i11} & n_{i21} & \cdots & n_{iJ1} & n_{i+1} \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ X=I & n_{I11} & n_{I21} & \cdots & n_{IJ1} & n_{I+1} \\ \hline & n_{+11} & n_{+21} & \cdots & n_{+J1} & n_{++1} \\ \hline Z=k\quad X=1 & n_{11k} & n_{12k} & \cdots & n_{1Jk} & n_{1+k} \\ X=2 & n_{21k} & n_{22k} & \cdots & n_{2Jk} & n_{2+k} \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ X=i & n_{i1k} & n_{i2k} & \cdots & n_{iJk} & n_{i+k} \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ X=I & n_{I1k} & n_{I2k} & \cdots & n_{IJk} & n_{I+k} \\ \hline & n_{+1k} & n_{+2k} & \cdots & n_{+Jk} & n_{++k} \\ \hline Z=K\quad X=1 & n_{11K} & n_{12K} & \cdots & n_{1JK} & n_{1+K} \\ X=2 & n_{21K} & n_{22K} & \cdots & n_{2JK} & n_{2+K} \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ X=i & n_{i1K} & n_{i2K} & \cdots & n_{iJK} & n_{i+K} \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ X=I & n_{I1K} & n_{I2K} & \cdots & n_{IJK} & n_{I+K} \\ \hline & n_{+1K} & n_{+2K} & \cdots & n_{+JK} & n_{++K} \\ \hline & n_{1++} & n_{2++} & \cdots & n_{J++} & n_{+++} \\ \hline \end{array}\]

7.2 Contoh

# Membuat tabel kontingensi
data <- matrix(c(45, 15, 60,
                 35, 10, 45,
                 25, 50, 75,
                 30, 40, 70,
                 15, 85, 100,
                 10, 90, 100), 
               nrow=6, byrow=TRUE)

colnames(data) <- c("Baik", "Buruk", "Total")
rownames(data) <- c("Rendah - Instagram", "Rendah - TikTok",
                    "Sedang - Instagram", "Sedang - TikTok",
                    "Tinggi - Instagram", "Tinggi - TikTok")
data

##                    Baik Buruk Total
## Rendah - Instagram   45    15    60
## Rendah - TikTok      35    10    45
## Sedang - Instagram   25    50    75
## Sedang - TikTok      30    40    70
## Tinggi - Instagram   15    85   100
## Tinggi - TikTok      10    90   100

Tabel Parsial

Hubungan antara dua variabel kategori dalam tabel kontingensi tiga arah dengan mempertahankan satu variabel sebagai kontrol

library(gmodels)

## Warning: package 'gmodels' was built under R version 4.3.3

# Membuat data tabel kontingensi 3 arah
data_medsos <- array(c(
  45, 15,  # Rendah - Instagram
  35, 10,  # Rendah - TikTok
  25, 50,  # Sedang - Instagram
  30, 40,  # Sedang - TikTok
  15, 85,  # Tinggi - Instagram
  10, 90   # Tinggi - TikTok
), dim = c(2,2,3), 
dimnames = list(
  "Kesehatan Mental" = c("Baik", "Buruk"),
  "Jenis Medsos" = c("Instagram", "TikTok"),
  "Durasi" = c("Rendah", "Sedang", "Tinggi")
))

# Menampilkan tabel kontingensi parsial untuk masing-masing durasi
for (d in dimnames(data_medsos)$Durasi) {
  cat("\nTabel Kontingensi untuk Durasi:", d, "\n")
  print(data_medsos[, , d])
}

## 
## Tabel Kontingensi untuk Durasi: Rendah 
##                 Jenis Medsos
## Kesehatan Mental Instagram TikTok
##            Baik         45     35
##            Buruk        15     10
## 
## Tabel Kontingensi untuk Durasi: Sedang 
##                 Jenis Medsos
## Kesehatan Mental Instagram TikTok
##            Baik         25     30
##            Buruk        50     40
## 
## Tabel Kontingensi untuk Durasi: Tinggi 
##                 Jenis Medsos
## Kesehatan Mental Instagram TikTok
##            Baik         15     10
##            Buruk        85     90

Tabel Marginal

Jumlah total observasi untuk setiap variabel dengan mengabaikan variabel lainnya dalam tabel kontingensi tiga arah

# Definisikan data dalam bentuk array (3-way contingency table)
data_medsos <- array(
  c(45, 15, 35, 10,  # Durasi = Rendah
    25, 50, 30, 40,  # Durasi = Sedang
    15, 85, 10, 90), # Durasi = Tinggi
  dim = c(2, 2, 3), 
  dimnames = list(
    "Kesehatan Mental" = c("Baik", "Buruk"),
    "Media Sosial" = c("Instagram", "TikTok"),
    "Durasi" = c("Rendah", "Sedang", "Tinggi")
  )
)

# 1️⃣ Tabel Marginal Kesehatan Mental (Mengabaikan Media Sosial & Durasi)
marginal_kesehatan <- apply(data_medsos, c(1), sum)
print("Tabel Marginal Kesehatan Mental:")

## [1] "Tabel Marginal Kesehatan Mental:"

print(marginal_kesehatan)

##  Baik Buruk 
##   160   290

# 2️⃣ Tabel Marginal Media Sosial + Kesehatan Mental (Mengabaikan Durasi)
marginal_medsos_kesehatan <- apply(data_medsos, c(1,2), sum)
print("Tabel Marginal Media Sosial + Kesehatan Mental:")

## [1] "Tabel Marginal Media Sosial + Kesehatan Mental:"

print(marginal_medsos_kesehatan)

##                 Media Sosial
## Kesehatan Mental Instagram TikTok
##            Baik         85     75
##            Buruk       150    140

# 4️⃣ Tabel Marginal Durasi Penggunaan + Kesehatan Mental (Mengabaikan Media Sosial)
marginal_durasi_kesehatan <- apply(data_medsos, c(1,3), sum)
print("Tabel Marginal Durasi Penggunaan + Kesehatan Mental:")

## [1] "Tabel Marginal Durasi Penggunaan + Kesehatan Mental:"

print(marginal_durasi_kesehatan)

##                 Durasi
## Kesehatan Mental Rendah Sedang Tinggi
##            Baik      80     55     25
##            Buruk     25     90    175

7.3 Distribusi Peluang

7.3.1 Menghitung Join Probability

# Hitung total responden
total_responden <- sum(data_medsos)

# Hitung Joint Probability dengan membagi setiap sel dengan total responden
joint_prob <- data_medsos / total_responden

# Cetak hasil
print("Joint Probability Table:")

## [1] "Joint Probability Table:"

print(joint_prob)

## , , Durasi = Rendah
## 
##                 Media Sosial
## Kesehatan Mental  Instagram     TikTok
##            Baik  0.10000000 0.07777778
##            Buruk 0.03333333 0.02222222
## 
## , , Durasi = Sedang
## 
##                 Media Sosial
## Kesehatan Mental  Instagram     TikTok
##            Baik  0.05555556 0.06666667
##            Buruk 0.11111111 0.08888889
## 
## , , Durasi = Tinggi
## 
##                 Media Sosial
## Kesehatan Mental  Instagram     TikTok
##            Baik  0.03333333 0.02222222
##            Buruk 0.18888889 0.20000000

7.3.2 Peluang Bersyarat

# Hitung total responden
total_responden <- sum(data_medsos)

# 1️⃣ P(Kesehatan Mental | Media Sosial, Durasi)
# Menjumlahkan berdasarkan Media Sosial & Durasi
sum_medsos_durasi <- apply(data_medsos, c(2,3), sum)

# Ubah agar dimensi sesuai untuk pembagian
sum_medsos_durasi <- array(sum_medsos_durasi, dim = c(2,2,3))  # Pastikan dimensi sama

# Hitung peluang bersyarat
cond_prob_kesehatan_given_medsos_durasi <- data_medsos / sum_medsos_durasi

# Cetak hasil
print("P(Kesehatan Mental | Media Sosial, Durasi):")

## [1] "P(Kesehatan Mental | Media Sosial, Durasi):"

print(cond_prob_kesehatan_given_medsos_durasi)

## , , Durasi = Rendah
## 
##                 Media Sosial
## Kesehatan Mental Instagram    TikTok
##            Baik  0.7500000 0.4666667
##            Buruk 0.3333333 0.1428571
## 
## , , Durasi = Sedang
## 
##                 Media Sosial
## Kesehatan Mental Instagram    TikTok
##            Baik       0.25 0.5000000
##            Buruk      0.50 0.8888889
## 
## , , Durasi = Tinggi
## 
##                 Media Sosial
## Kesehatan Mental Instagram TikTok
##            Baik   0.200000    0.1
##            Buruk  1.214286    0.9

# 2️⃣ P(Media Sosial | Kesehatan Mental, Durasi)
# Menjumlahkan berdasarkan Kesehatan Mental & Durasi
sum_kesehatan_durasi <- apply(data_medsos, c(1,3), sum)

# Ubah agar dimensi sesuai untuk pembagian
sum_kesehatan_durasi <- array(sum_kesehatan_durasi, dim = c(2,2,3))  # Sesuaikan dimensi

# Hitung peluang bersyarat
cond_prob_medsos_given_kesehatan_durasi <- data_medsos / sum_kesehatan_durasi

# Cetak hasil
print("P(Media Sosial | Kesehatan Mental, Durasi):")

## [1] "P(Media Sosial | Kesehatan Mental, Durasi):"

print(cond_prob_medsos_given_kesehatan_durasi)

## , , Durasi = Rendah
## 
##                 Media Sosial
## Kesehatan Mental Instagram    TikTok
##            Baik     0.5625 0.6363636
##            Buruk    0.6000 0.1111111
## 
## , , Durasi = Sedang
## 
##                 Media Sosial
## Kesehatan Mental Instagram TikTok
##            Baik  1.0000000  0.375
##            Buruk 0.2857143  1.600
## 
## , , Durasi = Tinggi
## 
##                 Media Sosial
## Kesehatan Mental Instagram    TikTok
##            Baik  0.2727273 0.4000000
##            Buruk 0.9444444 0.5142857

# 3️⃣ P(Durasi | Kesehatan Mental, Media Sosial)
# Menjumlahkan berdasarkan Kesehatan Mental & Media Sosial
sum_kesehatan_medsos <- apply(data_medsos, c(1,2), sum)

# Ubah agar dimensi sesuai untuk pembagian
sum_kesehatan_medsos <- array(sum_kesehatan_medsos, dim = c(2,2,3))  # Sesuaikan dimensi

# Hitung peluang bersyarat
cond_prob_durasi_given_kesehatan_medsos <- data_medsos / sum_kesehatan_medsos

# Cetak hasil
print("P(Durasi | Kesehatan Mental, Media Sosial):")

## [1] "P(Durasi | Kesehatan Mental, Media Sosial):"

print(cond_prob_durasi_given_kesehatan_medsos)

## , , Durasi = Rendah
## 
##                 Media Sosial
## Kesehatan Mental Instagram     TikTok
##            Baik  0.5294118 0.46666667
##            Buruk 0.1000000 0.07142857
## 
## , , Durasi = Sedang
## 
##                 Media Sosial
## Kesehatan Mental Instagram    TikTok
##            Baik  0.2941176 0.4000000
##            Buruk 0.3333333 0.2857143
## 
## , , Durasi = Tinggi
## 
##                 Media Sosial
## Kesehatan Mental Instagram    TikTok
##            Baik  0.1764706 0.1333333
##            Buruk 0.5666667 0.6428571

7.3.2 Perhitungan Ukuran Asosiasi

# Membuat tabel kontingensi untuk Instagram vs TikTok
tabel_medsos <- matrix(c(50, 25, 40, 30), nrow=2, byrow=TRUE)
colnames(tabel_medsos) <- c("Buruk", "Baik")
rownames(tabel_medsos) <- c("Instagram", "TikTok")


tabel_medsos

##           Buruk Baik
## Instagram    50   25
## TikTok       40   30

# Menentukan probabilitas kejadian
P_instagram <- tabel_medsos[1,1] / sum(tabel_medsos[1,])
P_tiktok <- tabel_medsos[2,1] / sum(tabel_medsos[2,])

# Menghitung Risk Difference
RD <- P_instagram - P_tiktok
RD

## [1] 0.0952381

# Menghitung Risk Ratio
RR <- P_instagram / P_tiktok
RR

## [1] 1.166667

# Odds Ratio
# Menentukan nilai dari tabel kontingensi
# Menentukan nilai dari tabel kontingensi
a <- tabel_medsos[1,1]  # Buruk, Instagram
b <- tabel_medsos[1,2]  # Baik, Instagram
c <- tabel_medsos[2,1]  # Buruk, TikTok
d <- tabel_medsos[2,2]  # Baik, TikTok

# Menghitung Odds Ratio secara manual
OR_manual <- (a * d) / (b * c)

# Menampilkan hasil
OR_manual

## [1] 1.5

INTERPRETASI :

1. Risk Difference (RD) = 0.0953 (9,53%)
   • Risiko kesehatan mental buruk pada pengguna Instagram 9,53% lebih tinggi dibandingkan pengguna TikTok.
2. Risk Ratio (RR) = 1.17
   • Pengguna Instagram memiliki 1,17 kali lipat risiko lebih tinggi mengalami kesehatan mental buruk dibandingkan pengguna TikTok.
3. Odds Ratio (OR) = 1.5
   • Peluang mengalami kesehatan mental buruk 1,5 kali lebih besar bagi pengguna Instagram dibandingkan pengguna TikTok.

7.3.2 UJI INDEPENDENSI

Independensi Bersyarat: Dua variabel, X dan Y , dikatakan independen bersyarat terhadap variabel ketiga, Z, jika rasio odds mereka dalam setiap strata Z sama dengan 1.

Metode Cochran-Mantel-Haenszel Tujuan Uji CMH Uji Cochran-Mantel-Haenszel (CMH) digunakan untuk menguji hubungan antara dua variabel kategori dengan mempertimbangkan efek dari variabel perancu.

# Memuat pustaka yang diperlukan
library(vcdExtra)

## Warning: package 'vcdExtra' was built under R version 4.3.3

## Loading required package: vcd

## Warning: package 'vcd' was built under R version 4.3.3

## Loading required package: grid

## Loading required package: gnm

## Warning: package 'gnm' was built under R version 4.3.3

# Membuat array 3 dimensi untuk data kontingensi
data_cmh <- array(c(
  45, 15,  35, 10,  # Rendah
  25, 50,  30, 40,  # Sedang
  15, 85,  10, 90   # Tinggi
), dim = c(2,2,3), 
dimnames = list(
  "Kesehatan Mental" = c("Baik", "Buruk"),
  "Media Sosial" = c("Instagram", "TikTok"),
  "Durasi" = c("Rendah", "Sedang", "Tinggi")
))

# Menampilkan tabel kontingensi
print(data_cmh)

## , , Durasi = Rendah
## 
##                 Media Sosial
## Kesehatan Mental Instagram TikTok
##            Baik         45     35
##            Buruk        15     10
## 
## , , Durasi = Sedang
## 
##                 Media Sosial
## Kesehatan Mental Instagram TikTok
##            Baik         25     30
##            Buruk        50     40
## 
## , , Durasi = Tinggi
## 
##                 Media Sosial
## Kesehatan Mental Instagram TikTok
##            Baik         15     10
##            Buruk        85     90

# Uji Cochran-Mantel-Haenszel (CMH)
cmh_test <- mantelhaen.test(data_cmh)

# Menampilkan hasil uji
print(cmh_test)

## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  data_cmh
## Mantel-Haenszel X-squared = 0.071932, df = 1, p-value = 0.7885
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.5819488 1.4391388
## sample estimates:
## common odds ratio 
##          0.915153

INTERPRETASI:

• Nilai common OR = 0.915 berarti bahwa peluang memiliki kesehatan mental baik bagi pengguna TikTok dibandingkan Instagram hampir sama, tetapi sedikit lebih rendah untuk TikTok. Penggunaan Instagram memiliki sedikit keunggulan dalam kesehatan mental dibandingkan TikTok, tetapi perbedaannya sangat kecil.

• p-value = 0.7885 (> 0.05) → Tidak ada hubungan signifikan antara media sosial yang digunakan (TikTok vs. Instagram) dengan kesehatan mental, setelah mengontrol durasi penggunaan.

-Interval Kepercayaan 95% odds ratio 0.582 ≤ � ≤ 1.439 Karena Confidence Interval mencakup 1, maka kita tidak dapat menyimpulkan adanya pengaruh yang signifikan dari jenis media sosial terhadap kesehatan mental.

7.3.5 Uji Odd Ratio Bersama

Jika nilai-nilai odds ratio bersyarat relatif sama (tidak berbeda secara ekstrem) dan memiliki arah yang sama, maka kita dapat menentukan sebuah nilai tunggal untuk odds ratio yang disebut odds ratio bersama.

• Odds ratio bersama ditaksir oleh statistik Mantel-Haenszel.

Perhitungan Odds Ratio Bersama

# Data dari tabel kontingensi 3 arah
data_medsos <- data.frame(
  Media_Sosial = rep(c("TikTok", "Instagram"), each = 2),
  Jenis_Kelamin = rep(c("Laki-laki", "Perempuan"), 2),
  A = c(40, 30, 50, 45),  # Kesehatan Baik
  B = c(60, 70, 50, 55),  # Kesehatan Buruk
  C = c(50, 45, 40, 30),  # Kontrol (baik)
  D = c(50, 55, 60, 70)   # Kontrol (buruk)
)

# Total N per strata
data_medsos$N <- data_medsos$A + data_medsos$B + data_medsos$C + data_medsos$D

# Hitung Common Odds Ratio (OR_MH)
OR_MH <- sum((data_medsos$A * data_medsos$D) / data_medsos$N) / 
         sum((data_medsos$B * data_medsos$C) / data_medsos$N)

# Log Odds Ratio
log_OR_MH <- log(OR_MH)

# Standard Error (SE) dari log(OR)
SE_log_OR <- sqrt(sum(1 / data_medsos$A + 1 / data_medsos$B + 
                      1 / data_medsos$C + 1 / data_medsos$D))

# Interval Kepercayaan 95%
CI_lower <- exp(log_OR_MH - 1.96 * SE_log_OR)
CI_upper <- exp(log_OR_MH + 1.96 * SE_log_OR)

# Cetak hasil
cat("Common Odds Ratio (OR_MH):", round(OR_MH, 2), "\n")

## Common Odds Ratio (OR_MH): 1

cat("Log(OR_MH):", round(log_OR_MH, 3), "\n")

## Log(OR_MH): 0

cat("Standard Error (SE):", round(SE_log_OR, 3), "\n")

## Standard Error (SE): 0.583

cat("95% Confidence Interval: [", round(CI_lower, 2), ",", round(CI_upper, 2), "]\n")

## 95% Confidence Interval: [ 0.32 , 3.13 ]

INTERPRETASI :

• OR_MH = 1 menunjukkan bahwa peluang memiliki kesehatan mental baik antara pengguna TikTok dan Instagram adalah sama.

• Log(OR) = 0 karena log dari 1 adalah 0, yang mengonfirmasi bahwa tidak ada hubungan antara variabel.

• SE = 0.583 menunjukkan adanya ketidakpastian dalam estimasi, tetapi tidak cukup besar untuk menyimpulkan efek yang kuat.

• -Interval Kepercayaan 95% odds ratio 0.32 ≤ � ≤ 3.13 Karena Confidence Interval mencakup 1, maka tidak dapat menyimpulkan adanya pengaruh yang signifikan dari jenis media sosial terhadap kesehatan mental.

7.3.6 Uji Homogentias Odds Ratio

Asosiasi homogen terjadi jika odds ratio pada setiap tabel parsial bernilai sama Jika odds ratio konstan di semua strata variabel kontrol (Z), maka tidak ada interaksi antara variabel x dan y dalam pengaruhnya terhadap Z.

perhitungan

# Data: 3-way contingency table (Media Sosial, Kesehatan Mental, Durasi)
data_medsos <- array(
  c(45, 15, 35, 10,  # Durasi = Rendah
    25, 50, 30, 40,  # Durasi = Sedang
    15, 85, 10, 90), # Durasi = Tinggi
  dim = c(2, 2, 3), 
  dimnames = list(
    "Kesehatan Mental" = c("Baik", "Buruk"),
    "Media Sosial" = c("Instagram", "TikTok"),
    "Durasi" = c("Rendah", "Sedang", "Tinggi")
  )
)
data_medsos

## , , Durasi = Rendah
## 
##                 Media Sosial
## Kesehatan Mental Instagram TikTok
##            Baik         45     35
##            Buruk        15     10
## 
## , , Durasi = Sedang
## 
##                 Media Sosial
## Kesehatan Mental Instagram TikTok
##            Baik         25     30
##            Buruk        50     40
## 
## , , Durasi = Tinggi
## 
##                 Media Sosial
## Kesehatan Mental Instagram TikTok
##            Baik         15     10
##            Buruk        85     90

# Fungsi Manual untuk Breslow-Day Test
breslow_day_test <- function(data) {
  K <- dim(data)[3]  # Jumlah strata (Durasi)
  OR_k <- numeric(K)  # Simpan OR per strata
  W_k <- numeric(K)   # Simpan weighting factor
  sum_WOR <- 0        # Numerator dari statistik BD
  sum_W <- 0          # Denominator dari statistik BD
  
  for (k in 1:K) {
    a <- data[1,1,k]
    b <- data[1,2,k]
    c <- data[2,1,k]
    d <- data[2,2,k]
    
    # Odds Ratio untuk tiap strata
    OR_k[k] <- (a * d) / (b * c)
    
    # Weighting Factor (Wald)
    W_k[k] <- (1/a + 1/b + 1/c + 1/d)^(-1)
    
    # Hitung numerator dan denominator
    sum_WOR <- sum_WOR + W_k[k] * log(OR_k[k])
    sum_W <- sum_W + W_k[k]
  }
  
  # Perhitungan Statistik Breslow-Day
  BD_stat <- sum((W_k * (log(OR_k) - sum_WOR/sum_W))^2) / sum_W
  
  # Approximate p-value (Chi-square dengan df = K - 1)
  p_value <- 1 - pchisq(BD_stat, df = K - 1)
  
  return(list("Breslow-Day Statistic" = BD_stat, "p-value" = p_value))
}

# Jalankan uji Breslow-Day
breslow_day_test(data_medsos)

## $`Breslow-Day Statistic`
## [1] 0.8541967
## 
## $`p-value`
## [1] 0.6523994

INTERPRETASI : p-value (0.6523994) > 0.05 → Odds Ratio dianggap homogen di seluruh durasi penggunaan media sosial.

8 Generalized Linear Model (GLM)

Model Linier Umum (GLM) merupakan perluasan dari regresi linier klasik. GLM memungkinkan kita memodelkan data ketika variabel respons tidak mengikuti distribusi normal dan/atau ketika hubungan antara variabel prediktor dengan ekspektasi respons bersifat non-linier.

GLM terdiri atas tiga komponen utama:

1. Distribusi dari keluarga eksponensial untuk variabel respons.
2. Fungsi link yang menghubungkan rata-rata respons dengan kombinasi linier dari prediktor.
3. Fungsi linear prediktor: \(\eta = X\beta\)

8.1 Keluarga Eksponensial

Distribusi termasuk dalam keluarga eksponensial jika dapat ditulis dalam bentuk:

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

Beberapa contoh distribusi yang termasuk dalam keluarga eksponensial antara lain:

• Normal
• Binomial
• Poisson
• Gamma

Contoh: Distribusi Binomial

Distribusi Binomial memiliki fungsi probabilitas:

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

Bentuk ini dapat ditulis ulang sebagai:

\[ 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 merupakan bagian dari keluarga eksponensial.

8.2 Regresi Logistik

Model regresi logistik mirip dengan regresi linier, tetapi output dibatasi pada nilai antara 0 dan 1 melalui fungsi aktivasi sigmoid. Model ini sering digunakan dalam klasifikasi biner.

Fungsi Sigmoid

Fungsi sigmoid dinyatakan sebagai:

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

Interpretasi hasil prediksi:

• Jika nilai fungsi sigmoid > 0.5, maka diklasifikasikan sebagai 1 (kelas positif).
• Jika < 0.5, maka diklasifikasikan sebagai 0 (kelas negatif).

Spesifikasi Model

Fungsi link (logit):

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

Estimasi parameter dilakukan dengan pendekatan Maximum Likelihood Estimation (MLE). Fungsi log-likelihood:

\[ 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^T \beta)}{1 + \exp(x_i^T \beta)} \]

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

8.3 Regresi Poisson

Model regresi Poisson digunakan ketika variabel respons adalah count data (bilangan bulat non-negatif). Model ini mengasumsikan bahwa variabel respons mengikuti distribusi Poisson.

Fungsi Probabilitas

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

Bentuk ini dapat ditulis sebagai keluarga eksponensial:

\[ 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!)\)

  Fungsi Link

Fungsi link kanonik:

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

Sehingga model menjadi:

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

Dan fungsi inverse:

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

Estimasi Parameter

Pendekatan MLE digunakan, dengan fungsi log-likelihood sebagai:

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

Parameter \(\beta\) diestimasi menggunakan iterasi numerik seperti Newton-Raphson.

Catatan Overdispersi

Jika nilai dispersi \(> 1\), kemungkinan terjadi overdispersi, sehingga model alternatif seperti Regresi Binomial Negatif dapat dipertimbangkan.

Model Poisson sangat bermanfaat untuk memodelkan data cacah, namun penting untuk mewaspadai potensi overdispersi dalam penerapannya.

9. Inferensi GLM

Generalized Linear Model (GLM) memerlukan pemahaman tentang ekspektasi dan varians dari estimator model. Hal ini menjadi dasar untuk melakukan inferensi statistik, seperti uji Wald, Likelihood Ratio Test, serta membentuk interval kepercayaan.

9.1 Ekspektasi dan Varians dalam GLM

9.1.1 Ekspektasi Estimator

Ekspektasi menunjukkan apakah suatu estimator bersifat tidak bias (unbiased), dengan definisi matematis:

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

Pada GLM, estimator \(\hat{\beta}\) yang diperoleh dari metode Maximum Likelihood Estimation (MLE) bersifat asymptotically unbiased.

9.1.2 Varians Estimator

Varians mengindikasikan presisi dari estimasi parameter. Dengan sampel besar, varians dari \(\hat{\beta}\) dapat didekati oleh:

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

Dengan \(\mathbf{W}\) merupakan matriks bobot yang tergantung pada distribusi dan fungsi link yang digunakan.

Estimator asimptotik mengikuti distribusi normal:

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

Distribusi ini digunakan dalam: - Uji Wald - Confidence Interval - Perhitungan P-value

9.1.3 Varians pada GLM Tidak Konstan

Berbeda dengan regresi linear biasa (OLS) yang mengasumsikan homoskedastisitas:

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

Dalam GLM, varian respon ditentukan oleh:

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

dengan \(\phi\) sebagai parameter dispersi dan \(V(\mu)\) sebagai fungsi varians.

---

9.1.4 Ekspektasi dan Varians Berdasarkan Fungsi Momen

Ekspektasi

Berdasarkan definisi fungsi momen:

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

Untuk distribusi dalam keluarga eksponensial:

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

Maka score function:

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

Dengan ekspektasi score:

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

Sehingga:

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

9.1.5 Varians

Turunan kedua log-likelihood:

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

Sehingga:

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

---

9.2 Metode Penaksiran Parameter

9.2.1 Maximum Likelihood Estimation (MLE)

MLE dilakukan dengan mencari nilai parameter yang memaksimalkan fungsi log-likelihood. Secara umum, langkah-langkahnya:

• Turunan pertama log-likelihood = 0
• Turunan kedua < 0 (maksimum)

Karena bentuk GLM tidak eksplisit, diperlukan metode numerik seperti:

9.2.2 Newton-Raphson

• Menggunakan vektor skor (gradien)
• Menggunakan matriks Hessian

Iterasi Newton-Raphson:

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

9.2.3 Fisher Scoring

• Modifikasi Newton-Raphson
• Mengganti Hessian dengan expected information matrix

9.2.4 IRLS (Iteratively Reweighted Least Square)

• Modifikasi Fisher Scoring
• Estimasi menyerupai metode Least Squares

---

9.3 Diagnostik Model GLM

9.3.1 Statistik Devians

Mengukur kecocokan model dibanding model saturated:

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

• Devians kecil → model cocok
• Devians besar → model tidak cocok

9.3.2 Statistik Chi-Kuadrat Pearson

Mengukur apakah model lebih baik dibandingkan tanpa model:

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

9.3.3 Analisis Residual

• Residual = selisih antara observasi dan prediksi
• Digunakan untuk mendeteksi penyimpangan sistematis
• Dapat digunakan untuk mengevaluasi asumsi model

9.4 Estimasi dan Inferensi Regresi Logistik

9.4.1 Fungsi Model Logistik

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

9.4.2 Log-likelihood:

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

9.4.3 Newton-Raphson:

1. Turunan Pertama (Score Function)

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

2. Turunan Kedua (Hessian Matrix)

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

3. Iterasi

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

---

9.4.4 Inferensi Parameter (Uji Wald)

Menggunakan uji Wald untuk menguji signifikansi parameter \(\beta_j\):

• \(H_0: \beta_j = 0\)
• \(H_1: \beta_j \neq 0\)

Jika \(\hat{\beta}_j\) mendekati distribusi normal:

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

Maka statistik Z:

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

Statistik Wald:

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

---

9.5 Estimasi dan Inferensi Regresi Poisson

9.5.1 Fungsi Distribusi Poisson

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

9.5.2 Model Regresi Poisson

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

9.5.3 Log-likelihood

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

9.5.4 Estimasi dengan IRLS

1. Definisi Model:

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

2. Log-likelihood: (sama seperti di atas)

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 - \mu}{\mu}\)
• \(\eta = \log(\lambda_i) = \mathbf{x}_i^\top \beta\)

10 KASUS DATA RILL

Data diambil dari Jurnal Dunia Kesmas tahun 2023

Data Survei Aspek Kehidupan Rumah Tangga Indonesia Tahun 2014-2015 (IFLS 5)

Ayah Perokok	Merokok (Ya)	Merokok (Tidak)	Total
Ya	54	42	96
Tidak	9	26	35

10.1 Tabel Kontingensi 2x2

tabel <- matrix(c(54, 42, 9, 26), nrow = 2, byrow = TRUE)
rownames(tabel) <- c("AyahPerokok_Ya", "AyahPerokok_Tidak")
colnames(tabel) <- c("Merokok_Ya", "Merokok_Tidak")
tabel

##                   Merokok_Ya Merokok_Tidak
## AyahPerokok_Ya            54            42
## AyahPerokok_Tidak          9            26

total <- sum(tabel)
total

## [1] 131

10.2 Peluang Bersama

peluang_bersama <- tabel / total
peluang_bersama

##                   Merokok_Ya Merokok_Tidak
## AyahPerokok_Ya    0.41221374     0.3206107
## AyahPerokok_Tidak 0.06870229     0.1984733

Peluang bahwa seseorang merokok dan memiliki ayah perokok adalah 0,4

Peluang bahwa seseorang tidak merokok dan memiliki ayah perokok adalah 0,06

Peluang bahwa seseorang merokok dan tidak memiliki ayah perokok adalah 0,32

Peluang bahwa seseorang tidak merokok dan tidak memiliki ayah perokok adalah 0,19

10.3 Peluang Marjinal

peluang_marjinal_baris <- rowSums(peluang_bersama)  # Ayah perokok
peluang_marjinal_kolom <- colSums(peluang_bersama)  # Kebiasaan merokok
peluang_marjinal_baris

##    AyahPerokok_Ya AyahPerokok_Tidak 
##         0.7328244         0.2671756

peluang_marjinal_kolom

##    Merokok_Ya Merokok_Tidak 
##      0.480916      0.519084

Peluang seseorang memiliki ayah perokok adalah 0,73

Peluang seseorang tidak memiliki ayah perokok adalah 0,26

Peluang seseorang merokok adalah 0,48

Peluang seseorang tidak merokok adalah 0,51

10.4 Peluang Bersyarat

a. Peluang Merokok berdasarkan Status Ayah

peluang_merokok_dengan_ayah <- tabel[ ,1] / rowSums(tabel)
names(peluang_merokok_dengan_ayah) <- c("AyahPerokok_Ya", "AyahPerokok_Tidak")
peluang_merokok_dengan_ayah

##    AyahPerokok_Ya AyahPerokok_Tidak 
##         0.5625000         0.2571429

Peluang seseorang merokok jika ayahnya perokok adalah 0,56

Peluang seseorang merokok jika ayahnya tidak merokok adalah 0,25

b. Peluang Ayah Perokok berdasarkan Status Merokok

peluang_ayah_dengan_merokok <- tabel[1, ] / colSums(tabel)
names(peluang_ayah_dengan_merokok) <- c("Merokok_Ya", "Merokok_Tidak")
peluang_ayah_dengan_merokok

##    Merokok_Ya Merokok_Tidak 
##     0.8571429     0.6176471

Peluang seseorang memiliki ayah perokok jika ia merokok adalah 0,85

Peluang seseorang memiliki ayah perokok jika ia tidak merokok adalah 0,61

10.5 Ukuran Asosiasi

# Data tabel kontingensi 2x2
tabel <- matrix(c(54, 42, 9, 26), nrow = 2, byrow = TRUE)
rownames(tabel) <- c("AyahPerokok_Ya", "AyahPerokok_Tidak")
colnames(tabel) <- c("Merokok_Ya", "Merokok_Tidak")
tabel

##                   Merokok_Ya Merokok_Tidak
## AyahPerokok_Ya            54            42
## AyahPerokok_Tidak          9            26

a <- tabel[1, 1]  # 54
b <- tabel[1, 2]  # 42
c <- tabel[2, 1]  # 9
d <- tabel[2, 2]  # 26

10.5.1 Risk Difference (RD)

risk_exposed <- a / (a + b)
risk_unexposed <- c / (c + d)
RD <- risk_exposed - risk_unexposed
RD

## [1] 0.3053571

kemungkinan seseorang menjadi perokok meningkat sebesar 40% jika ayahnya juga seorang perokok, dibandingkan dengan jika ayahnya tidak merokok.

10.5.2 Risk Ratio (RR)

RR <- risk_exposed / risk_unexposed
RR

## [1] 2.1875

risiko merokok bagi seseorang yang ayahnya perokok adalah tiga kali lebih besar dibandingkan dengan mereka yang ayahnya tidak merokok.

10.5.3 Odds Ratio (OR)

OR <- (a * d) / (b * c)
OR

## [1] 3.714286

Nilai Odds Ratio sebesar 6 menunjukkan bahwa peluang (odds) seseorang menjadi perokok adalah enam kali lebih tinggi jika ayahnya adalah perokok dibandingkan dengan jika ayahnya bukan perokok.

10.6 Estimasi

10.6.1 Estimasi Titik

# Data
x <- 63        # jumlah individu yang merokok
n <- 131       # total sampel
p_hat <- x / n # estimasi titik proporsi

# Estimasi titik
p_hat

## [1] 0.480916

Artinya, dari total 131 responden, sekitar 48,1% diketahui memiliki kebiasaan merokok.

10.6.2 Estimasi Interval

# Estimasi interval (confidence interval 95%)
z <- 1.96
se <- sqrt(p_hat * (1 - p_hat) / n)
lower <- p_hat - z * se
upper <- p_hat + z * se

# Hasil interval
c(lower, upper)

## [1] 0.3953554 0.5664766

Dengan tingkat kepercayaan 95%, kita dapat menyatakan bahwa proporsi sebenarnya dari populasi yang merokok diperkirakan berada antara 39,5% hingga 56,6%.

Artinya, jika kita mengambil banyak sampel acak yang berukuran sama, sekitar 95% dari interval yang dihitung seperti ini akan mengandung proporsi sebenarnya dari populasi.

10.7 Uji Hipotesis

# Data dari tabel kontingensi
n11 <- 54   # Merokok, Ayah perokok
n12 <- 42   # Tidak merokok, Ayah perokok
n21 <- 9    # Merokok, Ayah tidak perokok
n22 <- 26   # Tidak merokok, Ayah tidak perokok

# Ukuran sampel tiap kelompok
n1_ <- n11 + n12  # Total dengan ayah perokok
n2_ <- n21 + n22  # Total dengan ayah tidak perokok

10.7.1 Uji Proporsi masing-masing kelompok

p1_hat <- n11 / n1_  # proporsi merokok dengan ayah perokok
p2_hat <- n21 / n2_  # proporsi merokok dengan ayah tidak perokok

10.7.2 Proporsi gabungan

p_hat <- (n11 + n21) / (n1_ + n2_)
p_hat

## [1] 0.480916

10.7.3 Statistik uji Z

z <- (p1_hat - p2_hat) / sqrt(p_hat * (1 - p_hat) * (1/n1_ + 1/n2_))
z

## [1] 3.095199

Nilai statistik uji diperoleh sebesar

Z=3,10. Dengan menggunakan tingkat signifikansi α=0,05, nilai kritis dari distribusi normal standar adalah ±1,96. maka hipotesis nol yang menyatakan bahwa tidak ada perbedaan proporsi antara dua kelompok ditolak. Kesimpulan: Terdapat bukti statistik yang signifikan untuk menyatakan bahwa proporsi remaja yang merokok berbeda secara signifikan antara kelompok yang memiliki ayah perokok dan yang tidak.

10.8 Hipotesis Uji dalam Tabel Kontingensi 2x2

# Data kontingensi
a <- 54  # Ayah perokok, anak merokok
b <- 42  # Ayah perokok, anak tidak merokok
c <- 9   # Ayah tidak merokok, anak merokok
d <- 26  # Ayah tidak merokok, anak tidak merokok

# Ukuran sampel tiap grup
n1 <- a + b  # Total ayah perokok
n2 <- c + d  # Total ayah tidak perokok

# Proporsi masing-masing
p1 <- a / n1
p2 <- c / n2

10.8.1 Risk Difference (RD)

RD <- p1 - p2
SE_RD <- sqrt((p1 * (1 - p1)) / n1 + (p2 * (1 - p2)) / n2)
Z_RD <- RD / SE_RD

10.8.2 Relative Risk (RR)

RR <- p1 / p2
SE_lnRR <- sqrt((1/a) - (1/n1) + (1/c) - (1/n2))
Z_RR <- log(RR) / SE_lnRR

10.8.3 Odds Ratio (OR)

OR <- (a * d) / (b * c)
SE_lnOR <- sqrt(1/a + 1/b + 1/c + 1/d)
Z_OR <- log(OR) / SE_lnOR

Tampilkan hasil

list(
  Risk_Difference = RD,
  SE_RD = SE_RD,
  Z_RD = Z_RD,
  
  Relative_Risk = RR,
  SE_lnRR = SE_lnRR,
  Z_RR = Z_RR,
  
  Odds_Ratio = OR,
  SE_lnOR = SE_lnOR,
  Z_OR = Z_OR
)

## $Risk_Difference
## [1] 0.3053571
## 
## $SE_RD
## [1] 0.08956117
## 
## $Z_RD
## [1] 3.409482
## 
## $Relative_Risk
## [1] 2.1875
## 
## $SE_lnRR
## [1] 0.3010673
## 
## $Z_RR
## [1] 2.599948
## 
## $Odds_Ratio
## [1] 3.714286
## 
## $SE_lnOR
## [1] 0.4380647
## 
## $Z_OR
## [1] 2.995417

1. Risk Difference (RD) Nilai RD: 0.305 Z-statistik RD: 3.41

Interpretasi: Terdapat perbedaan absolut sebesar 30.5% dalam kemungkinan merokok antara anak yang memiliki ayah perokok dibandingkan yang tidak. Nilai Z_RD = 3.41 menunjukkan bahwa perbedaan ini signifikan secara statistik pada tingkat kepercayaan 95% (karena |Z| > 1.96), sehingga kita menolak hipotesis nol bahwa tidak ada perbedaan.

2. Relative Risk (RR) Nilai RR: 2.19 Z-statistik RR: 2.60

Interpretasi: Anak dari ayah perokok memiliki risiko merokok yang 2.19 kali lebih besar dibandingkan anak dari ayah non-perokok. Dengan nilai Z sebesar 2.60, hasil ini juga signifikan secara statistik pada taraf signifikansi 5%, mendukung adanya asosiasi kuat antara kebiasaan merokok ayah dan anak.

3. Odds Ratio (OR) Nilai OR: 3.71 Z-statistik OR: 3.00

Interpretasi: Peluang anak merokok jika ayahnya perokok adalah 3.71 kali lebih besar dibandingkan jika ayahnya tidak merokok. Z_OR = 3.00 > 1.96 menunjukkan bahwa OR juga signifikan secara statistik, yang memperkuat kesimpulan adanya hubungan yang kuat dan signifikan antara variabel.

10.9 Uji Independensi

10.9.1 Uji chi-square-test

# Data kontingensi
kontingensi <- matrix(c(54, 42, 9, 26), 
                      nrow = 2, 
                      byrow = TRUE,
                      dimnames = list("Ayah" = c("Perokok", "Tidak Perokok"),
                                      "Anak" = c("Merokok", "Tidak Merokok")))

# Menampilkan tabel
kontingensi

##                Anak
## Ayah            Merokok Tidak Merokok
##   Perokok            54            42
##   Tidak Perokok       9            26

# Uji Chi-Square
uji_chi <- chisq.test(kontingensi, correct = FALSE)
uji_chi

## 
##  Pearson's Chi-squared test
## 
## data:  kontingensi
## X-squared = 9.5803, df = 1, p-value = 0.001967

(p-value) = 0.001967 < 0,05 TOLAK H0 Terdapat hubungan yang signifikan secara statistik antara status merokok ayah dan kebiasaan merokok anak. Artinya, kebiasaan merokok seorang ayah mempengaruhi kemungkinan anaknya untuk merokok.

10.10 Partisi Chi-Square

Data Bersumber dari buku “An Introduction To Statistical Modelling”

Data Uji Coba pemberian obat dan plasebo

Treatment	Improved	No difference	Worse	Total
Placebo	18	17	15	50
Half dose	20	10	20	50
Full dose	25	13	12	50

# Tabel 1: Randomized controlled trial
trial_data <- matrix(c(18, 17, 15,
                       20, 10, 20,
                       25, 13, 12),
                     nrow = 3,
                     byrow = TRUE,
                     dimnames = list(
                       Treatment = c("Placebo", "Half dose", "Full dose"),
                       Response = c("Improved", "No difference", "Worse")))

trial_data

##            Response
## Treatment   Improved No difference Worse
##   Placebo         18            17    15
##   Half dose       20            10    20
##   Full dose       25            13    12

# Chi-Square
chisq.test(trial_data)

## 
##  Pearson's Chi-squared test
## 
## data:  trial_data
## X-squared = 5.1732, df = 4, p-value = 0.27

# G-squared
if (!require(DescTools)) install.packages("DescTools")

## Loading required package: DescTools

## Warning: package 'DescTools' was built under R version 4.3.3

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

library(DescTools)
GTest(trial_data)

## 
##  Log likelihood ratio (G-test) test of independence without correction
## 
## data:  trial_data
## G = 5.1291, X-squared df = 4, p-value = 0.2743

# Fisher test (tidak valid karena bukan tabel 2x2, tapi kita coba untuk lihat errornya)
tryCatch({
  fisher.test(trial_data)
}, error = function(e) e)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  trial_data
## p-value = 0.2885
## alternative hypothesis: two.sided

1. Pearson Chi-Square p-value = 0.27 > 0.05

Artinya: tidak ada cukup bukti untuk menolak hipotesis nol.

Kesimpulan: Tidak ada hubungan yang signifikan antara variabel-variabel dalam tabel

2. G-test p-value = 0.2743 > 0.05

Artinya: hasilnya sama dengan chi-squared test.

Kesimpulan: Tidak ada asosiasi yang signifikan antara variabel-variabel dalam tabel.

3. Fisher’s Exact Test p-value = 0.2885 > 0.05

Kesimpulan: Fisher’s Exact Test juga menunjukkan tidak ada hubungan signifikan antara variabel.

10.11 Uji Likelihood Ratio

# Data kontingensi
kontingensi <- matrix(c(54, 42, 9, 26),
                      nrow = 2,
                      byrow = TRUE,
                      dimnames = list(
                        "Ayah" = c("Perokok", "Tidak Perokok"),
                        "Anak" = c("Merokok", "Tidak Merokok")
                      ))

# Uji Likelihood Ratio (G-Test) menggunakan package 'DescTools'
# Install jika belum tersedia
# install.packages("DescTools")
library(DescTools)

# Lakukan uji Likelihood Ratio (G-Test)
uji_likelihood <- GTest(kontingensi)

# Tampilkan hasil
uji_likelihood

## 
##  Log likelihood ratio (G-test) test of independence without correction
## 
## data:  kontingensi
## G = 9.93, X-squared df = 1, p-value = 0.001626

(p-value) = 0.001626 < 0,05 TOLAK H0 Terdapat hubungan yang signifikan secara statistik antara status merokok ayah dan kebiasaan merokok anak. Artinya, anak dari ayah yang merokok memiliki kecenderungan lebih tinggi untuk merokok, dibandingkan dengan anak dari ayah yang tidak merokok.

10.12 Uji Exact Fisher

uji_fisher <- fisher.test(kontingensi)

# Tampilkan hasil
uji_fisher

## 
##  Fisher's Exact Test for Count Data
## 
## data:  kontingensi
## p-value = 0.002776
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.477752 9.917529
## sample estimates:
## odds ratio 
##   3.676977

(p-value) = 0.002776 < 0,05 TOLAK H0 maka H₀ ditolak pada tingkat signifikansi 5%. Artinya, terdapat hubungan yang signifikan secara statistik antara status merokok ayah dan kebiasaan merokok anak.

Anak dari ayah yang merokok memiliki peluang merokok sekitar 3,68 kali lebih tinggi dibandingkan dengan anak dari ayah yang tidak merokok.

Interval kepercayaan 95% untuk odds ratio (1.478 hingga 9.918) tidak mencakup angka 1, yang memperkuat bukti adanya asosiasi yang bermakna antara dua variabel tersebut.

10.13 Pearson Residual

# Total
N <- sum(kontingensi)

# Nilai harapan (expected values)
expected <- chisq.test(kontingensi)$expected

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

##                Anak
## Ayah              Merokok Tidak Merokok
##   Perokok        1.152672     -1.109486
##   Tidak Perokok -1.909007      1.837483

# Hitung probabilitas marginal
row_prop <- rowSums(kontingensi) / N
col_prop <- colSums(kontingensi) / N

# Buat matriks kosong untuk standardized residual
standardized_residual <- matrix(0, nrow = 2, ncol = 2)
rownames(standardized_residual) <- rownames(kontingensi)
colnames(standardized_residual) <- colnames(kontingensi)

# Hitung Standardized Residual (Adjusted Residual)
for (i in 1:2) {
  for (j in 1:2) {
    pij_row <- row_prop[i]
    pij_col <- col_prop[j]
    E_ij <- expected[i, j]
    O_ij <- kontingensi[i, j]
    standardized_residual[i, j] <- (O_ij - E_ij) / sqrt(E_ij * (1 - pij_row) * (1 - pij_col))
  }
}
standardized_residual

##                 Merokok Tidak Merokok
## Perokok        3.095199     -3.095199
## Tidak Perokok -3.095199      3.095199

Perokok - Merokok (3.095199): Residual positif besar (> 2) → jumlah perokok yang merokok jauh lebih banyak daripada yang diharapkan. Ini signifikan.

Perokok - Tidak Merokok (-3.095199): Residual negatif besar (< -2) → jumlah perokok yang tidak merokok jauh lebih sedikit dari yang diharapkan. Ini signifikan.

Tidak Perokok - Merokok (-3.095199): Residual negatif besar (< -2) → jumlah tidak perokok yang merokok jauh lebih sedikit dari yang diharapkan. Ini signifikan.

Tidak Perokok - Tidak Merokok (3.095199): Residual positif besar (> 2) → jumlah tidak perokok yang tidak merokok jauh lebih banyak dari yang diharapkan. Ini signifikan.

10.14 Tabel Kontingensi 3 Arah

Data diambil dari Jurnal Dunia Kesmas tahun 2023

Data Survei Aspek Kehidupan Rumah Tangga Indonesia Tahun 2014-2015 (IFLS 5)

\[ \begin{array}{|c|c|c|c|} \hline \textbf{Pengetahuan} & \textbf{Uang Saku} & \textbf{Merokok (Ya)} & \textbf{Merokok (Tidak)} \\ \hline \text{Rendah} & \text{Ya} & 15 & 10 \\ \text{Rendah} & \text{Tidak} & 10 & 10 \\ \text{Tinggi} & \text{Ya} & 21 & 15 \\ \text{Tinggi} & \text{Tidak} & 17 & 33 \\ \hline \end{array} \]

10.14.1 Tabel Kontingensi

# Membuat tabel kontingensi
library(knitr)

## Warning: package 'knitr' was built under R version 4.3.3

data <- matrix(c(15, 10, 
                 10, 10, 
                 21, 15, 
                 17, 33), 
               nrow = 4, byrow = TRUE)

rownames(data) <- c("Rendah - Ya", "Rendah - Tidak", "Tinggi - Ya", "Tinggi - Tidak")
colnames(data) <- c("Kebiasaan Merokok - Ya", "Kebiasaan Merokok - Tidak")

kable(data, caption = "Tabel Kontingensi Pengetahuan, Uang Saku, dan Kebiasaan Merokok")

Tabel Kontingensi Pengetahuan, Uang Saku, dan Kebiasaan Merokok

	Kebiasaan Merokok - Ya	Kebiasaan Merokok - Tidak
Rendah - Ya	15	10
Rendah - Tidak	10	10
Tinggi - Ya	21	15
Tinggi - Tidak	17	33

10.14.2 Tabel Frekuensi Parsial

# Menghitung frekuensi parsial (persentase per baris)
parsial <- prop.table(data, margin = 1)

kable(round(parsial, 3), caption = "Tabel Frekuensi Parsial (Per Baris)")

Tabel Frekuensi Parsial (Per Baris)

	Kebiasaan Merokok - Ya	Kebiasaan Merokok - Tidak
Rendah - Ya	0.600	0.400
Rendah - Tidak	0.500	0.500
Tinggi - Ya	0.583	0.417
Tinggi - Tidak	0.340	0.660

10.14.3 Tabel Frekuensi Marginal

# Menghitung frekuensi marginal
marginal_rows <- margin.table(data, margin = 1)
marginal_cols <- margin.table(data, margin = 2)

# Membuat tabel marginal
marginal_table <- addmargins(data)

kable(marginal_table, caption = "Tabel Frekuensi Marginal (Penjumlahan Baris dan Kolom)")

Tabel Frekuensi Marginal (Penjumlahan Baris dan Kolom)

	Kebiasaan Merokok - Ya	Kebiasaan Merokok - Tidak	Sum
Rendah - Ya	15	10	25
Rendah - Tidak	10	10	20
Tinggi - Ya	21	15	36
Tinggi - Tidak	17	33	50
Sum	63	68	131

# Peluang bersama
joint_prob <- prop.table(data)

kable(round(joint_prob, 3), caption = "Tabel Peluang Bersama")

Tabel Peluang Bersama

	Kebiasaan Merokok - Ya	Kebiasaan Merokok - Tidak
Rendah - Ya	0.115	0.076
Rendah - Tidak	0.076	0.076
Tinggi - Ya	0.160	0.115
Tinggi - Tidak	0.130	0.252

10.15 Peluang Marjinal

# Peluang marjinal baris
marginal_row <- prop.table(margin.table(data, margin = 1))

# Peluang marjinal kolom
marginal_col <- prop.table(margin.table(data, margin = 2))

# Menampilkan
kable(round(marginal_row, 3), caption = "Peluang Marjinal Baris")

Peluang Marjinal Baris

	x
Rendah - Ya	0.191
Rendah - Tidak	0.153
Tinggi - Ya	0.275
Tinggi - Tidak	0.382

kable(round(marginal_col, 3), caption = "Peluang Marjinal Kolom")

Peluang Marjinal Kolom

	x
Kebiasaan Merokok - Ya	0.481
Kebiasaan Merokok - Tidak	0.519

10.16 Peluang Bersyarat

# Peluang bersyarat: P(Kebiasaan Merokok | Pengetahuan dan Uang Saku)
conditional_prob_row <- prop.table(data, margin = 1)

# Peluang bersyarat: P(Pengetahuan dan Uang Saku | Kebiasaan Merokok)
conditional_prob_col <- prop.table(data, margin = 2)

kable(round(conditional_prob_row, 3), caption = "Peluang Bersyarat: P(Kebiasaan Merokok | Pengetahuan dan Uang Saku)")

Peluang Bersyarat: P(Kebiasaan Merokok | Pengetahuan dan Uang Saku)

	Kebiasaan Merokok - Ya	Kebiasaan Merokok - Tidak
Rendah - Ya	0.600	0.400
Rendah - Tidak	0.500	0.500
Tinggi - Ya	0.583	0.417
Tinggi - Tidak	0.340	0.660

kable(round(conditional_prob_col, 3), caption = "Peluang Bersyarat: P(Pengetahuan dan Uang Saku | Kebiasaan Merokok)")

Peluang Bersyarat: P(Pengetahuan dan Uang Saku | Kebiasaan Merokok)

	Kebiasaan Merokok - Ya	Kebiasaan Merokok - Tidak
Rendah - Ya	0.238	0.147
Rendah - Tidak	0.159	0.147
Tinggi - Ya	0.333	0.221
Tinggi - Tidak	0.270	0.485

10.17 Ukuran Asosiasi

# Membuat tabel kontingensi 2x2
library(knitr)

# Misal kita ambil: Tinggi - Tidak vs Rendah - Ya
#                  Merokok | Tidak Merokok
data2x2 <- matrix(c(17, 33, 
                    15, 10), 
                  nrow = 2, byrow = TRUE)

rownames(data2x2) <- c("Tinggi - Tidak", "Rendah - Ya")
colnames(data2x2) <- c("Kebiasaan Merokok - Ya", "Kebiasaan Merokok - Tidak")

kable(data2x2, caption = "Tabel 2x2 untuk Perhitungan Ukuran Asosiasi")

Tabel 2x2 untuk Perhitungan Ukuran Asosiasi

	Kebiasaan Merokok - Ya	Kebiasaan Merokok - Tidak
Tinggi - Tidak	17	33
Rendah - Ya	15	10

Perhitungan RD, RR, OR

# Ambil data
a <- data2x2[1,1] # Eksposur & Outcome (+/+)
b <- data2x2[1,2] # Eksposur (+), Outcome (-)
c <- data2x2[2,1] # Non-Eksposur (-), Outcome (+)
d <- data2x2[2,2] # Non-Eksposur (-), Outcome (-)

# Hitung proporsi
risk_exposed <- a / (a + b)
risk_unexposed <- c / (c + d)

# Risk Difference (RD)
RD <- risk_exposed - risk_unexposed

# Relative Risk (RR)
RR <- risk_exposed / risk_unexposed

# Odds Ratio (OR)
OR <- (a*d) / (b*c)

# Hasil
RD

## [1] -0.26

RR

## [1] 0.5666667

OR

## [1] 0.3434343

Interpretasi

Risk Difference (RD) = -0,26

Artinya, risiko merokok pada kelompok ‘Tinggi - Tidak’ lebih rendah 26% dibandingkan kelompok ‘Rendah - Ya’.

Karena RD bernilai negatif, kelompok ‘Tinggi - Tidak’ lebih terlindungi terhadap kebiasaan merokok dibandingkan kelompok ‘Rendah - Ya’.

Relative Risk (RR) = 0,57

Risiko merokok pada kelompok ‘Tinggi - Tidak’ adalah sekitar 56,7% dari risiko merokok pada kelompok ‘Rendah - Ya’.

Karena RR < 1, ini menunjukkan bahwa memiliki pengetahuan tinggi dan tidak memiliki uang saku berkaitan dengan penurunan risiko merokok.

Odds Ratio (OR) = 0,34

Peluang (odds) merokok pada kelompok ‘Tinggi - Tidak’ sekitar 34% dibandingkan dengan peluang merokok pada kelompok ‘Rendah - Ya’.

Karena OR < 1, ini juga menunjukkan bahwa kelompok ‘Tinggi - Tidak’ cenderung lebih kecil kemungkinannya untuk merokok dibandingkan ‘Rendah - Ya’.

10.18 Pengujian Conditional Independence

# Buat data
data <- data.frame(
  Pengetahuan = rep(c("Rendah", "Tinggi"), each = 4),
  Uang_Saku = rep(c("Ya", "Tidak"), each = 2, times = 2),
  Kebiasaan_Merokok = rep(c("Ya", "Tidak"), times = 4),
  Frekuensi = c(15, 10, 10, 10, 21, 15, 17, 33)
)

# Lihat data
print(data)

##   Pengetahuan Uang_Saku Kebiasaan_Merokok Frekuensi
## 1      Rendah        Ya                Ya        15
## 2      Rendah        Ya             Tidak        10
## 3      Rendah     Tidak                Ya        10
## 4      Rendah     Tidak             Tidak        10
## 5      Tinggi        Ya                Ya        21
## 6      Tinggi        Ya             Tidak        15
## 7      Tinggi     Tidak                Ya        17
## 8      Tinggi     Tidak             Tidak        33

# Bikin tabel kontingensi untuk setiap strata Uang Saku
library(dplyr)

## Warning: package 'dplyr' was built under R version 4.3.3

## 
## Attaching package: 'dplyr'

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

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

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

# Pisahkan berdasarkan Uang Saku
data_ya <- data %>% filter(Uang_Saku == "Ya")
data_tidak <- data %>% filter(Uang_Saku == "Tidak")

# Bentuk tabel kontingensi untuk 'Ya'
table_ya <- xtabs(Frekuensi ~ Pengetahuan + Kebiasaan_Merokok, data = data_ya)
print(table_ya)

##            Kebiasaan_Merokok
## Pengetahuan Tidak Ya
##      Rendah    10 15
##      Tinggi    15 21

# Bentuk tabel kontingensi untuk 'Tidak'
table_tidak <- xtabs(Frekuensi ~ Pengetahuan + Kebiasaan_Merokok, data = data_tidak)
print(table_tidak)

##            Kebiasaan_Merokok
## Pengetahuan Tidak Ya
##      Rendah    10 10
##      Tinggi    33 17

# Uji chi-square pada masing-masing strata
chi_ya <- chisq.test(table_ya, correct = FALSE)
chi_tidak <- chisq.test(table_tidak, correct = FALSE)

# Tampilkan hasil
chi_ya

## 
##  Pearson's Chi-squared test
## 
## data:  table_ya
## X-squared = 0.016944, df = 1, p-value = 0.8964

chi_tidak

## 
##  Pearson's Chi-squared test
## 
## data:  table_tidak
## X-squared = 1.5435, df = 1, p-value = 0.2141

Strata Uang Saku = “Ya”:

Nilai p-value = 0,8964 (>> 0,05).

Artinya, Pengetahuan dan Kebiasaan Merokok adalah independen secara kondisional terhadap Uang Saku = Ya.

Kesimpulan: Tidak ada hubungan yang signifikan antara tingkat pengetahuan dan kebiasaan merokok pada siswa yang memiliki uang saku.

Strata Uang Saku = “Tidak”:

Nilai p-value = 0,2141 (> 0,05).

Artinya, Pengetahuan dan Kebiasaan Merokok juga independen secara kondisional terhadap Uang Saku = Tidak.

Kesimpulan: Tidak ada hubungan yang signifikan antara tingkat pengetahuan dan kebiasaan merokok pada siswa yang tidak memiliki uang saku.

10.20 Uji Cochran-Mantel-Haenszel

# Data awal (berdasarkan tabel sebelumnya)
data <- array(
  c(15, 10,   # Pengetahuan Rendah, Uang Saku Ya
    10, 10,   # Pengetahuan Rendah, Uang Saku Tidak
    21, 15,   # Pengetahuan Tinggi, Uang Saku Ya
    17, 33),  # Pengetahuan Tinggi, Uang Saku Tidak
  dim = c(2, 2, 2),
  dimnames = list(
    "Pengetahuan" = c("Rendah", "Tinggi"),
    "Kebiasaan Merokok" = c("Ya", "Tidak"),
    "Uang Saku" = c("Ya", "Tidak")
  )
)

# Tampilkan data
data

## , , Uang Saku = Ya
## 
##            Kebiasaan Merokok
## Pengetahuan Ya Tidak
##      Rendah 15    10
##      Tinggi 10    10
## 
## , , Uang Saku = Tidak
## 
##            Kebiasaan Merokok
## Pengetahuan Ya Tidak
##      Rendah 21    17
##      Tinggi 15    33

# Uji Cochran-Mantel-Haenszel
uji_cmh <- mantelhaen.test(data)

# Tampilkan hasil
uji_cmh

## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  data
## Mantel-Haenszel X-squared = 4.0528, df = 1, p-value = 0.0441
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.084462 4.446890
## sample estimates:
## common odds ratio 
##          2.196015

Karena p-value = 0.0441 < 0.05, maka kita menolak hipotesis nol.

Ini berarti ada asosiasi yang signifikan antara Pengetahuan dan Kebiasaan Merokok, setelah mengontrol variabel Uang Saku.

Common Odds Ratio sebesar 2.196 menunjukkan bahwa individu dengan pengetahuan rendah memiliki kemungkinan sekitar 2,2 kali lebih besar untuk memiliki kebiasaan merokok dibandingkan individu dengan pengetahuan tinggi, setelah mengontrol pengaruh uang saku.

Karena interval kepercayaan (1.084, 4.447) tidak melewati 1, ini memperkuat bukti bahwa asosiasinya signifikan.

10.21 Uji Homogenitas

# Pastikan package DescTools terpasang
if (!require(DescTools)) install.packages("DescTools")
library(DescTools)

# Buat 2 strata berdasarkan Pengetahuan (Rendah dan Tinggi)

# Strata 1: Pengetahuan Rendah
# Baris: Uang Saku = Ya, Tidak
# Kolom: Merokok = Ya, Tidak
stratum1 <- matrix(c(15, 10, 10, 10), nrow = 2, byrow = TRUE)

# Strata 2: Pengetahuan Tinggi
stratum2 <- matrix(c(21, 15, 17, 33), nrow = 2, byrow = TRUE)

# Gabungkan jadi array 3D: [baris, kolom, strata]
data_array <- array(c(stratum1, stratum2),
  dim = c(2, 2, 2),
  dimnames = list(
    UangSaku = c("Ya", "Tidak"),
    Merokok = c("Ya", "Tidak"),
    Pengetahuan = c("Rendah", "Tinggi")
  )
)

# Lihat datanya
data_array

## , , Pengetahuan = Rendah
## 
##         Merokok
## UangSaku Ya Tidak
##    Ya    15    10
##    Tidak 10    10
## 
## , , Pengetahuan = Tinggi
## 
##         Merokok
## UangSaku Ya Tidak
##    Ya    21    15
##    Tidak 17    33

# Jalankan Uji Breslow-Day
breslow_result <- BreslowDayTest(data_array)

# Tampilkan hasil
print(breslow_result)

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

Nilai p-value = 0.4307 > 0.05

Artinya, tidak cukup bukti untuk menolak H₀

Kesimpulan: Odds ratio antara kelompok Pengetahuan Rendah dan Tinggi dapat dianggap homogen. Artinya, efek uang saku terhadap kebiasaan merokok konsisten di kedua strata pengetahuan.

10.22 Estimasi Parameter pada GLM

# Buat data frame berdasarkan tabel kontingensi
data <- data.frame(
  Merokok = c(rep(1,15), rep(0,10),  # Rendah, Uang Saku Ya
              rep(1,10), rep(0,10),  # Rendah, Uang Saku Tidak
              rep(1,21), rep(0,15),  # Tinggi, Uang Saku Ya
              rep(1,17), rep(0,33)), # Tinggi, Uang Saku Tidak
  UangSaku = c(rep(1,25), rep(0,20), rep(1,36), rep(0,50)),
  Pengetahuan = c(rep(0,45), rep(1,86))  # 0 = Rendah, 1 = Tinggi
)
data

# Fit regresi logistik
model <- glm(Merokok ~ UangSaku + Pengetahuan, data = data, family = binomial(link = "logit"))

# Lihat hasil estimasi beta
summary(model)

## 
## Call:
## glm(formula = Merokok ~ UangSaku + Pengetahuan, family = binomial(link = "logit"), 
##     data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.2086     0.3621  -0.576   0.5646  
## UangSaku      0.7907     0.3609   2.191   0.0284 *
## Pengetahuan  -0.3634     0.3792  -0.958   0.3379  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 181.41  on 130  degrees of freedom
## Residual deviance: 175.00  on 128  degrees of freedom
## AIC: 181
## 
## Number of Fisher Scoring iterations: 4

1. Uang Saku (p = 0.0284 < 0.05) Signifikan secara statistik.

Koefisien logit: 0.7907 → odds ratio = exp(0.7907) ≈ 2.20

Interpretasi: Siswa yang memiliki uang saku memiliki peluang 2,2 kali lebih besar untuk merokok dibanding yang tidak punya uang saku, dengan asumsi tingkat pengetahuan yang sama.

2. Pengetahuan (p = 0.3379 > 0.05) Tidak signifikan.

Koefisien negatif → indikasi bahwa siswa dengan pengetahuan tinggi cenderung lebih kecil peluangnya merokok, tapi secara statistik belum cukup kuat (p > 0.05).Odds ratio = exp(-0.3634) ≈ 0.70

Plot Probabilitas

new_data <- data.frame(
  UangSaku = seq(min(data$UangSaku), max(data$UangSaku), length.out = 100),
  Pengetahuan = mean(data$Pengetahuan, na.rm = TRUE)
)

# Prediksi probabilitas
new_data$prob <- predict(model, newdata = new_data, type = "response")

# Plot
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.3

ggplot(new_data, aes(x = UangSaku, y = prob)) +
  geom_line(color = "blue", size = 1.2) +
  labs(
    title = "Probabilitas Merokok berdasarkan Uang Saku",
    x = "Uang Saku",
    y = "Probabilitas Merokok"
  ) +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

10.23 Estimasi Regresi Poisson

# Model regresi Poisson
model_poisson <- glm(Merokok ~ UangSaku + Pengetahuan, 
                     family = poisson(link = "log"), data = data)

summary(model_poisson)

## 
## Call:
## glm(formula = Merokok ~ UangSaku + Pengetahuan, family = poisson(link = "log"), 
##     data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -0.8314     0.2608  -3.188  0.00143 **
## UangSaku      0.4031     0.2568   1.570  0.11646   
## Pengetahuan  -0.1742     0.2598  -0.671  0.50247   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 92.240  on 130  degrees of freedom
## Residual deviance: 88.969  on 128  degrees of freedom
## AIC: 220.97
## 
## Number of Fisher Scoring iterations: 5

μ=exp(−0.8314)≈0.435 Artinya, rata-rata jumlah yang merokok diperkirakan sekitar 0.435 orang jika uang saku dan pengetahuan = 0.

Signifikansi: Sangat signifikan (p = 0.00143, ditandai dengan **), artinya intercept ini secara statistik bermakna.

Plot Prediksi

# Buat data baru untuk prediksi
new_data <- data.frame(
  UangSaku = seq(min(data$UangSaku), max(data$UangSaku), length.out = 100),
  Pengetahuan = mean(data$Pengetahuan, na.rm = TRUE)
)

# Prediksi dari model Poisson
new_data$pred <- predict(model_poisson, newdata = new_data, type = "response")

# Visualisasi
library(ggplot2)
ggplot(new_data, aes(x = UangSaku, y = pred)) +
  geom_line(color = "darkorange", size = 1.2) +
  labs(
    title = "Prediksi Jumlah Merokok berdasarkan Uang Saku",
    x = "Uang Saku",
    y = "Prediksi Jumlah Merokok (μ̂)"
  ) +
  theme_minimal()

10.24 Regresi Poisson GLM

model_poisson <- glm(Merokok ~ UangSaku + Pengetahuan, 
                     family = poisson(link = "log"), data = data)

# Hasil ringkasan model
summary(model_poisson)

## 
## Call:
## glm(formula = Merokok ~ UangSaku + Pengetahuan, family = poisson(link = "log"), 
##     data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -0.8314     0.2608  -3.188  0.00143 **
## UangSaku      0.4031     0.2568   1.570  0.11646   
## Pengetahuan  -0.1742     0.2598  -0.671  0.50247   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 92.240  on 130  degrees of freedom
## Residual deviance: 88.969  on 128  degrees of freedom
## AIC: 220.97
## 
## Number of Fisher Scoring iterations: 5

1. UangSaku Estimate: 0.4031

Interpretasi: Setiap kenaikan 1 unit UangSaku → log(μ) naik 0.4031, atau:

μ=exp(−0.4031)≈1.497

Artinya, jumlah merokok diprediksi meningkat sekitar 49.7% untuk setiap tambahan 1 unit UangSaku, dengan asumsi Pengetahuan tetap.

Signifikansi: p-value = 0.11646 → tidak signifikan pada tingkat 5%. Jadi, meskipun ada kecenderungan arah hubungan positif, tidak cukup bukti statistik bahwa UangSaku berpengaruh secara nyata terhadap jumlah merokok.

2. Pengetahuan Estimate: -0.1742

Interpretasi: Setiap kenaikan 1 unit Pengetahuan → log(μ) turun 0.1742, atau:

μ=exp(0.17421)≈0.840 Artinya, peningkatan Pengetahuan diprediksi menurunkan jumlah merokok sekitar 16%, dengan asumsi variabel lain konstan.

Signifikansi: p-value = 0.50247 → sangat tidak signifikan, menunjukkan bahwa Pengetahuan tidak memiliki pengaruh nyata dalam model ini.

10.25 Ekspektasi dan Varians GLM

# Ekspektasi (mean prediksi μ̂)
mu_hat <- predict(model_poisson, type = "response")

# Varians (untuk Poisson, varians = mean)
var_hat <- mu_hat

# Tambahkan ke data
data$mu_hat <- mu_hat
data$var_hat <- var_hat

# Lihat beberapa hasil
data[, c("UangSaku", "Pengetahuan", "Merokok", "mu_hat", "var_hat")]

10.26 Metode Penaksiran Parameter

# X: Desain matriks (dengan intercept)
X <- as.matrix(cbind(1, data$UangSaku, data$Pengetahuan))

# y: respon
y <- data$Merokok

# Inisialisasi beta (semua nol)
beta <- matrix(0, nrow = ncol(X), ncol = 1)

# Fungsi iterasi Newton-Raphson
newton_raphson_poisson <- function(X, y, beta_init, tol = 1e-6, max_iter = 100) {
  beta <- beta_init
  for (i in 1:max_iter) {
    eta <- X %*% beta
    mu <- exp(eta)
    
    # Gradient (score function)
    score <- t(X) %*% (y - mu)
    
    # Hessian (second derivative)
    W <- diag(as.vector(mu))
    hessian <- - t(X) %*% W %*% X
    
    # Update beta
    beta_new <- beta - solve(hessian) %*% score
    
    # Cek konvergensi
    if (sum(abs(beta_new - beta)) < tol) {
      cat("Konvergen di iterasi ke-", i, "\n")
      break
    }
    beta <- beta_new
  }
  return(beta)
}

# Jalankan algoritma
beta_init <- matrix(0, nrow = ncol(X), ncol = 1)
beta_est <- newton_raphson_poisson(X, y, beta_init)

## Konvergen di iterasi ke- 6

# Tampilkan hasil
colnames(beta_est) <- "Estimate"
rownames(beta_est) <- c("Intercept", "UangSaku", "Pengetahuan")
beta_est

##               Estimate
## Intercept   -0.8313990
## UangSaku     0.4031486
## Pengetahuan -0.1742036

1. Intercept (-0.8314) μ=exp(−0.8314)≈0.435

Artinya, jika UangSaku = 0 dan Pengetahuan = 0, maka jumlah merokok yang diprediksi adalah sekitar 0.435 orang.

2. UangSaku (+0.4031)

Koefisien ini menunjukkan bahwa semakin besar uang saku, semakin besar kemungkinan orang tersebut merokok.

Dalam skala log: log(μ)↑0.4031⇒μ×exp(0.4031)≈μ×1.497 Artinya, setiap tambahan 1 unit uang saku, prediksi jumlah merokok naik sekitar 49.7% (dengan asumsi pengetahuan tetap).

3. Pengetahuan (–0.1742) Efek negatif → semakin tinggi pengetahuan, jumlah merokok diprediksi menurun.

Penurunan sekitar:

exp(−0.1742)≈0.84⇒turun sekitar 16%

10.27 Uji Wald

# Desain matriks X (dengan intercept)
X <- as.matrix(cbind(1, data$UangSaku, data$Pengetahuan))
y <- data$Merokok

# Hitung mu dan matriks W
eta <- X %*% beta_est
mu <- exp(eta)
W <- diag(as.vector(mu))

# Fisher Information (negatif Hessian)
fisher_info <- t(X) %*% W %*% X
vcov_matrix <- solve(fisher_info)

# Standard Error
se_beta <- sqrt(diag(vcov_matrix))

# Z-value dan P-value
z_values <- as.vector(beta_est) / se_beta
p_values <- 2 * (1 - pnorm(abs(z_values)))

# Tabel hasil Wald test
wald_results <- data.frame(
  Estimate = round(beta_est, 4),
  Std_Error = round(se_beta, 4),
  Z_value = round(z_values, 4),
  P_value = round(p_values, 4),
  row.names = c("Intercept", "UangSaku", "Pengetahuan")
)

wald_results

Baik UangSaku maupun Pengetahuan belum menunjukkan pengaruh signifikan terhadap merokok pada tingkat signifikansi 5%.

10.28 Evaluasi Kebaikan Model AIC dan BIC

# Hitung AIC
aic_val <- AIC(model_poisson)

# Hitung BIC
bic_val <- BIC(model_poisson)

# Tampilkan hasil
cat("AIC:", aic_val, "\n")

## AIC: 220.969

cat("BIC:", bic_val, "\n")

## BIC: 229.5946

10.29 Iteratively Reweighted Least Squares (IRLS)

# Data
y <- data$Merokok
X <- model.matrix(~ UangSaku + Pengetahuan, data = data)

# Inisialisasi beta
beta <- rep(0, ncol(X))

# Fungsi log-link Poisson
g_inv <- function(eta) exp(eta)

# Iterasi IRLS
max_iter <- 100
tol <- 1e-6

for (i in 1:max_iter) {
  eta <- X %*% beta
  mu <- g_inv(eta)
  W <- diag(as.vector(mu))  # Bobot
  z <- eta + (y - mu) / mu  # Working response
  
  # Update beta dengan solusi weighted least squares
  beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
  
  # Cek konvergensi
  if (sum(abs(beta_new - beta)) < tol) {
    cat("Konvergen di iterasi ke-", i, "\n")
    break
  }
  
  beta <- beta_new
}

## Konvergen di iterasi ke- 6

# Tampilkan hasil
rownames(beta) <- colnames(X)
colnames(beta) <- "Estimate"
beta

##               Estimate
## (Intercept) -0.8313990
## UangSaku     0.4031486
## Pengetahuan -0.1742036

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

Regresi logistik adalah metode statistik yang digunakan untuk memodelkan hubungan antara variabel respons biner dengan satu atau lebih variabel prediktor yang dapat memiliki skala pengukuran berbeda.

• 1. Prediktor Nominal

Variabel nominal adalah variabel kategori tanpa urutan logis antar kategori (misal: jenis kelamin, status pernikahan, wilayah). Dalam regresi logistik:

• Dummy Coding: Setiap kategori (kecuali satu sebagai referensi) diubah menjadi variabel dummy (0/1).

• Contoh: Jenis kelamin (Female = 1, Male = 0).

• Catatan: Jangan mengasumsikan urutan atau jarak antar kategori.

• 2. Prediktor Ordinal

Variabel ordinal memiliki kategori berurutan namun jarak antar kategori belum tentu sama (misal: tingkat pendidikan, skala Likert).

1. Dummy Coding

• Perlakukan seperti nominal dengan membuat dummy untuk tiap kategori.
• Kelebihan: Tidak mengasumsikan linearitas.
• Kekurangan: Tidak memanfaatkan informasi urutan.

2. Skor Ordinal

• Beri skor numerik (misal: SMA = 1, Sarjana = 2, Master = 3, Doktor = 4).

• Asumsi: Efek linear pada log-odds

• 3. Prediktor Rasio

Variabel rasio adalah variabel numerik kontinu dengan nol absolut (misal: pendapatan, usia).

• Penggunaan: Dapat langsung dimasukkan ke model tanpa transformasi khusus.
• Penskalaan: Bisa membagi dengan standar deviasi untuk interpretasi odds ratio per SD.
• Pemeriksaan Linearitas: Jika hubungan tidak linear, gunakan spline atau polynomial.

11. 1 Simulasi Data Rill

Data diambil dari penelitian:

Rahmadeni & Sarah Puspita (2023). “Aplikasi Regresi Logistik Ordinal dalam Pemodelan Status Gizi Balita” Puskesmas Limapuluh Kota Pekanbaru, 2018. Variabel:

Respon ordinal: Status Gizi Balita (kategori: gizi buruk, gizi kurang, gizi baik, gizi lebih)

Prediktor nominal: Jenis Kelamin (Laki-laki, Perempuan)

Prediktor ordinal: Usia Balita (kelompok umur)

Prediktor rasio: Berat Badan (kg)

library(MASS)

## 
## Attaching package: 'MASS'

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

library(ordinal)

## Warning: package 'ordinal' was built under R version 4.3.3

## 
## Attaching package: 'ordinal'

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

# Contoh data simulasi berdasarkan data jurnal
set.seed(123)
data <- data.frame(
  StatusGizi = factor(sample(c("Buruk", "Kurang", "Baik", "Lebih"), 190, replace = TRUE), ordered = TRUE),
  JenisKelamin = factor(sample(c("Laki-laki", "Perempuan"), 190, replace = TRUE)),
  UsiaBalita = ordered(sample(c("0-12", "13-24", "25-36", "37-48", "49-60"), 190, replace = TRUE)),
  BeratBadan = runif(190, 5, 20)  # berat badan dalam kg, skala rasio
)

# Fit model regresi logistik ordinal
model <- clm(StatusGizi ~ JenisKelamin + UsiaBalita + BeratBadan, data = data, link = "logit")

# Ringkasan hasil model
summary(model)

## formula: StatusGizi ~ JenisKelamin + UsiaBalita + BeratBadan
## data:    data
## 
##  link  threshold nobs logLik  AIC    niter max.grad cond.H 
##  logit flexible  190  -258.90 535.79 4(0)  3.27e-10 5.7e+03
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)
## JenisKelaminPerempuan -0.10988    0.26521  -0.414    0.679
## UsiaBalita.L          -0.22848    0.32691  -0.699    0.485
## UsiaBalita.Q          -0.43909    0.30760  -1.427    0.153
## UsiaBalita.C          -0.02148    0.29266  -0.073    0.942
## UsiaBalita^4           0.12511    0.27405   0.457    0.648
## BeratBadan             0.01791    0.03152   0.568    0.570
## 
## Threshold coefficients:
##              Estimate Std. Error z value
## Baik|Buruk    -0.7164     0.4383  -1.635
## Buruk|Kurang   0.2615     0.4358   0.600
## Kurang|Lebih   1.6439     0.4533   3.627

Interpretasi : Dataset berisi 190 observasi dengan variabel Status Gizi (ordinal), Jenis Kelamin (nominal), Usia Balita (ordinal), dan Berat Badan (rasio).

11.2 Eksplorasi Data

# Muat paket magrittr agar operator %>% tersedia
library(magrittr)
library(dplyr)
# Ringkasan distribusi variabel outcome (Status Gizi)
data %>% count(StatusGizi) %>% mutate(prop = n / sum(n))

# Rata-rata Berat Badan berdasarkan Status Gizi
data %>%
  group_by(StatusGizi) %>%
  summarise(rata_berat = mean(BeratBadan),
            sd_berat = sd(BeratBadan),
            n = n())

# Distribusi Status Gizi berdasarkan Jenis Kelamin
data %>%
  group_by(JenisKelamin, StatusGizi) %>%
  summarise(jumlah = n()) %>%
  mutate(prop = jumlah / sum(jumlah))

## `summarise()` has grouped output by 'JenisKelamin'. You can override using the
## `.groups` argument.

# Distribusi Status Gizi berdasarkan Usia Balita
data %>%
  group_by(UsiaBalita, StatusGizi) %>%
  summarise(jumlah = n()) %>%
  mutate(prop = jumlah / sum(jumlah))

## `summarise()` has grouped output by 'UsiaBalita'. You can override using the
## `.groups` argument.

Interpretasi :

Distribusi status gizi relatif merata, dengan proporsi terbesar pada kategori Baik dan Kurang (~29%), diikuti Buruk (~23%) dan Lebih (~19%).

Rata-rata berat badan balita berkisar antara 11.45 kg (kategori Buruk) hingga 12.63 kg (kategori Kurang), dengan variasi (SD) sekitar 4 kg, menunjukkan perbedaan berat badan yang moderat antar kategori status gizi.

Distribusi status gizi menurut jenis kelamin menunjukkan laki-laki sedikit lebih banyak pada kategori Lebih (23%) dibanding perempuan (15%), sedangkan perempuan sedikit lebih banyak pada kategori Kurang (31%) dibanding laki-laki (27%).

Distribusi status gizi menurut kelompok usia balita relatif konsisten, dengan sedikit penurunan proporsi kategori Lebih pada kelompok usia 0-12 bulan (13%) dibanding 13-24 bulan (17%).

11.3 Perlakuan Variabel Ordinal

# Mengubah StatusGizi menjadi faktor nominal (dummy coding)
data_nominal <- data %>%
  mutate(
    StatusGizi = factor(StatusGizi, levels = c("Buruk", "Kurang", "Baik", "Lebih")),
    JenisKelamin = factor(JenisKelamin),
    UsiaBalita = factor(UsiaBalita)  # ubah ordinal ke nominal untuk model ini
  )

# Membuat model regresi logistik biner contoh: prediksi StatusGizi "Baik" vs lainnya
# (karena regresi logistik biner, kita buat variabel target biner)
data_nominal <- data_nominal %>%
  mutate(StatusBaik = ifelse(StatusGizi == "Baik", 1, 0))

model_nominal <- glm(StatusBaik ~ JenisKelamin + UsiaBalita + BeratBadan,
                     data = data_nominal,
                     family = binomial)

# Menampilkan ringkasan model
summary(model_nominal)

## 
## Call:
## glm(formula = StatusBaik ~ JenisKelamin + UsiaBalita + BeratBadan, 
##     family = binomial, data = data_nominal)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)
## (Intercept)           -0.698421   0.521764  -1.339    0.181
## JenisKelaminPerempuan -0.202918   0.327576  -0.619    0.536
## UsiaBalita.L           0.266113   0.387296   0.687    0.492
## UsiaBalita.Q           0.281335   0.369226   0.762    0.446
## UsiaBalita.C          -0.051911   0.364452  -0.142    0.887
## UsiaBalita^4           0.065676   0.341006   0.193    0.847
## BeratBadan            -0.007817   0.038249  -0.204    0.838
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 228.64  on 189  degrees of freedom
## Residual deviance: 227.01  on 183  degrees of freedom
## AIC: 241.01
## 
## Number of Fisher Scoring iterations: 4

Interpretasi

1. Model dan Variabel Model memprediksi probabilitas balita memiliki status gizi “Baik” (StatusBaik = 1) berdasarkan prediktor:

JenisKelamin (Perempuan dibanding baseline Laki-laki)

UsiaBalita (variabel ordinal yang dikodekan polinomial: linear (L), kuadratik (Q), kubik (C), orde 4)

BeratBadan (variabel kontinu rasio)

Semua prediktor tidak signifikan pada level 0.05.

Tidak ada variabel prediktor yang berpengaruh signifikan terhadap probabilitas balita memiliki status gizi “Baik” dalam model ini.

Model ini belum mampu menjelaskan variasi status gizi “Baik” secara memadai.

11.4 Perlakuan sebagai Rasio (Numeric Rank)

# Mengubah StatusGizi menjadi variabel biner untuk contoh regresi logistik
data <- data %>%
  mutate(
    StatusBaik = ifelse(StatusGizi == "Baik", 1, 0),
    UsiaBalita_num = as.numeric(UsiaBalita),  # ordinal diubah ke numeric rank
    JenisKelamin = factor(JenisKelamin)
  )

# Membuat model regresi logistik dengan UsiaBalita sebagai numeric rank
model_rasio <- glm(StatusBaik ~ JenisKelamin + UsiaBalita_num + BeratBadan,
                   data = data,
                   family = binomial)

# Menampilkan ringkasan model
summary(model_rasio)

## 
## Call:
## glm(formula = StatusBaik ~ JenisKelamin + UsiaBalita_num + BeratBadan, 
##     family = binomial, data = data)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)
## (Intercept)           -0.93734    0.68549  -1.367    0.171
## JenisKelaminPerempuan -0.20032    0.32569  -0.615    0.539
## UsiaBalita_num         0.08662    0.12441   0.696    0.486
## BeratBadan            -0.01109    0.03770  -0.294    0.769
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 228.64  on 189  degrees of freedom
## Residual deviance: 227.61  on 186  degrees of freedom
## AIC: 235.61
## 
## Number of Fisher Scoring iterations: 4

Interpretasi :

Semua prediktor tidak signifikan pada level 0.05.

Penurunan deviance sangat kecil, menunjukkan model dengan prediktor hanya sedikit lebih baik dari model intercept saja.

Tidak ada variabel prediktor yang berpengaruh signifikan terhadap probabilitas balita memiliki status gizi “Baik” dalam model ini.

Model ini belum mampu menjelaskan variasi status gizi “Baik” secara memadai.

11.5 Goodness-of-Fit Model

# Model null (hanya intercept)
nullmod <- glm(StatusBaik ~ 1, data = data, family = binomial)

# Menghitung McFadden R2
r2_nominal <- 1 - (as.numeric(logLik(model_nominal)) / as.numeric(logLik(nullmod)))
r2_numeric <- 1 - (as.numeric(logLik(model_rasio)) / as.numeric(logLik(nullmod)))

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

## $McFadden_R2_Nominal
## [1] 0.007132963
## 
## $McFadden_R2_Numeric
## [1] 0.004488465

Interpretasi:

Nilai McFadden R² yang sangat kecil (mendekati nol) menunjukkan bahwa kedua model hanya sedikit lebih baik dibanding model null (intercept saja) dalam menjelaskan variasi data.

Model dengan prediktor nominal sedikit lebih baik daripada model dengan prediktor numeric rank, tetapi perbedaannya sangat kecil dan keduanya masih sangat rendah.

11.6 Plot untuk model Nominal

library(ggplot2)
#Menambahkan prediksi probabilitas pada data nominal
data_nominal <- data %>%
  mutate(predicted_nominal = predict(model_nominal, type = "response"))

# Menambahkan prediksi probabilitas pada data numeric rank
data_rasio <- data %>%
  mutate(predicted_rasio = predict(model_rasio, type = "response"))

# Visualisasi prediksi probabilitas model nominal berdasarkan BeratBadan dan JenisKelamin
ggplot(data_nominal, aes(x = BeratBadan, y = predicted_nominal, color = JenisKelamin)) +
  geom_point(alpha = 0.6) +
  labs(
    title = "Prediksi Probabilitas Status Gizi 'Baik' (Model Nominal)",
    x = "Berat Badan (kg)",
    y = "Prediksi Probabilitas"
  ) +
  theme_minimal()

11.6 Plot untuk model Rasio

# Menambahkan prediksi probabilitas ke data
class(data)

## [1] "data.frame"

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.3.3

## Warning: package 'readr' was built under R version 4.3.3

## Warning: package 'forcats' was built under R version 4.3.3

## Warning: package 'lubridate' was built under R version 4.3.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.0
## ✔ readr     2.1.5     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract()   masks magrittr::extract()
## ✖ dplyr::filter()    masks stats::filter()
## ✖ dplyr::lag()       masks stats::lag()
## ✖ MASS::select()     masks dplyr::select()
## ✖ purrr::set_names() masks magrittr::set_names()
## ✖ ordinal::slice()   masks dplyr::slice()
## ✖ dplyr::summarise() masks vcdExtra::summarise()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(dplyr)
data <- data %>% mutate(predicted = predict(model_rasio, type = "response"))

# Visualisasi prediksi probabilitas berdasarkan BeratBadan dan JenisKelamin
ggplot(data, aes(x = BeratBadan, y = predicted, color = JenisKelamin)) +
  geom_point(alpha = 0.6) +
  labs(
    title = "Prediksi Probabilitas Status Gizi 'Baik' (Model Rasio)",
    x = "Berat Badan (kg)",
    y = "Prediksi Probabilitas"
  ) +
  theme_minimal()

Interpretasi:

Model Nominal: Variasi prediksi probabilitas antar individu dan kelompok tampak lebih besar, dengan rentang prediksi yang sedikit lebih lebar (hingga sekitar 0.40).

Model Rasio: Prediksi probabilitas lebih terkonsentrasi dan cenderung lebih rendah secara umum, dengan rentang prediksi maksimal sekitar 0.36.

Model nominal (dummy coding) tampak lebih sensitif terhadap perbedaan kategori prediktor, sedangkan model rasio memberikan pola yang lebih halus namun dengan variasi yang lebih kecil.

11.7 Ringkasan Koefsien Model

library(broom)
library(kableExtra)

## Warning: package 'kableExtra' was built under R version 4.3.3

## 
## Attaching package: 'kableExtra'

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

library(dplyr)
#Model Nominal
tidy(model_nominal, conf.int = TRUE) %>%
  mutate(across(where(is.numeric), round, 3)) %>%
  kable(format = "html", 
        caption = "Ringkasan Koefisien Model Nominal (tanpa Odds Ratio)",
        booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(where(is.numeric), round, 3)`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
## 
##   # Previously
##   across(a:b, mean, na.rm = TRUE)
## 
##   # Now
##   across(a:b, \(x) mean(x, na.rm = TRUE))

Ringkasan Koefisien Model Nominal (tanpa Odds Ratio)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-0.698	0.522	-1.339	0.181	-1.741	0.313
JenisKelaminPerempuan	-0.203	0.328	-0.619	0.536	-0.853	0.436
UsiaBalita.L	0.266	0.387	0.687	0.492	-0.492	1.035
UsiaBalita.Q	0.281	0.369	0.762	0.446	-0.452	1.002
UsiaBalita.C	-0.052	0.364	-0.142	0.887	-0.774	0.663
UsiaBalita^4	0.066	0.341	0.193	0.847	-0.608	0.736
BeratBadan	-0.008	0.038	-0.204	0.838	-0.083	0.067

#Model Numeric
tidy(model_rasio, conf.int = TRUE) %>%
  mutate(across(where(is.numeric), round, 3)) %>%
  kable(format = "html", 
        caption = "Ringkasan Koefisien Model Numeric (tanpa Odds Ratio)",
        booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Ringkasan Koefisien Model Numeric (tanpa Odds Ratio)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-0.937	0.685	-1.367	0.171	-2.312	0.389
JenisKelaminPerempuan	-0.200	0.326	-0.615	0.539	-0.846	0.435
UsiaBalita_num	0.087	0.124	0.696	0.486	-0.157	0.333
BeratBadan	-0.011	0.038	-0.294	0.769	-0.085	0.063

12. Pemilihan Model Regresi Logistik dan Evaluasi

12.1 Membangun Model Regresi Logistik: Pendekatan Confrmatory dan Exploratory

Analisis pendekatan konfirmatori dan eksploratori dalam penelitian kuantitatif memiliki karakteristik dan tujuan berbeda. Pendekatan konfirmatori digunakan untuk menguji teori yang telah dirumuskan sebelumnya, sementara pendekatan eksploratori bertujuan mengeksplorasi hubungan antar variabel tanpa asumsi awal.

1. Pendekatan Konfirmatori (Confirmatory Approach)

Ciri utama:

Dibangun berdasarkan teori atau temuan penelitian sebelumnya untuk menguji validitas konstruk.

Menggunakan model persamaan struktural (SEM) dengan analisis faktor konfirmatori (CFA) untuk memverifikasi hubungan antara variabel laten dan indikator.

Fokus pada pengujian signifikansi efek melalui perbandingan model dengan uji rasio kemungkinan (Likelihood Ratio Test).

Contoh implementasi: Penelitian validasi skala psikologi menggunakan CFA untuk mengonfirmasi struktur faktor yang telah dihipotesiskan. Model dianggap fit jika memenuhi kriteria RMSEA <0.08 dan CFI >0.90.

2. Pendekatan Eksploratori (Exploratory Approach) Karakteristik khusus:

Digunakan ketika belum ada teori yang mapan untuk menjelaskan fenomena.

Mengandalkan analisis faktor eksploratori (EFA) untuk mengidentifikasi pola hubungan antar variabel indikator.

Melibatkan teknik seleksi variabel bertahap (forward/backward selection) berdasarkan kriteria statistik seperti AIC.

3. Aplikasi: Studi adaptasi instrumen psikometri menggunakan EFA untuk mengeksplorasi dimensi laten yang muncul secara alami dari data, sebelum dikonfirmasi dengan CFA

12.2 Simulasi Data

Data sekunder dari Dinas Pendidikan Provinsi Maluku, meliputi 21 SMA di Kota Ambon dengan variabel:

Respon: Peringkat akreditasi SMA (ordinal: C = cukup, B = baik, A = amat baik)

Prediktor: Status sekolah (negeri/swasta), lama berdiri, jumlah siswa, jumlah guru, status tanah bangunan, nilai rata-rata UN

12.3 Pemilihan Model

set.seed(123)
df_maluku <- data.frame(
  LamaBerdiri = round(runif(100, 1, 50)),
  JumlahSiswa = round(runif(100, 100, 1000)),
  JumlahGuru = round(runif(100, 10, 100)),
  NilaiUN = runif(100, 50, 100),
  StatusTanah = sample(0:1, 100, replace = TRUE),  # 0 = sewa, 1 = milik sendiri
  StatusSekolah = sample(0:1, 100, replace = TRUE) # 0 = negeri, 1 = swasta
)

library(lavaan)

## Warning: package 'lavaan' was built under R version 4.3.3

## This is lavaan 0.6-19
## lavaan is FREE software! Please report any bugs.

model_cfa_maluku <- '
  F1 =~ LamaBerdiri + JumlahSiswa + JumlahGuru
  F2 =~ NilaiUN + StatusTanah + StatusSekolah
  F1 ~~ F2
'

fit_cfa_maluku <- cfa(model_cfa_maluku, data = df_maluku)

## Warning: lavaan->lav_data_full():  
##    some observed variances are (at least) a factor 1000 times larger than 
##    others; use varTable(fit) to investigate

## Warning: lavaan->lav_lavaan_step11_estoptim():  
##    Model estimation FAILED! Returning starting values.

summary(fit_cfa_maluku, standardized = TRUE, fit.measures = TRUE)

## Warning: lavaan->lav_object_summary():  
##    fit measures not available if model did not converge

## lavaan 0.6-19 did NOT end normally after 552 iterations
## ** WARNING ** Estimates below are most likely unreliable
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        13
## 
##   Number of observations                           100
## 
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate    Std.Err  z-value  P(>|z|)   Std.lv    Std.all
##   F1 =~                                                                     
##     LamaBerdiri         1.000                                 0.009    0.001
##     JumlahSiswa     28771.581       NA                      265.713    1.125
##     JumlahGuru       -260.807       NA                       -2.409   -0.091
##   F2 =~                                                                     
##     NilaiUN             1.000                                 0.277    0.019
##     StatusTanah        -7.443       NA                       -2.063   -4.125
##     StatusSekolah      -0.044       NA                       -0.012   -0.025
## 
## Covariances:
##                    Estimate    Std.Err  z-value  P(>|z|)   Std.lv    Std.all
##   F1 ~~                                                                     
##     F2                 -0.000       NA                       -0.004   -0.004
## 
## Variances:
##                    Estimate    Std.Err  z-value  P(>|z|)   Std.lv    Std.all
##    .LamaBerdiri       194.349       NA                      194.349    1.000
##    .JumlahSiswa    -14801.599       NA                   -14801.599   -0.265
##    .JumlahGuru        688.446       NA                      688.446    0.992
##    .NilaiUN           213.646       NA                      213.646    1.000
##    .StatusTanah        -4.004       NA                       -4.004  -16.015
##    .StatusSekolah       0.245       NA                        0.245    0.999
##     F1                  0.000       NA                        1.000    1.000
##     F2                  0.077       NA                        1.000    1.000

12.4 Metode Stepwise: Forward, Backward, dan Kedua Arah

Data sekunder dari Dinas Pendidikan Provinsi Maluku, meliputi 21 SMA di Kota Ambon dengan variabel:

Respon: Peringkat akreditasi SMA (ordinal: C = cukup, B = baik, A = amat baik)

Prediktor: Status sekolah (negeri/swasta), lama berdiri, jumlah siswa, jumlah guru, status tanah bangunan, nilai rata-rata UN

library(MASS)

# Misal data simulasi mirip struktur jurnal
set.seed(123)
n <- 100
StatusSekolah <- factor(sample(c("Negeri", "Swasta"), n, replace = TRUE))
LamaBerdiri <- round(runif(n, 1, 50))
JumlahSiswa <- round(runif(n, 100, 1000))
JumlahGuru <- round(runif(n, 10, 100))
StatusTanah <- factor(sample(c("Milik Sendiri", "Sewa"), n, replace = TRUE))
NilaiUN <- runif(n, 50, 100)
Akreditasi <- ordered(sample(c("C", "B", "A"), n, replace = TRUE), levels = c("C", "B", "A"))

data <- data.frame(Akreditasi, StatusSekolah, LamaBerdiri, JumlahSiswa, JumlahGuru, StatusTanah, NilaiUN)

# Model penuh
full_model <- polr(Akreditasi ~ StatusSekolah + LamaBerdiri + JumlahSiswa + JumlahGuru + StatusTanah + NilaiUN, data = data, Hess=TRUE)

# Model kosong (intercept saja)
null_model <- polr(Akreditasi ~ 1, data = data, Hess=TRUE)

# Stepwise forward
step_forward <- step(null_model, scope = formula(full_model), direction = "forward")

## Start:  AIC=221.84
## Akreditasi ~ 1
## 
##                 Df    AIC
## + JumlahGuru     1 217.41
## + StatusTanah    1 221.39
## <none>             221.84
## + NilaiUN        1 222.00
## + JumlahSiswa    1 223.63
## + StatusSekolah  1 223.68
## + LamaBerdiri    1 223.72
## 
## Step:  AIC=217.41
## Akreditasi ~ JumlahGuru
## 
##                 Df    AIC
## <none>             217.41
## + NilaiUN        1 217.42
## + StatusTanah    1 217.64
## + JumlahSiswa    1 219.13
## + LamaBerdiri    1 219.37
## + StatusSekolah  1 219.40

# Stepwise backward
step_backward <- step(full_model, direction = "backward")

## Start:  AIC=223.07
## Akreditasi ~ StatusSekolah + LamaBerdiri + JumlahSiswa + JumlahGuru + 
##     StatusTanah + NilaiUN
## 
##                 Df    AIC
## - LamaBerdiri    1 221.19
## - StatusSekolah  1 221.25
## - JumlahSiswa    1 221.43
## - NilaiUN        1 222.74
## <none>             223.07
## - StatusTanah    1 223.22
## - JumlahGuru     1 226.55
## 
## Step:  AIC=221.19
## Akreditasi ~ StatusSekolah + JumlahSiswa + JumlahGuru + StatusTanah + 
##     NilaiUN
## 
##                 Df    AIC
## - StatusSekolah  1 219.36
## - JumlahSiswa    1 219.50
## - NilaiUN        1 220.85
## <none>             221.19
## - StatusTanah    1 221.30
## - JumlahGuru     1 224.82
## 
## Step:  AIC=219.36
## Akreditasi ~ JumlahSiswa + JumlahGuru + StatusTanah + NilaiUN
## 
##               Df    AIC
## - JumlahSiswa  1 217.68
## - NilaiUN      1 218.97
## - StatusTanah  1 219.34
## <none>           219.36
## - JumlahGuru   1 223.30
## 
## Step:  AIC=217.68
## Akreditasi ~ JumlahGuru + StatusTanah + NilaiUN
## 
##               Df    AIC
## - StatusTanah  1 217.42
## - NilaiUN      1 217.64
## <none>           217.68
## - JumlahGuru   1 221.59
## 
## Step:  AIC=217.42
## Akreditasi ~ JumlahGuru + NilaiUN
## 
##              Df    AIC
## - NilaiUN     1 217.41
## <none>          217.42
## - JumlahGuru  1 222.00
## 
## Step:  AIC=217.41
## Akreditasi ~ JumlahGuru
## 
##              Df    AIC
## <none>          217.41
## - JumlahGuru  1 221.84

# Stepwise kedua arah
step_both <- step(null_model, scope = formula(full_model), direction = "both")

## Start:  AIC=221.84
## Akreditasi ~ 1
## 
##                 Df    AIC
## + JumlahGuru     1 217.41
## + StatusTanah    1 221.39
## <none>             221.84
## + NilaiUN        1 222.00
## + JumlahSiswa    1 223.63
## + StatusSekolah  1 223.68
## + LamaBerdiri    1 223.72
## 
## Step:  AIC=217.41
## Akreditasi ~ JumlahGuru
## 
##                 Df    AIC
## <none>             217.41
## + NilaiUN        1 217.42
## + StatusTanah    1 217.64
## + JumlahSiswa    1 219.13
## + LamaBerdiri    1 219.37
## + StatusSekolah  1 219.40
## - JumlahGuru     1 221.84

summary(step_both)

## Call:
## polr(formula = Akreditasi ~ JumlahGuru, data = data, Hess = TRUE)
## 
## Coefficients:
##              Value Std. Error t value
## JumlahGuru 0.01825   0.007328    2.49
## 
## Intercepts:
##     Value   Std. Error t value
## C|B -0.0587  0.4330    -0.1356
## B|A  1.5483  0.4626     3.3473
## 
## Residual Deviance: 211.4146 
## AIC: 217.4146

12.5 Evaluasi Model: ROC dan AUC

library(pROC)

## Warning: package 'pROC' was built under R version 4.3.3

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

## 
## Attaching package: 'pROC'

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

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

library(pROC)

# Simulasi data biner outcome (misal sukses akreditasi A vs bukan A)
set.seed(123)
n <- 100
StatusSekolah <- factor(sample(c("Negeri", "Swasta"), n, replace = TRUE))
LamaBerdiri <- round(runif(n, 1, 50))
JumlahSiswa <- round(runif(n, 100, 1000))
JumlahGuru <- round(runif(n, 10, 100))
StatusTanah <- factor(sample(c("Milik Sendiri", "Sewa"), n, replace = TRUE))
NilaiUN <- runif(n, 50, 100)
Akreditasi <- factor(sample(c("C", "B", "A"), n, replace = TRUE), ordered = TRUE)

# Membuat variabel biner: 1 jika Akreditasi A, 0 jika bukan
y <- ifelse(Akreditasi == "A", 1, 0)

data <- data.frame(y = factor(y), StatusSekolah, LamaBerdiri, JumlahSiswa, JumlahGuru, StatusTanah, NilaiUN)

# Model regresi logistik biner
model_logistic <- glm(y ~ StatusSekolah + LamaBerdiri + JumlahSiswa + JumlahGuru + StatusTanah + NilaiUN,
                      data = data, family = binomial)

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

# Membuat objek ROC
roc_obj <- roc(data$y, pred_prob)

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

## Setting direction: controls < cases

# Plot kurva ROC
plot(roc_obj, main = "ROC Curve Model Regresi Logistik Akreditasi SMA Maluku")

# Menghitung AUC
auc_value <- auc(roc_obj)
print(paste("AUC:", round(auc_value, 3)))

## [1] "AUC: 0.683"

12.6 Pseudo R-Squared

set.seed(123)
n <- 100
# Simulasi variabel prediktor sesuai data Maluku
StatusSekolah <- factor(sample(c("Negeri", "Swasta"), n, replace = TRUE))
LamaBerdiri <- round(runif(n, 1, 50))
JumlahSiswa <- round(runif(n, 100, 1000))
JumlahGuru <- round(runif(n, 10, 100))
StatusTanah <- factor(sample(c("Milik Sendiri", "Sewa"), n, replace = TRUE))
NilaiUN <- runif(n, 50, 100)

# Outcome biner: Akreditasi A (1) vs selainnya (0)
Akreditasi <- factor(sample(c("C", "B", "A"), n, replace = TRUE), ordered = TRUE)
y <- ifelse(Akreditasi == "A", 1, 0)

data_maluku <- data.frame(
  y = factor(y),
  StatusSekolah,
  LamaBerdiri,
  JumlahSiswa,
  JumlahGuru,
  StatusTanah,
  NilaiUN
)

# Model penuh
model_full <- glm(y ~ StatusSekolah + LamaBerdiri + JumlahSiswa + JumlahGuru + StatusTanah + NilaiUN,
                  data = data_maluku, family = binomial)

# Model nol (intercept saja)
model_null <- glm(y ~ 1, data = data_maluku, family = binomial)

# Stepwise regression (arah kedua)
model_start <- glm(y ~ 1, data = data_maluku, family = binomial)
model_stepwise <- step(model_start, scope = list(lower = model_start, upper = model_full), direction = "both", trace = 0)

# Fungsi hitung Pseudo R2
logisticPseudoR2s <- function(model) {
  dev <- model$deviance
  nullDev <- model$null.deviance
  n <- length(model$fitted.values)
  
  R_l <- 1 - dev / nullDev                       # Hosmer and Lemeshow R2 (McFadden)
  R_cs <- 1 - exp(-(nullDev - dev) / n)         # Cox and Snell R2
  R_n <- R_cs / (1 - exp(-nullDev / n))         # Nagelkerke R2
  
  cat("Pseudo R^2 for logistic regression\n")
  cat("Hosmer and Lemeshow R^2:", round(R_l, 4), "\n")
  cat("Cox and Snell R^2:", round(R_cs, 4), "\n")
  cat("Nagelkerke R^2:", round(R_n, 4), "\n")
}

# Hitung Pseudo R2 untuk model stepwise
logisticPseudoR2s(model_stepwise)

## Pseudo R^2 for logistic regression
## Hosmer and Lemeshow R^2: 0.0813 
## Cox and Snell R^2: 0.1016 
## Nagelkerke R^2: 0.1388

# Lihat ringkasan model stepwise
summary(model_stepwise)

## 
## Call:
## glm(formula = y ~ JumlahGuru + StatusTanah + NilaiUN, family = binomial, 
##     data = data_maluku)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)     -2.958783   1.379619  -2.145   0.0320 *
## JumlahGuru       0.020057   0.008598   2.333   0.0197 *
## StatusTanahSewa -0.751322   0.439795  -1.708   0.0876 .
## NilaiUN          0.021682   0.015326   1.415   0.1571  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 131.79  on 99  degrees of freedom
## Residual deviance: 121.07  on 96  degrees of freedom
## AIC: 129.07
## 
## Number of Fisher Scoring iterations: 4

12.7 Tabel Klasifkasi dan Evaluasi

# Fit model regresi logistik biner
model <- glm(y ~ StatusSekolah + LamaBerdiri + JumlahSiswa + JumlahGuru + StatusTanah + NilaiUN,
             data = data_maluku, family = binomial)

# Prediksi probabilitas
prob <- predict(model, type = "response")

# Tentukan cutoff threshold 0.5 untuk klasifikasi
pred_class <- ifelse(prob >= 0.5, 1, 0)

# Buat tabel klasifikasi (confusion matrix)
table_klasifikasi <- table(Predicted = pred_class, Actual = as.numeric(as.character(data_maluku$y)))
print(table_klasifikasi)

##          Actual
## Predicted  0  1
##         0 50 24
##         1 13 13

# Menghitung metrik evaluasi
TP <- table_klasifikasi["1", "1"]  # True Positive
TN <- table_klasifikasi["0", "0"]  # True Negative
FP <- table_klasifikasi["1", "0"]  # False Positive
FN <- table_klasifikasi["0", "1"]  # False Negative

accuracy <- (TP + TN) / sum(table_klasifikasi)
sensitivity <- TP / (TP + FN)  # Recall
specificity <- TN / (TN + FP)
precision <- TP / (TP + FP)

cat("Evaluasi Model:\n")

## Evaluasi Model:

cat(sprintf("Accuracy    : %.3f\n", accuracy))

## Accuracy    : 0.630

cat(sprintf("Sensitivity : %.3f\n", sensitivity))

## Sensitivity : 0.351

cat(sprintf("Specificity : %.3f\n", specificity))

## Specificity : 0.794

cat(sprintf("Precision   : %.3f\n", precision))

## Precision   : 0.500

Kesimpulan:

Accuracy (0.630): Model secara keseluruhan benar mengklasifikasikan 63% dari total data. Ini menunjukkan performa model sedang, tidak sangat baik tapi juga tidak buruk.

Sensitivity (Recall) (0.351): Model hanya mampu mendeteksi 35.1% dari kasus positif sebenarnya (kelas 1). Ini berarti banyak kasus positif yang terlewat (false negative cukup tinggi), sehingga model kurang sensitif dalam mengenali kelas positif.

Specificity (0.794): Model mampu mengenali 79.4% kasus negatif dengan benar (kelas 0). Ini menunjukkan model cukup baik dalam mengidentifikasi kelas negatif.

Precision (0.500): Dari semua prediksi positif yang dibuat model, 50% benar-benar positif. Ini menunjukkan tingkat kesalahan prediksi positif cukup tinggi (false positive juga signifikan).

12.8 Metode Perbandingan Model dalam Regresi Logistik

Tujuannya menyajikan cara membandingkan model regresi logistik menggunakan ukuran: - Deviance - AIC (Akaike Information Criterion) - Likelihood-Ratio serta menjelaskan prinsip Parsimony dalam pemilihan model.

# Model 1: hanya dua prediktor
model1 <- glm(y ~ LamaBerdiri + StatusSekolah, data = data_maluku, family = binomial)

# Model 2: semua prediktor
model2 <- glm(y ~ LamaBerdiri + StatusSekolah + JumlahSiswa + JumlahGuru + StatusTanah + NilaiUN,
              data = data_maluku, family = binomial)

# Ringkasan model
summary(model1)

## 
## Call:
## glm(formula = y ~ LamaBerdiri + StatusSekolah, family = binomial, 
##     data = data_maluku)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)         -0.310277   0.486539  -0.638    0.524
## LamaBerdiri         -0.009189   0.016257  -0.565    0.572
## StatusSekolahSwasta  0.039748   0.421088   0.094    0.925
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 131.79  on 99  degrees of freedom
## Residual deviance: 131.47  on 97  degrees of freedom
## AIC: 137.47
## 
## Number of Fisher Scoring iterations: 4

summary(model2)

## 
## Call:
## glm(formula = y ~ LamaBerdiri + StatusSekolah + JumlahSiswa + 
##     JumlahGuru + StatusTanah + NilaiUN, family = binomial, data = data_maluku)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)  
## (Intercept)         -2.3400869  1.6145367  -1.449   0.1472  
## LamaBerdiri         -0.0093758  0.0175548  -0.534   0.5933  
## StatusSekolahSwasta  0.0170079  0.4592096   0.037   0.9705  
## JumlahSiswa         -0.0004745  0.0008714  -0.545   0.5860  
## JumlahGuru           0.0202637  0.0087685   2.311   0.0208 *
## StatusTanahSewa     -0.7973998  0.4557165  -1.750   0.0802 .
## NilaiUN              0.0201053  0.0155111   1.296   0.1949  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 131.79  on 99  degrees of freedom
## Residual deviance: 120.54  on 93  degrees of freedom
## AIC: 134.54
## 
## Number of Fisher Scoring iterations: 4

# Perbandingan AIC
cat("AIC Model 1:", AIC(model1), "\n")

## AIC Model 1: 137.4686

cat("AIC Model 2:", AIC(model2), "\n")

## AIC Model 2: 134.5409

# Perbandingan Deviance
cat("Deviance Model 1:", deviance(model1), "\n")

## Deviance Model 1: 131.4686

cat("Deviance Model 2:", deviance(model2), "\n")

## Deviance Model 2: 120.5409

# Uji Likelihood Ratio Test (LRT) antara model 1 dan model 2
lrt <- anova(model1, model2, test = "Chisq")
print(lrt)

## Analysis of Deviance Table
## 
## Model 1: y ~ LamaBerdiri + StatusSekolah
## Model 2: y ~ LamaBerdiri + StatusSekolah + JumlahSiswa + JumlahGuru + 
##     StatusTanah + NilaiUN
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1        97     131.47                       
## 2        93     120.54  4   10.928  0.02739 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

12.10 Prinsip Parsimony

Jika penurunan AIC tidak signifkan, pilih model lebih sederhana. Parsimony mencegah overftting.

1. Akaike Information Criterion (AIC)

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

AIC adalah kriteria untuk memilih model terbaik dengan mempertimbangkan keseimbangan antara goodness-of-fit dan kompleksitas model (jumlah parameter \(k\)). Model dengan nilai AIC lebih kecil dianggap lebih baik karena memberikan fit yang baik tanpa overfitting.

2. Deviance

\[ D = -2 \left[ \log L(\text{model}) - \log L(\text{model saturasi}) \right] \]

Deviance mengukur seberapa jauh model yang diestimasi dari model saturasi (model sempurna). Nilai deviance yang lebih kecil menunjukkan model yang lebih baik dalam menjelaskan data.

2. Likelihood Ratio Test (LRT)

\[ G^2 = -2 \left( \log L_0 - \log L_1 \right) \]

LRT digunakan untuk membandingkan dua model bersarang: model sederhana (\(L_0\)) dan model kompleks (\(L_1\)). Statistik \(G^2\) mengikuti distribusi chi-square dengan derajat kebebasan sesuai selisih jumlah parameter. Jika p-value kecil, model kompleks secara signifikan lebih baik.

12.11 Evaluasi Tabel Klasifkasi dan Akurasi Model

# Outcome biner: kasus perceraian terdeteksi (1) atau tidak (0)

jumlah_pulau <- rpois(n, lambda = 3)
jarak_ibu_kota <- runif(n, 10, 200)  # dalam km
luas_kabupaten <- runif(n, 500, 5000) # dalam km2

lin_pred <- -2 + 0.5 * jumlah_pulau - 0.01 * jarak_ibu_kota + 0.0002 * luas_kabupaten
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)

df <- data.frame(
  y = factor(y),
  jumlah_pulau,
  jarak_ibu_kota,
  luas_kabupaten
)

# Bangun model regresi logistik
model <- glm(y ~ jumlah_pulau + jarak_ibu_kota + luas_kabupaten, data = df, family = binomial)

# Prediksi probabilitas
prob <- predict(model, type = "response")

# Tentukan cutoff threshold 0.5 untuk klasifikasi
pred_class <- ifelse(prob >= 0.5, 1, 0)

# Buat tabel klasifikasi (confusion matrix)
conf_matrix <- table(Predicted = pred_class, Actual = as.numeric(as.character(df$y)))
print(conf_matrix)

##          Actual
## Predicted  0  1
##         0 73 18
##         1  4  5

# Hitung akurasi
accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix)
cat(sprintf("Akurasi Model: %.3f\n", accuracy))

## Akurasi Model: 0.780

12.2 Sensitivitas dan Spesifsitas

1. Sensitivitas (True Positive Rate)

Sensitivitas mengukur kemampuan model untuk mendeteksi dengan benar kasus positif.
Rumusnya:

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

di mana:
- \(TP\) = True Positive (kasus positif yang benar diprediksi positif)
- \(FN\) = False Negative (kasus positif yang salah diprediksi negatif)

Sensitivitas tinggi penting untuk meminimalkan kasus positif yang terlewat, misalnya dalam diagnosis penyakit.

2. Spesifisitas (True Negative Rate)

Spesifisitas mengukur kemampuan model untuk mendeteksi dengan benar kasus negatif.
Rumusnya:

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

di mana:
- \(TN\) = True Negative (kasus negatif yang benar diprediksi negatif)
- \(FP\) = False Positive (kasus negatif yang salah diprediksi positif)

Spesifisitas tinggi penting untuk menghindari kesalahan mengklasifikasikan kasus negatif sebagai positif.

Sumber data : BPS Maluku 2019 Data perceraian di Maluku dapat diakses dari Badan Pusat Statistik Provinsi Maluku

set.seed(2025)
n <- 150

# Simulasi variabel berdasarkan faktor yang dilaporkan berpengaruh
jumlah_perceraian <- rpois(n, lambda = 5)           # Jumlah perceraian per kabupaten/kota (indikator outcome)
jumlah_pulau <- rpois(n, lambda = 3)                # Jumlah pulau di kabupaten/kota
jarak_ibu_kota <- runif(n, 10, 200)                  # Jarak ke ibu kota kabupaten (km)
luas_kabupaten <- runif(n, 500, 5000)                # Luas wilayah kabupaten (km2)

# Outcome biner: kasus perceraian tinggi (1) jika jumlah perceraian > median, 0 jika tidak
median_perceraian <- median(jumlah_perceraian)
y <- ifelse(jumlah_perceraian > median_perceraian, 1, 0)

df_maluku <- data.frame(
  y = factor(y),
  jumlah_pulau,
  jarak_ibu_kota,
  luas_kabupaten
)

# Model regresi logistik sederhana berdasarkan variabel yang tersedia
model <- glm(y ~ jumlah_pulau + jarak_ibu_kota + luas_kabupaten, data = df_maluku, family = binomial)

# Ringkasan model
summary(model)

## 
## Call:
## glm(formula = y ~ jumlah_pulau + jarak_ibu_kota + luas_kabupaten, 
##     family = binomial, data = df_maluku)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -2.435e-01  5.626e-01  -0.433    0.665
## jumlah_pulau    1.368e-01  9.764e-02   1.401    0.161
## jarak_ibu_kota -1.590e-03  3.027e-03  -0.525    0.599
## luas_kabupaten -8.323e-05  1.293e-04  -0.644    0.520
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 205.78  on 149  degrees of freedom
## Residual deviance: 203.15  on 146  degrees of freedom
## AIC: 211.15
## 
## Number of Fisher Scoring iterations: 4

# Prediksi probabilitas
prob <- predict(model, type = "response")

# Klasifikasi dengan threshold 0.5
pred_class <- ifelse(prob >= 0.5, 1, 0)

# Tabel klasifikasi (confusion matrix)
conf_matrix <- table(Predicted = pred_class, Actual = as.numeric(as.character(df_maluku$y)))
print(conf_matrix)

##          Actual
## Predicted  0  1
##         0 69 50
##         1 15 16

# Hitung sensitivitas dan spesifisitas
TP <- conf_matrix["1", "1"]
TN <- conf_matrix["0", "0"]
FP <- conf_matrix["1", "0"]
FN <- conf_matrix["0", "1"]

sensitivitas <- TP / (TP + FN)
spesifisitas <- TN / (TN + FP)

cat(sprintf("Sensitivitas: %.3f\n", sensitivitas))

## Sensitivitas: 0.242

cat(sprintf("Spesifisitas: %.3f\n", spesifisitas))

## Spesifisitas: 0.821

12.3 Detail ROC Penjelasan Kurva ROC (Receiver Operating Characteristic)

1. Definisi Kurva ROC (Receiver Operating Characteristic)

Kurva ROC adalah grafik yang digunakan untuk mengevaluasi performa model klasifikasi biner dengan memplot:

• Sumbu Y: Sensitivitas (True Positive Rate, TPR) = \(\frac{TP}{TP + FN}\)
  
• Sumbu X: 1 - Spesifisitas (False Positive Rate, FPR) = \(\frac{FP}{FP + TN}\)

Garis diagonal dari kiri bawah ke kanan atas menunjukkan performa model acak (random guess). Kurva yang mendekati pojok kiri atas menunjukkan performa klasifikasi yang lebih baik.

2. Pergerakan Kurva ROC Berdasarkan Threshold

• Saat threshold menurun, model mengklasifikasikan lebih banyak pengamatan sebagai positif:
  • Sensitivitas naik
    
  • Spesifisitas turun
• Saat threshold naik, model menjadi lebih konservatif:
  • Sensitivitas turun
    
  • Spesifisitas naik

3. Kurva ROC Ideal dan Interpretasi AUC

Kurva ROC ideal memiliki bentuk naik tajam ke sensitivitas = 1, lalu bergerak horizontal ke FPR = 1.
Area Under the Curve (AUC) mengukur luas area di bawah kurva ROC:
- AUC = 0.5 → model setara tebakan acak
- AUC > 0.7 → model cukup baik
- AUC > 0.9 → model sangat baik

AUC juga dikenal sebagai concordance index, yaitu probabilitas model memberikan skor lebih tinggi untuk kasus positif dibanding negatif.

4. Kegunaan Kurva ROC

• Membandingkan performa beberapa model klasifikasi
  
• Memilih threshold optimal berdasarkan trade-off sensitivitas dan spesifisitas
  
• Sangat berguna pada dataset dengan kelas tidak seimbang

5. Visualisasi data dalam R

# Cek distribusi kelas
table(df$y)

## 
##  0  1 
## 77 23

# Jika salah satu kelas kosong, coba buat ulang data simulasi dengan parameter yang menyeimbangkan kelas
set.seed(123)
n <- 150
jumlah_pulau <- rpois(n, lambda = 3)
jarak_ibu_kota <- runif(n, 10, 200)
luas_kabupaten <- runif(n, 500, 5000)

# Modifikasi linear predictor agar probabilitas tidak terlalu ekstrim
lin_pred <- 0.5 + 0.2 * jumlah_pulau - 0.01 * jarak_ibu_kota + 0.0005 * luas_kabupaten
p <- 1 / (1 + exp(-lin_pred))

# Outcome biner
y_binary <- rbinom(n, 1, p)

df <- data.frame(y = factor(y_binary, levels = c(0,1)), jumlah_pulau, jarak_ibu_kota, luas_kabupaten)

# Cek distribusi kelas lagi
table(df$y)

## 
##   0   1 
##  35 115

# Fit model regresi logistik
model <- glm(y ~ jumlah_pulau + jarak_ibu_kota + luas_kabupaten, data = df, family = binomial)
prob <- predict(model, type = "response")

# Hitung ROC dan AUC
roc_obj <- roc(response = df$y, predictor = prob, levels = c("0", "1"))

## Setting direction: controls < cases

plot(roc_obj, col = "blue", lwd = 2, main = "Kurva ROC Model Regresi Logistik Perceraian Maluku")
abline(a=0, b=1, lty=2, col="gray")

12.4 Simulasi Pemilihan Threshold Optimal

set.seed(2025)
n <- 150

# Simulasi variabel prediktor mirip data perceraian Maluku
jumlah_pulau <- rpois(n, lambda = 3)
jarak_ibu_kota <- runif(n, 10, 200)
luas_kabupaten <- runif(n, 500, 5000)

# Linear predictor dengan koefisien simulasi agar probabilitas seimbang
lin_pred <- -0.2 + 1.1 * (jumlah_pulau / max(jumlah_pulau)) - 0.9 * (jarak_ibu_kota / max(jarak_ibu_kota))
p <- 1 / (1 + exp(-lin_pred))

# Outcome biner: kasus perceraian (1) atau tidak (0)
y <- rbinom(n, 1, p)

data <- data.frame(y = factor(y), jumlah_pulau, jarak_ibu_kota, luas_kabupaten)

# Fit model regresi logistik dengan dua prediktor utama
model <- glm(y ~ jumlah_pulau + jarak_ibu_kota, data = data, family = binomial)

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

# Thresholds yang akan diuji
thresholds <- seq(0.1, 0.9, by = 0.05)

# Hitung Sensitivitas untuk tiap threshold
Sensitivity <- sapply(thresholds, function(t) {
  pred_class <- ifelse(pred >= t, 1, 0)
  cm <- table(Pred = pred_class, Obs = as.numeric(as.character(data$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)
})

# Hitung Spesifisitas untuk tiap threshold
Specificity <- sapply(thresholds, function(t) {
  pred_class <- ifelse(pred >= t, 1, 0)
  cm <- table(Pred = pred_class, Obs = as.numeric(as.character(data$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)
})

# Gabungkan hasil ke data frame
results <- data.frame(Threshold = thresholds, Sensitivity = Sensitivity, Specificity = Specificity)

print(results)

##    Threshold Sensitivity Specificity
## 1       0.10  1.00000000      0.0000
## 2       0.15  1.00000000      0.0000
## 3       0.20  1.00000000      0.0000
## 4       0.25  1.00000000      0.0000
## 5       0.30  1.00000000      0.0375
## 6       0.35  0.91428571      0.2000
## 7       0.40  0.77142857      0.3875
## 8       0.45  0.62857143      0.5375
## 9       0.50  0.42857143      0.6875
## 10      0.55  0.28571429      0.8375
## 11      0.60  0.17142857      0.9000
## 12      0.65  0.08571429      0.9750
## 13      0.70  0.04285714      1.0000
## 14      0.75  0.01428571      1.0000
## 15      0.80  0.00000000      1.0000
## 16      0.85  0.00000000      1.0000
## 17      0.90  0.00000000      1.0000

12.5 Precision-Recall Curve (PR Curve)

1. Definisi Precision-Recall Curve

Precision-Recall Curve (PR Curve) adalah alat evaluasi performa model klasifikasi, khususnya berguna untuk data dengan ketidakseimbangan kelas (class imbalance).

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

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

2. Interpretasi PR Curve

• PR Curve menunjukkan bagaimana presisi berubah saat recall meningkat.
  
• Idealnya, presisi dan recall keduanya tinggi, namun biasanya ada trade-off antara keduanya.
  
• Model dengan performa baik memiliki PR Curve yang melengkung ke pojok kanan atas, menunjukkan presisi dan recall yang tinggi secara bersamaan.

3. Area Under Precision-Recall Curve (AUPRC)

• Luas area di bawah kurva PR (AUPRC) mendekati 1 menunjukkan model sangat baik dalam memprediksi kelas positif.
  
• Baseline AUPRC adalah prevalensi kelas positif dalam data, yang menjadi titik acuan sederhana.

4. Perbandingan PR Curve dan ROC Curve

Aspek	ROC Curve	Precision-Recall Curve
Fokus	Semua kelas	Kelas positif saja
Kekuatan	Data seimbang	Data tidak seimbang
Sumbu Y	Sensitivitas (Recall)	Precision
Sumbu X	1 - Spesifisitas (FPR)	Recall

12.6 Visualisasi PR Curve di R

# Install dan load paket PRROC jika belum ada
if(!require(PRROC)) install.packages("PRROC")

## Loading required package: PRROC

## Warning: package 'PRROC' was built under R version 4.3.3

## Loading required package: rlang

## 
## Attaching package: 'rlang'

## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice

## The following object is masked from 'package:magrittr':
## 
##     set_names

library(PRROC)

set.seed(2025)
n <- 250

# Simulasi variabel prediktor mirip data perceraian Maluku
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.4)
x3 <- rnorm(n)

# Linear predictor dengan koefisien simulasi
lin_pred <- -0.8 + 1.3 * x1 - 0.9 * x2 + 0.7 * x3
p <- 1 / (1 + exp(-lin_pred))

# Outcome biner berdasarkan probabilitas p
y <- rbinom(n, 1, p)

data <- data.frame(y = y, x1, x2, x3)

# Fit model regresi logistik
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)

# Prediksi probabilitas
prob <- predict(model, type = "response")

# Pisahkan skor probabilitas berdasarkan kelas aktual
scores_positive <- prob[data$y == 1]  # skor kelas positif
scores_negative <- prob[data$y == 0]  # skor kelas negatif

# Hitung dan plot Precision-Recall Curve
pr <- pr.curve(scores.class0 = scores_positive,
               scores.class1 = scores_negative,
               curve = TRUE)

plot(pr, main = "Precision-Recall Curve (PR Curve)")

# Tampilkan nilai AUPRC
cat(sprintf("Area Under Precision-Recall Curve (AUPRC): %.3f\n", pr$auc.integral))

## Area Under Precision-Recall Curve (AUPRC): 0.773

12.7 Pseudo R-squared pada Regresi Logistik

Dokumen ini menjelaskan dan menghitung pseudo R-squared pada regresi logistik dengan data simulasi baru, mengacu pada Applied Logistic Regression oleh Hosmer, Lemeshow & Sturdivant (2013) dan Statistical Methods for the Social Sciences oleh Agresti (2018).

• \(R^2_{CoxSnell}\)
• \(R^2_{McFadden}\)

12.8 Simulasi Data

set.seed(2025)
n <- 250
z1 <- rnorm(n)                               # Variabel numerik (misal faktor sosial-ekonomi)
z2 <- rbinom(n, 1, 0.3)                      # Variabel biner (misal status tertentu)
z3 <- rnorm(n)                               # Variabel numerik lain

# Linear predictor dengan koefisien simulasi yang mencerminkan pengaruh variabel
lin_pred <- -1.3 + 1.5 * z1 - 0.7 * z2 + 0.5 * z3

# Probabilitas outcome biner
p <- 1 / (1 + exp(-lin_pred))

# Outcome biner z (misal kasus perceraian: 1 = terjadi, 0 = tidak)
z <- rbinom(n, 1, p)

# Data frame lengkap
data <- data.frame(z = factor(z), z1, z2, z3)

# Tampilkan ringkasan data dan distribusi kelas outcome
summary(data)

##  z             z1                 z2              z3          
##  0:184   Min.   :-2.84924   Min.   :0.000   Min.   :-2.64475  
##  1: 66   1st Qu.:-0.70558   1st Qu.:0.000   1st Qu.:-0.67554  
##          Median :-0.04942   Median :0.000   Median : 0.03484  
##          Mean   :-0.01558   Mean   :0.336   Mean   : 0.05985  
##          3rd Qu.: 0.68185   3rd Qu.:1.000   3rd Qu.: 0.62409  
##          Max.   : 2.85278   Max.   :1.000   Max.   : 3.27079

table(data$z)

## 
##   0   1 
## 184  66

12.9 Model Logistik dan Null Model

model_full <- glm(z ~ z1 + z2 + z3, data = data, family = binomial) model_null <- glm(z ~ 1, data = data, family = binomial)

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

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

Dengan: - \(L_0\): likelihood model null (tanpa prediktor) - \(L_M\): likelihood model penuh

set.seed(2025)
n <- 200
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.4)
x3 <- rnorm(n)
lin_pred <- -0.7 + 1.3 * x1 - 0.9 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
data <- data.frame(y = factor(y), x1, x2, x3)

# Fit model regresi logistik
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)

# Pastikan paket pscl terpasang dan dimuat
if (!require(pscl)) install.packages("pscl")

## Loading required package: pscl

## Warning: package 'pscl' was built under R version 4.3.3

## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.

library(pscl)

# Hitung berbagai pseudo R-squared dengan pscl::pR2
pseudo_r2_pscl <- pR2(model)

## fitting null model for pseudo-r2

print(pseudo_r2_pscl)

##          llh      llhNull           G2     McFadden         r2ML         r2CU 
##  -99.7067372 -128.8592752   58.3050760    0.2262355    0.2528769    0.3491128

# Pastikan paket DescTools terpasang dan dimuat
if (!require(DescTools)) install.packages("DescTools")
library(DescTools)

# Hitung pseudo R-squared dengan DescTools::PseudoR2 (semua jenis)
pseudo_r2_desc <- PseudoR2(model, which = "all")
print(pseudo_r2_desc)

##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.2262355       0.1951938       0.2528769       0.3491128       0.2257218 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.4008910       0.2733239       0.3704361       0.2721155     207.4134743 
##             BIC          logLik         logLik0              G2 
##     220.6067438     -99.7067372    -128.8592752      58.3050760

Interpretasi :

llh (log-likelihood model): -99.71, nilai log-likelihood dari model yang di-fit.

llhNull (log-likelihood null model): -128.86, nilai log-likelihood model dengan intercept saja (tanpa prediktor).

G2 (deviance): 58.31, pengukuran perbedaan antara model penuh dan null model, semakin besar menunjukkan model penuh lebih baik.

McFadden R² (0.226): Rasio penurunan deviance, nilai sekitar 0.2–0.4 dianggap model cukup baik dalam konteks regresi logistik.

McFaddenAdj (0.195): McFadden R² yang disesuaikan untuk jumlah prediktor dan ukuran sampel.

CoxSnell R² (0.253) dan Nagelkerke R² (0.349): Ukuran pseudo R² yang mirip R² pada regresi linear, Nagelkerke sudah disesuaikan agar maksimum 1. Nilai ini menunjukkan model menjelaskan sekitar 25–35% variasi outcome.

AldrichNelson, VeallZimmermann, Efron, McKelveyZavoina, Tjur: Variasi ukuran pseudo R² lain yang memberikan perspektif berbeda tentang kekuatan model, dengan nilai antara 0.22 sampai 0.40, mengindikasikan model memiliki kekuatan prediksi sedang.

AIC (207.41) dan BIC (220.61): Ukuran kecocokan model yang memperhitungkan kompleksitas model; nilai lebih kecil lebih baik.

13. Distribusi Multinomial

Distribusi multinomial adalah generalisasi dari distribusi binomial untuk percobaan dengan lebih dari dua kategori hasil. Probabilitas untuk mendapatkan kombinasi hasil tertentu dalam \(n\) percobaan independen dengan \(k\) kategori adalah:

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

dengan syarat:

• \(n = x_1 + x_2 + \ldots + x_k\)
• \(p_1 + p_2 + \ldots + p_k = 1\)
• \(x_i\) adalah jumlah kejadian kategori ke-\(i\)
• \(p_i\) adalah peluang kategori ke-\(i\) pada satu percobaan

13.1 Contoh Kasus Nyata

Sumber : Rachmat, M. T., & Sari, R. P. (2022). Analisa perhitungan dalam penerapan konsep distribusi peluang diskrit dan statistika deskriptif. Jurnal Pendidikan Matematika: Judika Education

Sebuah klinik mengelompokkan hasil pemeriksaan pasien menjadi tiga kategori dengan peluang:

• Tidak sesuai target (\(p_1 = 0.47\))
• Sesuai target (\(p_2 = 0.13\))
• Melampaui target (\(p_3 = 0.40\))

Jika diambil sampel 6 pasien, hitung probabilitas bahwa:

• 2 pasien tidak sesuai target
• 1 pasien sesuai target
• 3 pasien melampaui target

Perhitungan Manual

\[ P = \frac{6!}{2! \times 1! \times 3!} \times (0.47)^2 \times (0.13)^1 \times (0.40)^3 \]

\[ = 60 \times 0.2209 \times 0.13 \times 0.064 \approx 0.110 \]

Jadi, probabilitasnya sekitar 11%.

di R

library(combinat)

## 
## Attaching package: 'combinat'

## The following object is masked from 'package:utils':
## 
##     combn

#Parameter
n <- 6
x <- c(2, 1, 3)
p <- c(0.47, 0.13, 0.40)

#Hitung probabilitas distribusi multinomial
prob <- dmnom(x, n, p)
print(prob)

## [1] 0.1102733

Interpretasi :

Dalam sebuah klinik, probabilitas untuk mendapatkan 6 pasien dengan distribusi hasil pemeriksaan sebagai berikut: 2 pasien tidak sesuai target, 1 pasien sesuai target, dan 3 pasien melampaui target, adalah sebesar 0,1103 atau sekitar 11%. Artinya, jika kita mengambil sampel acak sebanyak 6 pasien dari populasi dengan peluang kategori masing-masing 47% tidak sesuai target, 13% sesuai target, dan 40% melampaui target, maka kemungkinan tepat terjadi kombinasi tersebut adalah sekitar 11%

14. Multinomial Logistic Regression

Multinomial logistic regression adalah metode statistik untuk memodelkan variabel dependen kategorikal dengan lebih dari dua kategori tanpa urutan (nominal). Model ini memperluas regresi logistik biner dengan membandingkan setiap kategori dengan kategori referensi.

Model ini menggunakan fungsi prediktor linear:

\[ f(k,i) = \beta_{0,k} + \beta_{1,k} x_{1,i} + \beta_{2,k} x_{2,i} + \cdots + \beta_{M,k} x_{M,i} \]

atau dalam bentuk vektor:

\[ f(k,i) = \boldsymbol{\beta}_k \cdot \mathbf{x}_i \]

dengan \(k = 1, \ldots, K-1\) kategori (kategori ke-\(K\) menjadi referensi).

Probabilitas kategori \(k\) untuk observasi \(i\) adalah:

\[ P(Y_i = k) = \frac{\exp(f(k,i))}{1 + \sum_{j=1}^{K-1} \exp(f(j,i))}, \quad k = 1, \ldots, K-1 \]

dan untuk kategori referensi \(K\):

\[ P(Y_i = K) = \frac{1}{1 + \sum_{j=1}^{K-1} \exp(f(j,i))} \]

Log-rasio peluang antara kategori \(k\) dan referensi \(K\) dinyatakan sebagai:

\[ \log \frac{P(Y_i = k)}{P(Y_i = K)} = \boldsymbol{\beta}_k \cdot \mathbf{x}_i \]

Koefisien \(\boldsymbol{\beta}_k\) mengukur perubahan log-rasio peluang kategori \(k\) terhadap referensi untuk setiap kenaikan satu unit pada variabel prediktor.

14.1. Baseline-Category Logit Model

Baseline-category logit model adalah model regresi logistik multinomial yang digunakan untuk memodelkan variabel dependen kategorikal nominal dengan lebih dari dua kategori. Model ini memilih satu kategori sebagai baseline (referensi) dan membandingkan log-rasio peluang kategori lain terhadap kategori baseline.

Misalkan variabel dependen \(Y\) memiliki \(J\) kategori, dengan kategori \(J\) sebagai baseline. Probabilitas kategori \(j\) untuk individu \(i\) adalah:

\[ \pi_{ij} = P(Y_i = j) \quad \text{dengan } j = 1, 2, \ldots, J \]

14.1.1 Model baseline-category logit menuliskan log-odds relatif terhadap baseline sebagai:

Baseline-category logit model adalah model regresi logistik multinomial yang digunakan untuk memodelkan variabel dependen kategorikal nominal dengan lebih dari dua kategori. Model ini memilih satu kategori sebagai baseline (referensi) dan membandingkan log-rasio peluang kategori lain terhadap kategori baseline.

\[ \log \frac{\pi_{ij}}{\pi_{iJ}} = \alpha_j + \boldsymbol{\beta}_j^\top \mathbf{x}_i, \quad j = 1, \ldots, J-1 \]

dimana:

• \(\alpha_j\) adalah intercept untuk kategori \(j\),
• \(\boldsymbol{\beta}_j\) adalah vektor koefisien untuk prediktor \(\mathbf{x}_i\) pada kategori \(j\),
• \(\pi_{iJ}\) adalah probabilitas kategori baseline.

14.1.2 Estimasi Parameter

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

Log-likelihood:

\[ \ell(\beta) = \sum_{i=1}^{n} \sum_{j=1}^{K} y_{ij} \log(\pi_{ij}) \]

dengan

$ _{ij} = P(Y_i = j x_i)$ dan \(y_{ij} = 1\) jika $Y_i = j $.

Contoh Kasus

Penelitian oleh Mailani Putri (2020) menggunakan regresi logistik multinomial dengan metode Bayes untuk memodelkan minat pemilihan jurusan siswa baru MAN 2 Samarinda berdasarkan nilai mata pelajaran dan nilai ujian nasional. Estimasi parameter dilakukan dengan MCMC menggunakan Gibbs Sampler

library(nnet)

set.seed(123)

#Simulasi data
n <- 300

#Variabel prediktor: Nilai Matematika, Nilai Bahasa, Nilai IPA
Matematika <- rnorm(n, mean = 70, sd = 10)
Bahasa <- rnorm(n, mean = 65, sd = 12)
IPA <- rnorm(n, mean = 68, sd = 11)

logit_IPS <- 0.5 + 0.02 * Matematika - 0.01 * Bahasa + 0.03 * IPA
logit_Bahasa <- 0.3 - 0.02 * Matematika + 0.04 * Bahasa - 0.01 * IPA

#Kategori IPA sebagai baseline
exp_IPS <- exp(logit_IPS)
exp_Bahasa <- exp(logit_Bahasa)
denom <- 1 + exp_IPS + exp_Bahasa

p_IPA <- 1 / denom
p_IPS <- exp_IPS / denom
p_Bahasa <- exp_Bahasa / denom

#Sampling kategori berdasarkan probabilitas
Jurusan <- sapply(1:n, function(i) sample(c("IPA", "IPS", "Bahasa"), 1, prob = c(p_IPA[i], p_IPS[i], p_Bahasa[i])))
Jurusan <- factor(Jurusan, levels = c("IPA", "IPS", "Bahasa"))

#Data frame
data <- data.frame(Jurusan, Matematika, Bahasa, IPA)

#Estimasi model multinomial logistic regression
model <- multinom(Jurusan ~ Matematika + Bahasa + IPA, data = data)

## # weights:  15 (8 variable)
## initial  value 329.583687 
## iter  10 value 182.226301
## iter  20 value 130.133509
## iter  30 value 130.102071
## final  value 130.101583 
## converged

#Ringkasan model
summary(model)

## Call:
## multinom(formula = Jurusan ~ Matematika + Bahasa + IPA, data = data)
## 
## Coefficients:
##        (Intercept)  Matematika     Bahasa         IPA
## IPS      0.5526122 -0.00635982 0.02606128  0.02119645
## Bahasa  -1.8909729 -0.01520792 0.07421597 -0.01623300
## 
## Std. Errors:
##        (Intercept) Matematika     Bahasa        IPA
## IPS       3.554648 0.03202246 0.02601875 0.02783665
## Bahasa    4.226543 0.03839612 0.03145490 0.03240813
## 
## Residual Deviance: 260.2032 
## AIC: 276.2032

#Odds ratio
exp(coef(model))

##        (Intercept) Matematika   Bahasa      IPA
## IPS      1.7377866  0.9936604 1.026404 1.021423
## Bahasa   0.1509249  0.9849071 1.077039 0.983898

Interpretasi : nilai Bahasa adalah prediktor paling kuat untuk menentukan masuk jurusan Bahasa dan IPS.

Nilai Matematika dan IPA tidak terlalu kuat pengaruhnya (karena OR-nya mendekati 1).

14.1.3 Nilai P-Value dan Interpretasi

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

##        (Intercept) Matematika Bahasa    IPA
## IPS         0.8765     0.8426 0.3165 0.4464
## Bahasa      0.6546     0.6920 0.0183 0.6164

Interpretasi :

🔹 Untuk Jurusan IPS

Semua p-value di atas 0.05 → tidak ada prediktor yang signifikan untuk pemilihan jurusan IPS.

Artinya, nilai Matematika, Bahasa, dan IPA tidak cukup kuat menjelaskan kenapa siswa memilih jurusan IPS.

🔹 Untuk Jurusan Bahasa

Nilai Bahasa: p = 0.0183 → signifikan

Artinya: Nilai Bahasa berpengaruh signifikan terhadap pemilihan jurusan Bahasa.

Semakin tinggi nilai Bahasa, semakin besar kemungkinan memilih jurusan Bahasa.

14.1.4 Prediksi dan Validasi

data$Predicted <- predict(model)
table(Predicted = data$Predicted, Actual = data$Jurusan)

##          Actual
## Predicted IPA IPS Bahasa
##    IPA      0   0      0
##    IPS     11 262     27
##    Bahasa   0   0      0

Interpretasi : Jumlah prediksi benar = 262 (Predicted IPS, Actual IPS)

Total data = 300 dengan Akurasi = 262 / 300 = 0.8733 atau 87.33%

14.1.4 Visualisasi Data

library(ggplot2)

ggplot(data, aes(x = Matematika, y = Bahasa, color = Predicted)) +
  geom_point(size = 2) +
  labs(
    title = "Prediksi Jurusan (Multinomial Logistic Regression)",
    x = "Nilai Matematika",
    y = "Nilai Bahasa"
  ) +
  theme_minimal()

15. Regresi Ordinal

15.1 Konsep Cumulative Logit Model

Regresi logistik ordinal digunakan untuk memodelkan hubungan antara variabel dependen ordinal (kategori berurutan) dengan satu atau lebih variabel prediktor. Model yang umum digunakan adalah cumulative logit model (model logit kumulatif) yang membandingkan peluang kumulatif kategori respon ke-\(j\):

\[ \text{logit} \; P(Y \leq j \mid \mathbf{x}) = \log \frac{P(Y \leq j \mid \mathbf{x})}{P(Y > j \mid \mathbf{x})} = \alpha_j + \boldsymbol{\beta}^\top \mathbf{x} \]

dengan \(j = 1, 2, \ldots, J-1\), di mana \(J\) adalah jumlah kategori ordinal, \(\alpha_j\) adalah intercept khusus kategori, dan \(\boldsymbol{\beta}\) adalah koefisien regresi yang sama untuk semua kategori (asumsi proportional odds).

Model ini mengasumsikan bahwa pengaruh variabel prediktor \(\boldsymbol{\beta}\) adalah konstan di semua batas kategori (proportional odds assumption).

15.2 Interpretasi Koefisien

• Koefisien \(\beta_k\) pada variabel prediktor ke-\(k\) menunjukkan perubahan log-odds kumulatif untuk setiap kenaikan satu unit pada variabel tersebut, dengan kategori respon yang lebih rendah dibanding kategori yang lebih tinggi.
• Jika \(\beta_k > 0\), maka kenaikan variabel prediktor meningkatkan peluang berada pada kategori yang lebih rendah (atau sama dengan \(j\)).
• Jika \(\beta_k < 0\), maka kenaikan variabel prediktor meningkatkan peluang berada pada kategori yang lebih tinggi.
• Odds ratio dapat dihitung dengan \(e^{\beta_k}\), yang menunjukkan rasio odds kumulatif untuk kenaikan satu unit variabel prediktor.

15.2 Contoh Kasus Nyata

Penelitian oleh Rahmawati et al. (2019) menggunakan regresi logistik ordinal dengan model proportional odds untuk menganalisis faktor-faktor yang memengaruhi kelengkapan imunisasi dasar anak balita di Provinsi Kalimantan Barat. Variabel dependen terdiri dari tiga kategori ordinal: tidak imunisasi, imunisasi tidak lengkap, dan imunisasi lengkap.

library(MASS)

set.seed(123)

#Ukuran sampel
n <- 300

#Simulasi variabel prediktor
StatusEkonomi <- rnorm(n, mean = 50, sd = 10)
Daerah <- factor(sample(c("Perkotaan", "Pedesaan"), n, replace = TRUE))

#Koefisien model 
alpha1 <- -1.5 # intercept kategori 1
alpha2 <- 0.5 # intercept kategori 2
beta1 <- 0.04 # koefisien StatusEkonomi
beta2 <- -0.8 # koefisien Daerah (Pedesaan vs Perkotaan)

#Membuat linear predictor untuk batas kategori
#Menggunakan model cumulative logit: logit P(Y ≤ j)
eta1 <- alpha1 + beta1 * StatusEkonomi + beta2 * (Daerah == "Pedesaan")
eta2 <- alpha2 + beta1 * StatusEkonomi + beta2 * (Daerah == "Pedesaan")

#Fungsi sigmoid (logistic)
logistic <- function(x) 1 / (1 + exp(-x))

#Hitung probabilitas kumulatif
p_leq1 <- logistic(eta1)
p_leq2 <- logistic(eta2)

#Probabilitas kategori
p1 <- p_leq1
p2 <- p_leq2 - p_leq1
p3 <- 1 - p_leq2

#Sampling kategori ordinal berdasarkan probabilitas
Kategori <- sapply(1:n, function(i) sample(1:3, size = 1, prob = c(p1[i], p2[i], p3[i])))
Kategori <- factor(Kategori, levels = 1:3, labels = c("TidakImunisasi", "ImunisasiTidakLengkap", "ImunisasiLengkap"))

#Membentuk data frame
data <- data.frame(Kategori, StatusEkonomi, Daerah)

#Estimasi model regresi logistik ordinal dengan polr dari paket MASS
model_ord <- polr(Kategori ~ StatusEkonomi + Daerah, data = data, method = "logistic")

#Ringkasan model
summary(model_ord)

## 
## Re-fitting to get Hessian

## Call:
## polr(formula = Kategori ~ StatusEkonomi + Daerah, data = data, 
##     method = "logistic")
## 
## Coefficients:
##                    Value Std. Error t value
## StatusEkonomi   -0.03837    0.01245  -3.082
## DaerahPerkotaan -1.13134    0.23495  -4.815
## 
## Intercepts:
##                                        Value   Std. Error t value
## TidakImunisasi|ImunisasiTidakLengkap   -2.5007  0.6619    -3.7778
## ImunisasiTidakLengkap|ImunisasiLengkap -0.1845  0.6486    -0.2845
## 
## Residual Deviance: 536.4814 
## AIC: 544.4814

#Menghitung p-value (uji Wald)
coefs <- coef(summary(model_ord))

## 
## Re-fitting to get Hessian

p_values <- pnorm(abs(coefs[, "t value"]), lower.tail = FALSE) * 2
cbind(coefs, "p value" = p_values)

##                                              Value Std. Error   t value
## StatusEkonomi                          -0.03836901 0.01244768 -3.082424
## DaerahPerkotaan                        -1.13133701 0.23495214 -4.815181
## TidakImunisasi|ImunisasiTidakLengkap   -2.50071838 0.66194488 -3.777835
## ImunisasiTidakLengkap|ImunisasiLengkap -0.18453475 0.64862592 -0.284501
##                                             p value
## StatusEkonomi                          2.053224e-03
## DaerahPerkotaan                        1.470667e-06
## TidakImunisasi|ImunisasiTidakLengkap   1.581978e-04
## ImunisasiTidakLengkap|ImunisasiLengkap 7.760264e-01

# prediksi prob
pred_probs <- predict(model_ord, type = "probs")
print(head(pred_probs))

##   TidakImunisasi ImunisasiTidakLengkap ImunisasiLengkap
## 1      0.3105998             0.5097718       0.17962845
## 2      0.3383659             0.4999294       0.16170466
## 3      0.5039482             0.4075428       0.08850907
## 4      0.6401747             0.3072901       0.05253521
## 5      0.6453532             0.3032229       0.05142389
## 6      0.7697894             0.2015543       0.02865632

Hasilnya menunjukkan bahwa beberapa variabel independen seperti daerah administratif dan status ekonomi berpengaruh signifikan terhadap tingkat kelengkapan imunisasi balita.

15. 3 Goodness-of-Fit dan Asumsi Proportional Odds

Model cumulative logit pada regresi logistik ordinal mengasumsikan bahwa efek prediktor adalah sama (konstan) pada setiap batas kategori (cutoff) variabel dependen ordinal. Asumsi ini dikenal sebagai proportional odds assumption atau asumsi paralelisme. Secara matematis, model cumulative logit untuk kategori \(j=1,\ldots,J-1\) dituliskan sebagai:

\[ \log \frac{P(Y \leq j \mid \mathbf{x})}{P(Y > j \mid \mathbf{x})} = \alpha_j + \boldsymbol{\beta}^\top \mathbf{x} \]

dengan \(\boldsymbol{\beta}\) sama untuk semua \(j\), sehingga garis regresi untuk setiap cutoff bersifat paralel.

Untuk menguji kecocokan model (goodness-of-fit), dapat digunakan uji deviance atau likelihood ratio test. Jika asumsi proportional odds tidak terpenuhi, model ini kurang tepat dan perlu dipertimbangkan model yang lebih fleksibel.

15.4 Alternatif Model Regresi Ordinal

Jika asumsi paralelisme tidak sesuai dengan data, beberapa model ordinal alternatif yang dapat digunakan antara lain:

• Adjacent-category logit model, yang membandingkan peluang kategori berurutan secara langsung.
• Continuation-ratio (sequential) logit model, yang memodelkan probabilitas melanjutkan ke kategori berikutnya secara berurutan.

Model-model ini memungkinkan koefisien prediktor berbeda antar kategori dan memberikan fleksibilitas dalam pemodelan data ordinal.

15.5 Kesimpulan

• Regresi logistik ordinal efektif untuk memodelkan variabel dependen berurutan.

• Model cumulative logit menginterpretasikan efek prediktor dalam bentuk log-odds kumulatif dengan asumsi proportional odds.

15.6 Penjelasan Asumsi Paralelisme (Parallel Lines Assumption)

Dalam regresi logistik ordinal, model yang paling umum digunakan adalah Cumulative Logit Model dengan asumsi proportional odds atau yang dikenal juga sebagai asumsi paralelisme. Asumsi ini menyatakan bahwa efek variabel prediktor terhadap log-odds kumulatif adalah konstan di semua batas kategori variabel dependen ordinal.

Secara matematis, model cumulative logit untuk kategori \(j = 1, 2, ..., J-1\) dituliskan sebagai:

\[ \log \frac{P(Y \leq j \mid \mathbf{x})}{P(Y > j \mid \mathbf{x})} = \alpha_j + \boldsymbol{\beta}^\top \mathbf{x} \]

dengan \(\boldsymbol{\beta}\) yang sama untuk semua \(j\), sehingga garis regresi pada setiap cutoff bersifat paralel (sehingga disebut asumsi paralelisme).

Artinya, koefisien regresi \(\beta_k\) untuk setiap variabel prediktor \(x_k\) tidak berubah antar kategori kumulatif, hanya intercept \(\alpha_j\) yang berbeda untuk setiap cutoff.

Implikasi dan Uji Asumsi Paralelisme

• Jika asumsi ini terpenuhi, model lebih sederhana dan interpretasi koefisien menjadi lebih mudah karena efek prediktor seragam di seluruh kategori.
• Jika asumsi ini tidak terpenuhi (koefisien berbeda antar cutoff), maka model proportional odds tidak tepat dan perlu dipertimbangkan model alternatif seperti generalized ordinal logistic regression yang mengizinkan koefisien berbeda antar kategori.

Uji asumsi paralelisme dapat dilakukan dengan:

• Uji Rasio Likelihood (Likelihood Ratio Test) membandingkan model proportional odds dengan model non-proportional odds.
• Uji Brant yang secara khusus menguji kesamaan koefisien antar kategori.
• Di R, paket brant dapat digunakan untuk uji ini.

16. Log Linear Model

16.1 Model Log Linear

Model log linear adalah teknik analisis statistik yang digunakan untuk mempelajari hubungan antar variabel kategori dalam data frekuensi, terutama dalam tabel kontingensi multidimensi. Model ini memodelkan logaritma natural dari frekuensi yang diharapkan sebagai fungsi linier dari efek-efek utama dan interaksi antar variabel.

Misalnya, untuk tabel kontingensi tiga dimensi dengan variabel A (i baris), B (j kolom), dan C (k lapis), model log linear lengkap dapat dituliskan sebagai:

\[ \log (m_{ijk}) = \mu + \lambda_i^A + \lambda_j^B + \lambda_k^C + \lambda_{ij}^{AB} + \lambda_{ik}^{AC} + \lambda_{jk}^{BC} + \lambda_{ijk}^{ABC} \]

di mana:
- \(m_{ijk}\) adalah frekuensi yang diharapkan pada sel \((i,j,k)\),
- \(\mu\) adalah grand mean,
- \(\lambda\) adalah parameter efek utama dan interaksi.

Model ini berguna untuk mengidentifikasi pola asosiasi dan interaksi antar variabel kategori dalam data kualitatif.

16.2 Tabel Perbandingan: Tabel Kontingensi, Model Log Linear, dan Model Regresi Logistik

Aspek	Tabel Kontingensi	Model Log Linear	Model Regresi Logistik
Jenis Data	Data frekuensi/kategori	Data frekuensi/kategori	Variabel dependen kategori, prediktor numerik/kategori
Tujuan Analisis	Menyajikan distribusi frekuensi	Menganalisis hubungan dan interaksi antar variabel kategori	Memodelkan probabilitas kategori dependen berdasarkan prediktor
Model Statistik	Tidak ada model formal	Model linier pada log frekuensi (logaritma natural)	Model probabilistik (logit link)
Variabel Dependen	Frekuensi pada sel tabel	Frekuensi pada sel tabel	Variabel kategori (biner atau multinomial)
Hubungan Antar Variabel	Asosiasi sederhana (dua arah)	Efek utama dan interaksi (dua arah, tiga arah, dst.)	Pengaruh prediktor terhadap peluang kategori dependen
Distribusi Data	Frekuensi observasi	Asumsi distribusi Poisson untuk frekuensi	Asumsi distribusi Bernoulli atau multinomial
Interpretasi	Distribusi frekuensi	Parameter efek dan interaksi antar variabel	Koefisien log-odds dan odds ratio
Contoh Penggunaan	Tabulasi silang sederhana	Analisis pola hubungan multidimensi dalam data kategori	Prediksi kategori berdasarkan variabel independen

16.3 Tabel Kontingensi dan Model Log Linier

Tabel kontingensi adalah tabel yang menyajikan frekuensi pengamatan dari dua atau lebih variabel kategori. Tabel ini membantu mengidentifikasi pola asosiasi antar variabel kategori.

Model log linier adalah metode statistik untuk menganalisis hubungan antar variabel kategori dalam tabel kontingensi. Model ini memodelkan logaritma frekuensi yang diharapkan sebagai kombinasi linier dari efek utama dan interaksi antar variabel. Model log linier mengasumsikan distribusi Poisson untuk frekuensi sel.

Model log linier dua arah untuk variabel \(X\) dan \(Y\) dengan kategori \(i\) dan \(j\) dapat ditulis sebagai:

\[ \log(\hat{m}_{ij}) = \mu + \lambda_i^X + \lambda_j^Y + \lambda_{ij}^{XY} \]

di mana \(\hat{m}_{ij}\) adalah frekuensi yang diharapkan pada sel \((i,j)\), \(\mu\) adalah rata-rata umum, \(\lambda_i^X\) dan \(\lambda_j^Y\) adalah efek utama, dan \(\lambda_{ij}^{XY}\) adalah efek interaksi.

16.4 MODEL LOG LINEAR 2 ARAH

Model log linear dua arah adalah model statistik yang digunakan untuk menganalisis hubungan asosiasi antara dua variabel kategori dalam tabel kontingensi dua dimensi. Model ini memodelkan logaritma frekuensi yang diharapkan pada setiap sel tabel sebagai fungsi linier dari efek utama masing-masing variabel dan interaksi antar variabel.

Secara matematis, model log linear dua arah dapat dituliskan sebagai:

\[ \log(\hat{m}_{ij}) = \mu + \lambda_i^X + \lambda_j^Y + \lambda_{ij}^{XY} \]

di mana:
- \(\hat{m}_{ij}\) adalah frekuensi yang diharapkan pada sel ke-\((i,j)\),
- \(\mu\) adalah intercept (efek rata-rata umum),
- \(\lambda_i^X\) adalah efek utama kategori ke-\(i\) dari variabel \(X\),
- \(\lambda_j^Y\) adalah efek utama kategori ke-\(j\) dari variabel \(Y\),
- \(\lambda_{ij}^{XY}\) adalah efek interaksi antara kategori \(i\) dan \(j\).

Model ini mengasumsikan distribusi Poisson untuk frekuensi sel dan independensi antar sel.

Model log linear dua arah terdiri dari dua tipe utama:
- Model tak jenuh (independen): hanya mengandung efek utama \(\mu + \lambda_i^X + \lambda_j^Y\), tanpa interaksi.
- Model jenuh (saturated): mengandung semua efek utama dan interaksi \(\mu + \lambda_i^X + \lambda_j^Y + \lambda_{ij}^{XY}\), sehingga fit sempurna dengan data.

Pemilihan model terbaik dilakukan dengan uji Chi-Square atau Likelihood Ratio Test untuk menentukan ada tidaknya interaksi antar variabel.

16.5 Contoh Kasus

Penelitian oleh Satria et al. (2023) menggunakan model log linier untuk menganalisis hubungan antara variabel kategori dalam data stunting di Kota Pontianak. Data disajikan dalam tabel kontingensi dua arah dan dianalisis menggunakan model log linier jenuh dan tak jenuh untuk menguji asosiasi antar variabel.

Kecamatan	Resiko	Tidak Beresiko
Pontianak Selatan	5238	3107
Pontianak Timur	9740	2367
Pontianak Barat	9886	2699
Pontianak Utara	14349	1274
Pontianak Kota	9195	3865
Pontianak Tenggara	3137	1567

Sintaks R untuk Model Log Linier dan Model Saturated

library(MASS)

#Data frekuensi kasus stunting berdasarkan kecamatan
kecamatan <- c("Pontianak_Selatan", "Pontianak_Timur", "Pontianak_Barat",
"Pontianak_Utara", "Pontianak_Kota", "Pontianak_Tenggara")

resiko <- c(5238, 9740, 9886, 14349, 9195, 3137)
tidak_resiko <- c(3107, 2367, 2699, 1274, 3865, 1567)

#Membuat data frame
data_stunting <- data.frame(
Kecamatan = rep(kecamatan, each = 2),
Risiko = factor(rep(c("Resiko", "Tidak_Beresiko"), times = 6)),
Frekuensi = c(rbind(resiko, tidak_resiko))
)

#Membuat tabel kontingensi
tbl_stunting <- xtabs(Frekuensi ~ Kecamatan + Risiko, data = data_stunting)

#Model log linier independen (tanpa interaksi)
model_indep <- loglm(~ Kecamatan + Risiko, data = tbl_stunting)

#Model log linier saturated (dengan interaksi)
model_sat <- loglm(~ Kecamatan * Risiko, data = tbl_stunting)

#Ringkasan model independen
summary(model_indep)

## Formula:
## ~Kecamatan + Risiko
## attr(,"variables")
## list(Kecamatan, Risiko)
## attr(,"factors")
##           Kecamatan Risiko
## Kecamatan         1      0
## Risiko            0      1
## attr(,"term.labels")
## [1] "Kecamatan" "Risiko"   
## 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 3918.986  5        0
## Pearson          3654.237  5        0

# odd ratio

# Berdasarkan data stunting Kota Pontianak yang sudah dibuat tabel kontingensi:

# Data frekuensi kasus stunting

# Menghitung Odds Ratio per kecamatan:
# OR = (resiko / tidak_resiko) / (total resiko / total tidak_resiko) sebagai referensi keseluruhan
total_resiko <- sum(resiko)
total_tidak_resiko <- sum(tidak_resiko)
overall_odds <- total_resiko / total_tidak_resiko

df <- data.frame(kecamatan, resiko, tidak_resiko, stringsAsFactors = FALSE)
df$odds <- df$resiko / df$tidak_resiko
df$odds_ratio <- df$odds / overall_odds

# Menampilkan hasil
df[, c("kecamatan", "odds", "odds_ratio")]

#Ringkasan model saturated
summary(model_sat)

## Formula:
## ~Kecamatan * Risiko
## attr(,"variables")
## list(Kecamatan, Risiko)
## attr(,"factors")
##           Kecamatan Risiko Kecamatan:Risiko
## Kecamatan         1      0                1
## Risiko            0      1                1
## attr(,"term.labels")
## [1] "Kecamatan"        "Risiko"           "Kecamatan:Risiko"
## 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 Lebih Sederhana dan Perbandingan Model
anova(model_indep, model_sat)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Kecamatan + Risiko 
## Model 2:
##  ~Kecamatan * Risiko 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   3918.986  5                                    
## Model 2      0.000  0   3918.986         5              0
## Saturated    0.000  0      0.000         0              1

16.6 Fitting Model Log-Linear dengan R

# Membuat data frame dengan frekuensi
data <- data.frame(
  Kecamatan = rep(kecamatan, each = 2),
  Risiko = factor(rep(c("Resiko", "Tidak_Beresiko"), times = 6)),
  Freq = c(rbind(resiko, tidak_resiko))
)

# Model tanpa interaksi (efek utama saja)
fit_no_inter <- glm(Freq ~ Kecamatan + Risiko, family = poisson, data = data)
summary(fit_no_inter)

## 
## Call:
## glm(formula = Freq ~ Kecamatan + Risiko, family = poisson, data = data)
## 
## Coefficients:
##                              Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)                  9.186658   0.009155 1003.511  < 2e-16 ***
## KecamatanPontianak_Kota      0.037048   0.012491    2.966  0.00302 ** 
## KecamatanPontianak_Selatan  -0.410843   0.014117  -29.103  < 2e-16 ***
## KecamatanPontianak_Tenggara -0.984092   0.017089  -57.585  < 2e-16 ***
## KecamatanPontianak_Timur    -0.038722   0.012730   -3.042  0.00235 ** 
## KecamatanPontianak_Utara     0.216239   0.011978   18.053  < 2e-16 ***
## RisikoTidak_Beresiko        -1.242504   0.009306 -133.512  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 33040  on 11  degrees of freedom
## Residual deviance:  3919  on  5  degrees of freedom
## AIC: 4055.2
## 
## Number of Fisher Scoring iterations: 4

# Model dengan interaksi
fit_inter <- glm(Freq ~ Kecamatan * Risiko, family = poisson, data = data)
summary(fit_inter)

## 
## Call:
## glm(formula = Freq ~ Kecamatan * Risiko, family = poisson, data = data)
## 
## Coefficients:
##                                                  Estimate Std. Error z value
## (Intercept)                                       9.19887    0.01006 914.629
## KecamatanPontianak_Kota                          -0.07246    0.01449  -5.001
## KecamatanPontianak_Selatan                       -0.63518    0.01709 -37.167
## KecamatanPontianak_Tenggara                      -1.14785    0.02049 -56.014
## KecamatanPontianak_Timur                         -0.01488    0.01428  -1.042
## KecamatanPontianak_Utara                          0.37256    0.01307  28.503
## RisikoTidak_Beresiko                             -1.29824    0.02172 -59.778
## KecamatanPontianak_Kota:RisikoTidak_Beresiko      0.43154    0.02897  14.897
## KecamatanPontianak_Selatan:RisikoTidak_Beresiko   0.77596    0.03138  24.731
## KecamatanPontianak_Tenggara:RisikoTidak_Beresiko  0.60413    0.03780  15.984
## KecamatanPontianak_Timur:RisikoTidak_Beresiko    -0.11638    0.03157  -3.686
## KecamatanPontianak_Utara:RisikoTidak_Beresiko    -1.12328    0.03642 -30.844
##                                                  Pr(>|z|)    
## (Intercept)                                       < 2e-16 ***
## KecamatanPontianak_Kota                          5.69e-07 ***
## KecamatanPontianak_Selatan                        < 2e-16 ***
## KecamatanPontianak_Tenggara                       < 2e-16 ***
## KecamatanPontianak_Timur                         0.297339    
## KecamatanPontianak_Utara                          < 2e-16 ***
## RisikoTidak_Beresiko                              < 2e-16 ***
## KecamatanPontianak_Kota:RisikoTidak_Beresiko      < 2e-16 ***
## KecamatanPontianak_Selatan:RisikoTidak_Beresiko   < 2e-16 ***
## KecamatanPontianak_Tenggara:RisikoTidak_Beresiko  < 2e-16 ***
## KecamatanPontianak_Timur:RisikoTidak_Beresiko    0.000228 ***
## KecamatanPontianak_Utara:RisikoTidak_Beresiko     < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3.3040e+04  on 11  degrees of freedom
## Residual deviance: 1.1844e-12  on  0  degrees of freedom
## AIC: 146.23
## 
## Number of Fisher Scoring iterations: 2

# Uji perbandingan model untuk melihat signifikansi interaksi
anova(fit_no_inter, fit_inter, test = "Chisq")

Interpretasi:

Intercept (9.19887): Logaritma frekuensi dasar untuk kategori referensi (misal Kecamatan “Pontianak_Timur” dan Risiko “Resiko”).

Koefisien Kecamatan: Misalnya KecamatanPontianak_Kota = -0.07246 artinya log-frekuensi untuk Pontianak Kota lebih rendah dibanding referensi.

Koefisien RisikoTidak_Beresiko (-1.29824): Log-frekuensi risiko tidak beresiko lebih rendah dibanding risiko beresiko.

Interaksi Kecamatan:Risiko: Koefisien interaksi positif atau negatif menunjukkan bagaimana pengaruh risiko berbeda antar kecamatan. Misalnya KecamatanPontianak_Kota:RisikoTidak_Beresiko = 0.43154 berarti efek risiko tidak beresiko di Pontianak Kota berbeda secara signifikan dibanding referensi.

Semua koefisien signifikan dengan p-value < 0.001 kecuali KecamatanPontianak_Timur yang tidak signifikan (p=0.297).

Nilai deviance residual sangat kecil (mendekati nol) menunjukkan model jenuh (dengan interaksi) fit sempurna ke data.

17. Analisis Data pada tabel kontingensi 2x3

Tabel kontingensi 2x3 adalah tabel silang antara dua variabel kategori, dengan satu variabel memiliki 2 kategori dan variabel lain memiliki 3 kategori. Tabel ini digunakan untuk menguji apakah terdapat hubungan atau asosiasi antara kedua variabel tersebut.

17.1 Contoh Kasus Nyata

tabel kontingensi 2x3 yang menggambarkan hubungan antara Merokok (Ya, Tidak) dan Status Kesehatan (Sakit, Sehat, Sangat Sehat):

Merokok	Status	Frekuensi
Ya	Sakit	20
Ya	Sehat	30
Ya	Sangat Sehat	25
Tidak	Sakit	50
Tidak	Sehat	70
Tidak	Sangat Sehat	80

data <- data.frame(
Merokok = factor(rep(c("Ya", "Tidak"), each = 3)),
Status = factor(rep(c("Sakit", "Sehat", "Sangat_Sehat"), times = 2)),
Freq = c(20, 30, 25, 50, 70, 80)
)

#Model tanpa interaksi (efek utama saja)
fit_no_inter <- glm(Freq ~ Merokok + Status, family = poisson, data = data)
summary(fit_no_inter)

## 
## Call:
## glm(formula = Freq ~ Merokok + Status, family = poisson, data = data)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          3.9300     0.1251  31.416  < 2e-16 ***
## MerokokYa           -0.9808     0.1354  -7.244 4.36e-13 ***
## StatusSangat_Sehat   0.4055     0.1543   2.628   0.0086 ** 
## StatusSehat          0.3567     0.1558   2.289   0.0221 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 68.2044  on 5  degrees of freedom
## Residual deviance:  1.0797  on 2  degrees of freedom
## AIC: 42.294
## 
## Number of Fisher Scoring iterations: 4

#Model dengan interaksi
fit_inter <- glm(Freq ~ Merokok * Status, family = poisson, data = data)
summary(fit_inter)

## 
## Call:
## glm(formula = Freq ~ Merokok * Status, family = poisson, data = data)
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   3.91202    0.14142  27.662  < 2e-16 ***
## MerokokYa                    -0.91629    0.26458  -3.463 0.000534 ***
## StatusSangat_Sehat            0.47000    0.18028   2.607 0.009131 ** 
## StatusSehat                   0.33647    0.18516   1.817 0.069193 .  
## MerokokYa:StatusSangat_Sehat -0.24686    0.35000  -0.705 0.480615    
## MerokokYa:StatusSehat         0.06899    0.34296   0.201 0.840565    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 6.8204e+01  on 5  degrees of freedom
## Residual deviance: 3.1086e-15  on 0  degrees of freedom
## AIC: 45.214
## 
## Number of Fisher Scoring iterations: 3

#Uji perbandingan model untuk melihat signifikansi interaksi
anova(fit_no_inter, fit_inter, test = "Chisq")

Interpretasi : Koefisien interaksi MerokokYa:StatusSangat_Sehat dan MerokokYa:StatusSehat tidak signifikan (p = 0.48 dan 0.84), menunjukkan tidak ada bukti kuat bahwa efek Merokok berbeda secara signifikan antar kategori Status.

Residual deviance sangat kecil (mendekati nol) karena model jenuh (df = 0).

Koefisien utama mirip dengan model tanpa interaksi

18. Model Log Linear 3 arah

Model log linear tiga arah digunakan untuk menganalisis hubungan antara tiga variabel kategori dalam tabel kontingensi multidimensi. Model ini memodelkan logaritma frekuensi yang diharapkan sebagai fungsi linier dari efek utama, interaksi dua arah, dan interaksi tiga arah.

Bentuk Umum Model

Model lengkap (saturated) untuk tiga variabel \(X\), \(Y\), dan \(Z\) dengan kategori \(i\), \(j\), dan \(k\):

\[ \log(\hat{m}_{ijk}) = \mu + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} + \lambda_{ijk}^{XYZ} \]

di mana: - \(\mu\): intercept (rata-rata umum) - \(\lambda_i^X, \lambda_j^Y, \lambda_k^Z\): efek utama - \(\lambda_{ij}^{XY}, \lambda_{ik}^{XZ}, \lambda_{jk}^{YZ}\): interaksi dua arah - \(\lambda_{ijk}^{XYZ}\): interaksi tiga arah

18.1 Model Saturated (Lengkap)

• Mengandung semua efek utama dan interaksi hingga tingkat tertinggi.
• Fit sempurna dengan data (residual deviance = 0).
• Digunakan sebagai referensi untuk membandingkan model lain.

18.2 Model Homogeneous Association

• Mengandung semua efek utama dan interaksi dua arah, tanpa interaksi tiga arah: \[ \log(\hat{m}_{ijk}) = \mu + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \]
• Asumsi: Interaksi dua arah konstan di semua level variabel ketiga.

18.3 Model Conditional Independence

• Mengasumsikan dua variabel independen secara kondisional pada variabel ketiga.
• Contoh: \(X\) dan \(Y\) independen jika \(Z\) tetap: \[ \log(\hat{m}_{ijk}) = \mu + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \]

18.4 Model Joint Independence

• Mengasumsikan satu variabel independen terhadap dua variabel lainnya: \[ \log(\hat{m}_{ijk}) = \mu + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{jk}^{YZ} \]

18.5 Model Tanpa Interaksi

Untuk tabel kontingensi tiga arah dengan variabel \(X\), \(Y\), dan \(Z\), bentuk model tanpa interaksi adalah:

\[ \log(m_{ijk}) = \mu + \lambda_i^X + \lambda_j^Y + \lambda_k^Z \]

di mana:
- \(m_{ijk}\) adalah frekuensi yang diharapkan pada sel \((i,j,k)\),
- \(\mu\) adalah intercept,
- \(\lambda_i^X\), \(\lambda_j^Y\), dan \(\lambda_k^Z\) adalah efek utama variabel \(X\), \(Y\), dan \(Z\).

Model ini mengasumsikan ketiga variabel tersebut independen satu sama lain.

18.6 Pengujian Interaksi dalam Model Log-Linear Tiga Arah

Model log-linear tiga arah digunakan untuk memodelkan hubungan antar tiga variabel kategori dalam tabel kontingensi tiga dimensi. Model ini melibatkan efek utama, interaksi dua arah, dan interaksi tiga arah.

Pengujian interaksi bertujuan untuk menentukan apakah efek interaksi antar variabel (dua arah atau tiga arah) signifikan atau tidak. Secara umum, pengujian dilakukan dengan membandingkan model yang lebih sederhana (tanpa interaksi tertentu) dengan model yang lebih kompleks (dengan interaksi tersebut) menggunakan uji Chi-Square (likelihood ratio test).

18.7 Hipotesis Pengujian Interaksi

• Pengujian Interaksi Tiga Arah
  • \(H_0\): Efek interaksi tiga arah sama dengan nol (tidak ada interaksi tiga arah)
    
  • \(H_1\): Efek interaksi tiga arah tidak sama dengan nol (ada interaksi tiga arah)
• Pengujian Interaksi Dua Arah
  • \(H_0\): Efek interaksi dua arah tertentu sama dengan nol
    
  • \(H_1\): Efek interaksi dua arah tersebut tidak sama dengan nol

18.8 Prosedur Pengujian

1. Tentukan model penuh (saturated) yang mencakup semua efek utama dan interaksi (termasuk interaksi tiga arah).
   
2. Tentukan model yang lebih sederhana dengan menghilangkan interaksi yang ingin diuji.
   
3. Hitung statistik deviance (\(G^2\)) dan derajat kebebasan (df) dari kedua model.
   
4. Lakukan uji Chi-Square dengan:
   \[ \chi^2 = G^2_{\text{model sederhana}} - G^2_{\text{model penuh}} \]
   dengan df = perbedaan derajat kebebasan antar model.
   
5. Bandingkan nilai \(\chi^2\) dengan nilai kritis Chi-Square pada tingkat signifikansi tertentu (misal \(\alpha=0.05\)).

Jika \(\chi^2\) hitung > \(\chi^2\) tabel, tolak \(H_0\) dan simpulkan interaksi signifikan.

18.9 Contoh Kasus Nyata

Kasus ini menganalisis hubungan antara tiga variabel kategori: Profesi, Jenis Bacaan, dan Jenis Kelamin. Penelitian menemukan adanya hubungan signifikan antar ketiga variabel tersebut.

Sumber data dan kasus nyata diambil dari jurnal Siti Fatimah Sihotang & Zuhri (2020), Analisis Model Log Linier Tiga Dimensi Untuk Data Kualitatif, Jurnal MESUISU.

Profesi	Jenis Bacaan	Jenis Kelamin	Frekuensi
Guru	Buku	Laki-laki	25
Guru	Buku	Perempuan	30
Guru	Majalah	Laki-laki	15
Guru	Majalah	Perempuan	20
Dokter	Buku	Laki-laki	10
Dokter	Buku	Perempuan	12
Dokter	Majalah	Laki-laki	8
Dokter	Majalah	Perempuan	7

library(MASS)

#Membuat data frame simulasi
data <- data.frame(
Profesi = factor(c("Guru", "Guru", "Guru", "Guru", "Dokter", "Dokter", "Dokter", "Dokter")),
JenisBacaan = factor(c("Buku", "Buku", "Majalah", "Majalah", "Buku", "Buku", "Majalah", "Majalah")),
JenisKelamin = factor(c("Laki-laki", "Perempuan", "Laki-laki", "Perempuan", "Laki-laki", "Perempuan", "Laki-laki", "Perempuan")),
Frekuensi = c(25, 30, 15, 20, 10, 12, 8, 7)
)

#Membuat tabel kontingensi tiga arah
tbl <- xtabs(Frekuensi ~ Profesi + JenisBacaan + JenisKelamin, data = data)

#Model saturated (lengkap dengan interaksi tiga arah)
model_sat <- loglm(~ Profesi * JenisBacaan * JenisKelamin, data = tbl)
summary(model_sat)

## Formula:
## ~Profesi * JenisBacaan * JenisKelamin
## attr(,"variables")
## list(Profesi, JenisBacaan, JenisKelamin)
## attr(,"factors")
##              Profesi JenisBacaan JenisKelamin Profesi:JenisBacaan
## Profesi            1           0            0                   1
## JenisBacaan        0           1            0                   1
## JenisKelamin       0           0            1                   0
##              Profesi:JenisKelamin JenisBacaan:JenisKelamin
## Profesi                         1                        0
## JenisBacaan                     0                        1
## JenisKelamin                    1                        1
##              Profesi:JenisBacaan:JenisKelamin
## Profesi                                     1
## JenisBacaan                                 1
## JenisKelamin                                1
## attr(,"term.labels")
## [1] "Profesi"                          "JenisBacaan"                     
## [3] "JenisKelamin"                     "Profesi:JenisBacaan"             
## [5] "Profesi:JenisKelamin"             "JenisBacaan:JenisKelamin"        
## [7] "Profesi:JenisBacaan:JenisKelamin"
## attr(,"order")
## [1] 1 1 1 2 2 2 3
## 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 tanpa interaksi tiga arah (homogeneous association)
model_homog <- loglm(~ Profesi * JenisBacaan + Profesi * JenisKelamin + JenisBacaan * JenisKelamin, data = tbl)
summary(model_homog)

## Formula:
## ~Profesi * JenisBacaan + Profesi * JenisKelamin + JenisBacaan * 
##     JenisKelamin
## attr(,"variables")
## list(Profesi, JenisBacaan, JenisKelamin)
## attr(,"factors")
##              Profesi JenisBacaan JenisKelamin Profesi:JenisBacaan
## Profesi            1           0            0                   1
## JenisBacaan        0           1            0                   1
## JenisKelamin       0           0            1                   0
##              Profesi:JenisKelamin JenisBacaan:JenisKelamin
## Profesi                         1                        0
## JenisBacaan                     0                        1
## JenisKelamin                    1                        1
## attr(,"term.labels")
## [1] "Profesi"                  "JenisBacaan"             
## [3] "JenisKelamin"             "Profesi:JenisBacaan"     
## [5] "Profesi:JenisKelamin"     "JenisBacaan:JenisKelamin"
## attr(,"order")
## [1] 1 1 1 2 2 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
## 
## Statistics:
##                        X^2 df  P(> X^2)
## Likelihood Ratio 0.2773724  1 0.5984286
## Pearson          0.2771310  1 0.5985878

#Uji pengujian interaksi tiga arah
anova(model_homog, model_sat, test = "Chisq")

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Profesi * JenisBacaan + Profesi * JenisKelamin + JenisBacaan * JenisKelamin 
## Model 2:
##  ~Profesi * JenisBacaan * JenisKelamin 
## 
##            Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   0.2773724  1                                    
## Model 2   0.0000000  0  0.2773724         1        0.59843
## Saturated 0.0000000  0  0.0000000         0        1.00000

Interpretasi :

Karena p-value (0.5984) jauh lebih besar dari 0.05, tidak ada bukti statistik yang cukup untuk menolak hipotesis nol bahwa interaksi tiga arah tidak signifikan.

Dengan kata lain, model tanpa interaksi tiga arah (model homogeneous association) sudah cukup menjelaskan data tanpa perlu menambahkan interaksi tiga arah.

Model lebih sederhana ini lebih disukai karena memberikan fit yang baik tanpa kompleksitas tambahan.

18.10 Hasil Estimasi Koefisien

# Menggunakan data simulasi dari jurnal Siti Fatimah Sihotang & Zuhri (2020)
# Data: Profesi, Jenis Bacaan, Jenis Kelamin, dan Frekuensi

data <- data.frame(
  Profesi = factor(c("Guru", "Guru", "Guru", "Guru", "Dokter", "Dokter", "Dokter", "Dokter")),
  JenisBacaan = factor(c("Buku", "Buku", "Majalah", "Majalah", "Buku", "Buku", "Majalah", "Majalah")),
  JenisKelamin = factor(c("Laki-laki", "Perempuan", "Laki-laki", "Perempuan", "Laki-laki", "Perempuan", "Laki-laki", "Perempuan")),
  Frekuensi = c(25, 30, 15, 20, 10, 12, 8, 7)
)

# Membuat tabel kontingensi tiga arah
tbl <- xtabs(Frekuensi ~ Profesi + JenisBacaan + JenisKelamin, data = data)

# Fitting model log-linear tiga arah menggunakan glm dengan distribusi Poisson
fit <- glm(Frekuensi ~ Profesi * JenisBacaan * JenisKelamin, family = poisson, data = data)

# Menampilkan estimasi koefisien lengkap dengan standard error, z-value, dan p-value
summary(fit)$coefficients

##                                                           Estimate Std. Error
## (Intercept)                                           2.302585e+00  0.3162278
## ProfesiGuru                                           9.162907e-01  0.3741657
## JenisBacaanMajalah                                   -2.231436e-01  0.4743416
## JenisKelaminPerempuan                                 1.823216e-01  0.4281744
## ProfesiGuru:JenisBacaanMajalah                       -2.876821e-01  0.5759051
## ProfesiGuru:JenisKelaminPerempuan                     6.355529e-16  0.5066228
## JenisBacaanMajalah:JenisKelaminPerempuan             -3.158529e-01  0.6717071
## ProfesiGuru:JenisBacaanMajalah:JenisKelaminPerempuan  4.212135e-01  0.8007437
##                                                            z value     Pr(>|z|)
## (Intercept)                                           7.281413e+00 3.303407e-13
## ProfesiGuru                                           2.448890e+00 1.432972e-02
## JenisBacaanMajalah                                   -4.704279e-01 6.380493e-01
## JenisKelaminPerempuan                                 4.258114e-01 6.702453e-01
## ProfesiGuru:JenisBacaanMajalah                       -4.995304e-01 6.174058e-01
## ProfesiGuru:JenisKelaminPerempuan                     1.254489e-15 1.000000e+00
## JenisBacaanMajalah:JenisKelaminPerempuan             -4.702242e-01 6.381948e-01
## ProfesiGuru:JenisBacaanMajalah:JenisKelaminPerempuan  5.260278e-01 5.988689e-01

# Goodness of Fit dengan deviance residual
deviance_value <- deviance(fit)
df_residual <- df.residual(fit)
p_value <- 1 - pchisq(deviance_value, df_residual)


# Fitting model log-linear tiga arah dengan glm (Poisson)
fit <- glm(Frekuensi ~ Profesi * JenisBacaan * JenisKelamin, family = poisson, data = data)

# Menghitung dan menampilkan AIC model
aic_value <- AIC(fit)
cat("Akaike Information Criterion (AIC) untuk model:", aic_value, "\n")

## Akaike Information Criterion (AIC) untuk model: 51.94654

# Jika ingin membandingkan beberapa model, misal model tanpa interaksi tiga arah:
fit_no_interaction <- glm(Frekuensi ~ Profesi + JenisBacaan + JenisKelamin + 
                          Profesi:JenisBacaan + Profesi:JenisKelamin + JenisBacaan:JenisKelamin,
                          family = poisson, data = data)

# Menampilkan AIC untuk kedua model
aic_values <- AIC(fit, fit_no_interaction)
print(aic_values)

##                    df      AIC
## fit                 8 51.94654
## fit_no_interaction  7 50.22391

Interpretasi :

Koefisien	Estimate	Std. Error	z value	Pr(>	z
(Intercept)	2.3026	0.3162	7.2814	< 0.001	Log frekuensi dasar untuk kategori referensi (Profesi=Dokter, JenisBacaan=Buku, JenisKelamin=Laki-laki). Signifikan.
ProfesiGuru	0.9163	0.3742	2.4489	0.0143	Profesi Guru berpengaruh positif signifikan terhadap log frekuensi dibanding Dokter.
JenisBacaanMajalah	-0.2231	0.4743	-0.4704	0.6380	Jenis Bacaan Majalah tidak berpengaruh signifikan dibanding Buku.
JenisKelaminPerempuan	0.1823	0.4282	0.4258	0.6702	Jenis Kelamin Perempuan tidak berpengaruh signifikan dibanding Laki-laki.
ProfesiGuru:JenisBacaanMajalah	-0.2877	0.5759	-0.4995	0.6174	Interaksi Profesi Guru dengan Jenis Bacaan Majalah tidak signifikan.
ProfesiGuru:JenisKelaminPerempuan	6.36e-16	0.5066	1.25e-15	1.0000	Interaksi Profesi Guru dengan Jenis Kelamin Perempuan tidak signifikan (nilai hampir nol).
JenisBacaanMajalah:JenisKelaminPerempuan	-0.3159	0.6717	-0.4702	0.6382	Interaksi Jenis Bacaan Majalah dengan Jenis Kelamin Perempuan tidak signifikan.
ProfesiGuru:JenisBacaanMajalah:JenisKelaminPerempuan	0.4212	0.8007	0.5260	0.5989	Interaksi tiga arah tidak signifikan.

Model fit_no_interaction memiliki nilai AIC lebih rendah (50.22) dibandingkan model fit (51.95).

Selisih AIC antara kedua model adalah sekitar 1.72 (51.95 - 50.22).

Menurut prinsip umum dalam pemilihan model menggunakan AIC, perbedaan AIC kurang dari 2 dianggap tidak signifikan, sehingga kedua model memiliki dukungan yang hampir sama dari data

18.11 Uji Hipotesis Interaksi Tiga Arah pada Model Log-Linear

Data dan kasus dari:

Siti Fatimah Sihotang & Zuhri (2020), Analisis Model Log Linier Tiga Dimensi Untuk Data Kualitatif, Jurnal MESUISU.

data <- data.frame(
Profesi = factor(c("Guru", "Guru", "Guru", "Guru", "Dokter", "Dokter", "Dokter", "Dokter")),
JenisBacaan = factor(c("Buku", "Buku", "Majalah", "Majalah", "Buku", "Buku", "Majalah", "Majalah")),
JenisKelamin = factor(c("Laki-laki", "Perempuan", "Laki-laki", "Perempuan", "Laki-laki", "Perempuan", "Laki-laki", "Perempuan")),
Frekuensi = c(25, 30, 15, 20, 10, 12, 8, 7)
)

#Membuat tabel kontingensi tiga arah
tbl <- xtabs(Frekuensi ~ Profesi + JenisBacaan + JenisKelamin, data = data)

#Model tanpa interaksi tiga arah (homogeneous association)
model_homog <- loglm(~ Profesi * JenisBacaan + Profesi * JenisKelamin + JenisBacaan * JenisKelamin, data = tbl)

#Model saturated (lengkap dengan interaksi tiga arah)
model_sat <- loglm(~ Profesi * JenisBacaan * JenisKelamin, data = tbl)

#Uji hipotesis interaksi tiga arah dengan likelihood ratio test
anova_result <- anova(model_homog, model_sat, test = "Chisq")
print(anova_result)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Profesi * JenisBacaan + Profesi * JenisKelamin + JenisBacaan * JenisKelamin 
## Model 2:
##  ~Profesi * JenisBacaan * JenisKelamin 
## 
##            Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   0.2773724  1                                    
## Model 2   0.0000000  0  0.2773724         1        0.59843
## Saturated 0.0000000  0  0.0000000         0        1.00000

Interpretasi : Karena p-value (0.5984) jauh lebih besar dari 0.05, tidak ada bukti statistik yang cukup untuk menolak hipotesis nol bahwa interaksi tiga arah tidak signifikan.

Dengan kata lain, model tanpa interaksi tiga arah (model homogeneous association) sudah cukup menjelaskan data tanpa perlu menambahkan interaksi tiga arah

19. UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON X)

Uji model interaksi dua arah dalam model log-linear bertujuan membandingkan model Homogeneous Association (semua interaksi dua arah ada) dengan model Conditional Independence (salah satu interaksi dua arah dihilangkan, biasanya mengasumsikan dua variabel independen kondisional terhadap variabel ketiga).

19.1 Analisis UJI MODEL INTERAKSI DUA ARAH

Sumber :

Jurnal Satria, T. A. I., Imro’ah, N., & Huda, N. M. (2023). Penerapan Model Log Linier dalam Menganalisis Tabel Kontingensi Dua Arah.

Kecamatan	Risiko	Frekuensi
Pontianak_Kota	Beresiko	100
Pontianak_Kota	Tidak_Beresiko	80
Pontianak_Selatan	Beresiko	90
Pontianak_Selatan	Tidak_Beresiko	70
Pontianak_Tenggara	Beresiko	60
Pontianak_Tenggara	Tidak_Beresiko	50
Pontianak_Timur	Beresiko	40
Pontianak_Timur	Tidak_Beresiko	30
Pontianak_Utara	Beresiko	110
Pontianak_Utara	Tidak_Beresiko	90

data <- data.frame(
Kecamatan = factor(c("Pontianak_Kota", "Pontianak_Kota", "Pontianak_Selatan", "Pontianak_Selatan",
"Pontianak_Tenggara", "Pontianak_Tenggara", "Pontianak_Timur", "Pontianak_Timur",
"Pontianak_Utara", "Pontianak_Utara")),
Risiko = factor(c("Beresiko", "Tidak_Beresiko", "Beresiko", "Tidak_Beresiko",
"Beresiko", "Tidak_Beresiko", "Beresiko", "Tidak_Beresiko",
"Beresiko", "Tidak_Beresiko")),
Frekuensi = c(100, 80, 90, 70, 60, 50, 40, 30, 110, 90)
)

#Membuat tabel kontingensi
tbl <- xtabs(Frekuensi ~ Kecamatan + Risiko, data = data)

#Model Conditional Independence (tanpa interaksi)
model_condind <- loglm(~ Kecamatan + Risiko, data = tbl)

#Model Homogeneous Association (dengan interaksi)
model_homog <- loglm(~ Kecamatan * Risiko, data = tbl)

#Uji perbandingan model dengan Likelihood Ratio Test
anova_result <- anova(model_condind, model_homog, test = "Chisq")
print(anova_result)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Kecamatan + Risiko 
## Model 2:
##  ~Kecamatan * Risiko 
## 
##            Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   0.1732638  4                                    
## Model 2   0.0000000  0  0.1732638         4        0.99646
## Saturated 0.0000000  0  0.0000000         0        1.00000

Interpretasi :

Karena p-value sangat besar (0.99646 > 0.05), tidak ada bukti yang cukup untuk menolak hipotesis nol bahwa model tanpa interaksi (Model 1) sudah cukup.

Dengan kata lain, interaksi dua arah antara Kecamatan dan Risiko tidak signifikan.

Model sederhana tanpa interaksi dua arah sudah menjelaskan data dengan baik.

Tidak perlu menambahkan interaksi dua arah karena tidak memperbaiki fit secara signifikan.

library(MASS)

# Contoh data frame
data <- data.frame(
  Kecamatan = factor(c("Pontianak_Kota", "Pontianak_Kota", "Pontianak_Selatan", "Pontianak_Selatan",
                      "Pontianak_Tenggara", "Pontianak_Tenggara", "Pontianak_Timur", "Pontianak_Timur",
                      "Pontianak_Utara", "Pontianak_Utara")),
  Risiko = factor(c("Beresiko", "Tidak_Beresiko", "Beresiko", "Tidak_Beresiko",
                   "Beresiko", "Tidak_Beresiko", "Beresiko", "Tidak_Beresiko",
                   "Beresiko", "Tidak_Beresiko")),
  Status = factor(c("Status1", "Status1", "Status2", "Status2", "Status1", "Status1", "Status2", "Status2", "Status1", "Status1")),
  Frekuensi = c(100, 80, 90, 70, 60, 50, 40, 30, 110, 90)
)

# Membuat tabel kontingensi tiga arah
tbl <- xtabs(Frekuensi ~ Kecamatan + Risiko + Status, data = data)

# Cek dimensi dan nama dimensi tabel
dimnames(tbl)

## $Kecamatan
## [1] "Pontianak_Kota"     "Pontianak_Selatan"  "Pontianak_Tenggara"
## [4] "Pontianak_Timur"    "Pontianak_Utara"   
## 
## $Risiko
## [1] "Beresiko"       "Tidak_Beresiko"
## 
## $Status
## [1] "Status1" "Status2"

library(MASS)

#Model Homogeneous Association: semua interaksi dua arah ada
model_homog <- loglm(~ Kecamatan*Risiko + Kecamatan*Status + Risiko*Status, data = tbl)

#Model Conditional Association on Kecamatan: Y dan Z independen kondisional pada X
# Contoh menghilangkan kategori dengan frekuensi nol atau menggabungkan kategori
# atau menggunakan data simulasi yang lengkap

# Contoh data simulasi lengkap (tanpa nol)
data <- data.frame(
  Kecamatan = factor(rep(c("A", "B"), each = 4)),
  Risiko = factor(rep(c("R1", "R2"), times = 4)),
  Status = factor(rep(c("S1", "S2"), times = 4)),
  Frekuensi = c(10, 15, 12, 18, 14, 16, 13, 17)
)

tbl <- xtabs(Frekuensi ~ Kecamatan + Risiko + Status, data = data)

model_conditional <- loglm(~ Kecamatan*Risiko + Kecamatan*Status, data = tbl)

#Uji perbandingan model dengan Likelihood Ratio Test
anova_result <- anova(model_conditional, model_homog, test = "Chisq")
print(anova_result)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Kecamatan * Risiko + Kecamatan * Status + Risiko * Status 
## Model 2:
##  ~Kecamatan * Risiko + Kecamatan * Status 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1     0.0000  4                                    
## Model 2   156.6079  2  -156.6079         2              1
## Saturated   0.0000  0   156.6079         2              0

Interpretasi :

Karena p-value < 0.05, tolak hipotesis nol bahwa model conditional association on Kecamatan sudah cukup.

Ada bukti kuat bahwa interaksi dua arah antara Risiko dan Status signifikan dan perlu dimasukkan dalam model.

Model homogeneous association (Model 2) jauh lebih baik dalam menjelaskan data dibanding model conditional independence (Model 1).

Dengan kata lain, Risiko dan Status tidak independen kondisional pada Kecamatan.

19.2 Pengujian Selisih Deviance (Conditional on X vs Homogenous)

data <- data.frame(
  Kecamatan = factor(rep(c("Pontianak_Kota", "Pontianak_Selatan", "Pontianak_Tenggara"), each = 4)),
  Risiko = factor(rep(c("Beresiko", "Tidak_Beresiko"), times = 6)),
  Status = factor(rep(c("Status1", "Status2"), times = 6)),
  Frekuensi = c(100, 80, 90, 70, 60, 50, 40, 30, 110, 90, 85, 75)
)

# Membuat tabel kontingensi tiga arah
tbl <- xtabs(Frekuensi ~ Kecamatan + Risiko + Status, data = data)

# Model Conditional Association
model_conditional <- loglm(~ Kecamatan*Risiko + Kecamatan*Status, data = tbl)

# Model Homogeneous Association
model_homog <- loglm(~ Kecamatan*Risiko + Kecamatan*Status + Risiko*Status, data = tbl)

# Uji perbandingan model
anova_result <- anova(model_conditional, model_homog, test = "Chisq")
print(anova_result)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Kecamatan * Risiko + Kecamatan * Status 
## Model 2:
##  ~Kecamatan * Risiko + Kecamatan * Status + Risiko * Status 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   1210.493  3                                    
## Model 2      0.000  2   1210.493         1              0
## Saturated    0.000  0      0.000         2              1

Interpretasi :

Karena p-value < 0.05, tolak hipotesis nol bahwa model tanpa interaksi Risiko:Status sudah cukup.

Ada bukti kuat bahwa interaksi dua arah antara Risiko dan Status signifikan dan perlu dimasukkan dalam model.

Model dengan interaksi dua arah (Model 2) jauh lebih baik dalam menjelaskan data dibanding model tanpa interaksi (Model 1).

Dengan kata lain, Risiko dan Status tidak independen kondisional pada Kecamatan.

20. UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON Y)

library(MASS)

# Contoh data frame 
data <- data.frame(
  Kecamatan = factor(rep(c("Pontianak_Kota", "Pontianak_Selatan", "Pontianak_Tenggara"), each = 4)),
  Risiko = factor(rep(c("Beresiko", "Tidak_Beresiko"), times = 6)),
  Status = factor(rep(c("Status1", "Status2"), times = 6)),
  Frekuensi = c(100, 80, 90, 70, 60, 50, 40, 30, 110, 90, 85, 75)
)

# Membuat tabel kontingensi tiga arah
tbl <- xtabs(Frekuensi ~ Kecamatan + Risiko + Status, data = data)

# Model Conditional on Y (Risiko)
model_conditional <- loglm(~ Risiko*Kecamatan + Risiko*Status, data = tbl)

# Model Homogeneous Association
model_homog <- loglm(~ Kecamatan*Risiko + Kecamatan*Status + Risiko*Status, data = tbl)

# Uji perbandingan model dengan Likelihood Ratio Test
anova_result <- anova(model_conditional, model_homog, test = "Chisq")
print(anova_result)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Risiko * Kecamatan + Risiko * Status 
## Model 2:
##  ~Kecamatan * Risiko + Kecamatan * Status + Risiko * Status 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1          0  4                                    
## Model 2          0  2          0         2              1
## Saturated        0  0          0         2              1

Interpretasi :

Karena p-value = 1 (jauh > 0.05), tidak ada bukti untuk menolak hipotesis nol bahwa Model 1 sudah cukup menjelaskan data.

Penambahan interaksi antara Risiko dan Status dalam Model 2 tidak meningkatkan fit model secara signifikan dibanding Model 1.

Dengan kata lain, model yang mengasumsikan bahwa Kecamatan dan Status independen kondisional pada Risiko (Model 1) sudah memadai.

Model yang lebih kompleks (Model 2) tidak memberikan perbaikan yang berarti.

21 Pengujian Ada Tidaknya Interaksi antara X dan Y (Homogenous model vs Conditional Association on Z)

library(MASS)

# Contoh data frame 
data <- data.frame(
  Kecamatan = factor(rep(c("Pontianak_Kota", "Pontianak_Selatan", "Pontianak_Tenggara"), each = 4)),
  Risiko = factor(rep(c("Beresiko", "Tidak_Beresiko"), times = 6)),
  Status = factor(rep(c("Status1", "Status2"), times = 6)),
  Frekuensi = c(100, 80, 90, 70, 60, 50, 40, 30, 110, 90, 85, 75)
)

# Membuat tabel kontingensi tiga arah
tbl <- xtabs(Frekuensi ~ Kecamatan + Risiko + Status, data = data)

# Model Conditional Association on Status
model_conditional <- loglm(~ Status*Kecamatan + Status*Risiko, data = tbl)

# Model Homogeneous Association
model_homog <- loglm(~ Kecamatan*Risiko + Kecamatan*Status + Risiko*Status, data = tbl)

# Uji perbandingan model dengan Likelihood Ratio Test
anova_result <- anova(model_conditional, model_homog, test = "Chisq")
print(anova_result)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Status * Kecamatan + Status * Risiko 
## Model 2:
##  ~Kecamatan * Risiko + Kecamatan * Status + Risiko * Status 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1          0  4                                    
## Model 2          0  2          0         2              1
## Saturated        0  0          0         2              1

Interpretasi :

Karena p-value = 1 (jauh > 0.05), tidak ada bukti untuk menolak hipotesis nol bahwa Model 1 sudah cukup menjelaskan data.

Penambahan interaksi antara Kecamatan dan Risiko serta interaksi Risiko dan Status dalam Model 2 tidak meningkatkan fit model secara signifikan dibanding Model 1.

Dengan kata lain, model yang mengasumsikan interaksi antara Status dengan Kecamatan dan Risiko (Model 1) sudah memadai.

Model yang lebih kompleks (Model 2) tidak memberikan perbaikan yang berarti. x
# 22. PEMILIHAN MODEL TERBAIK

library(MASS)

# Data simulasi berdasarkan data sebelumnya
data <- data.frame(
  Kecamatan = factor(rep(c("Pontianak_Kota", "Pontianak_Selatan", "Pontianak_Tenggara"), each = 4)),
  Risiko = factor(rep(c("Beresiko", "Tidak_Beresiko"), times = 6)),
  Status = factor(rep(c("Status1", "Status2"), times = 6)),
  Frekuensi = c(100, 80, 90, 70, 60, 50, 40, 30, 110, 90, 85, 75)
)

# Ubah data menjadi format data frame dengan variabel faktor dan frekuensi
tbl_df <- data

# Buat model Poisson regresi dengan interaksi antara Kecamatan dan Risiko
bestmodel <- glm(Frekuensi ~ Kecamatan + Risiko + Status + Kecamatan:Risiko,
                 family = poisson(link = "log"),
                 data = tbl_df)

# Tampilkan ringkasan model
summary(bestmodel)

## 
## Call:
## glm(formula = Frekuensi ~ Kecamatan + Risiko + Status + Kecamatan:Risiko, 
##     family = poisson(link = "log"), data = tbl_df)
## 
## Coefficients: (1 not defined because of singularities)
##                                                  Estimate Std. Error z value
## (Intercept)                                       4.55388    0.07255  62.771
## KecamatanPontianak_Selatan                       -0.64185    0.12354  -5.195
## KecamatanPontianak_Tenggara                       0.02598    0.10194   0.255
## RisikoTidak_Beresiko                             -0.23639    0.10922  -2.164
## StatusStatus2                                          NA         NA      NA
## KecamatanPontianak_Selatan:RisikoTidak_Beresiko   0.01325    0.18555   0.071
## KecamatanPontianak_Tenggara:RisikoTidak_Beresiko  0.06933    0.15205   0.456
##                                                  Pr(>|z|)    
## (Intercept)                                       < 2e-16 ***
## KecamatanPontianak_Selatan                       2.04e-07 ***
## KecamatanPontianak_Tenggara                        0.7989    
## RisikoTidak_Beresiko                               0.0304 *  
## StatusStatus2                                          NA    
## KecamatanPontianak_Selatan:RisikoTidak_Beresiko    0.9431    
## KecamatanPontianak_Tenggara:RisikoTidak_Beresiko   0.6484    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 96.338  on 11  degrees of freedom
## Residual deviance: 14.854  on  6  degrees of freedom
## AIC: 99.732
## 
## Number of Fisher Scoring iterations: 4

Interpretasi :

Variabel KecamatanPontianak_Selatan dan RisikoTidak_Beresiko berpengaruh signifikan terhadap frekuensi (count).

Variabel KecamatanPontianak_Tenggara dan interaksi Kecamatan:Risiko tidak signifikan.

Variabel Status perlu diperiksa ulang karena koefisiennya tidak terestimasi.

Model secara keseluruhan fit dengan baik (penurunan deviance besar dan residual deviance kecil).

23. CONTOH KASUS NYATA dengan DATA RILL

MODEL LOG LINEAR 2 ARAH

Model log linear dua arah adalah model statistik yang digunakan untuk menganalisis hubungan asosiasi antara dua variabel kategori dalam tabel kontingensi dua dimensi. Model ini memodelkan logaritma frekuensi yang diharapkan pada setiap sel tabel sebagai fungsi linier dari efek utama masing-masing variabel dan interaksi antar variabel.

Secara matematis, model log linear dua arah dapat dituliskan sebagai:

\[ \log(\hat{m}_{ij}) = \mu + \lambda_i^X + \lambda_j^Y + \lambda_{ij}^{XY} \]

Model log linear dua arah terdiri dari dua tipe utama:
- Model tak jenuh (independen): hanya mengandung efek utama \(\mu + \lambda_i^X + \lambda_j^Y\), tanpa interaksi.
- Model jenuh (saturated): mengandung semua efek utama dan interaksi \(\mu + \lambda_i^X + \lambda_j^Y + \lambda_{ij}^{XY}\), sehingga fit sempurna dengan data.

23.1 Contoh Kasus

Penelitian oleh Satria et al. (2023) menggunakan model log linier untuk menganalisis hubungan antara variabel kategori dalam data stunting di Kota Pontianak. Data disajikan dalam tabel kontingensi dua arah dan dianalisis menggunakan model log linier jenuh dan tak jenuh untuk menguji asosiasi antar variabel.

Data diambil dari hasil rekapitulasi Pendataan Keluarga 2021 (PK21) berupa jumlah kelurga beresiko stunting di Kota Pontianak pada tahun 2021 yang disajikan pada tabel berikut.

Kecamatan	Resiko	Tidak Beresiko
Pontianak Selatan	5238	3107
Pontianak Timur	9740	2367
Pontianak Barat	9886	2699
Pontianak Utara	14349	1274
Pontianak Kota	9195	3865
Pontianak Tenggara	3137	1567

Sumber : Badan Kependudukan dan Keluarga Berencana Nasional (BKKBN) perwakilan Kalimantan Barat

Sintaks R untuk Model Log Linier dan Model Saturated

23.2 Input Data & Buat Data Frame

library(MASS)

#Data frekuensi kasus stunting berdasarkan kecamatan
kecamatan <- c("Pontianak_Selatan", "Pontianak_Timur", "Pontianak_Barat",
"Pontianak_Utara", "Pontianak_Kota", "Pontianak_Tenggara")

resiko <- c(5238, 9740, 9886, 14349, 9195, 3137)
tidak_resiko <- c(3107, 2367, 2699, 1274, 3865, 1567)

#Membuat data frame
data_stunting <- data.frame(
Kecamatan = rep(kecamatan, each = 2),
Risiko = factor(rep(c("Resiko", "Tidak_Beresiko"), times = 6)),
Frekuensi = c(rbind(resiko, tidak_resiko))
)

print(data_stunting)

##             Kecamatan         Risiko Frekuensi
## 1   Pontianak_Selatan         Resiko      5238
## 2   Pontianak_Selatan Tidak_Beresiko      3107
## 3     Pontianak_Timur         Resiko      9740
## 4     Pontianak_Timur Tidak_Beresiko      2367
## 5     Pontianak_Barat         Resiko      9886
## 6     Pontianak_Barat Tidak_Beresiko      2699
## 7     Pontianak_Utara         Resiko     14349
## 8     Pontianak_Utara Tidak_Beresiko      1274
## 9      Pontianak_Kota         Resiko      9195
## 10     Pontianak_Kota Tidak_Beresiko      3865
## 11 Pontianak_Tenggara         Resiko      3137
## 12 Pontianak_Tenggara Tidak_Beresiko      1567

23.3 Membuat tabel kontingensi

tbl_stunting <- xtabs(Frekuensi ~ Kecamatan + Risiko, data = data_stunting)
print(tbl_stunting)

##                     Risiko
## Kecamatan            Resiko Tidak_Beresiko
##   Pontianak_Barat      9886           2699
##   Pontianak_Kota       9195           3865
##   Pontianak_Selatan    5238           3107
##   Pontianak_Tenggara   3137           1567
##   Pontianak_Timur      9740           2367
##   Pontianak_Utara     14349           1274

23.4 Model log linier Tanpa Interaksi (Independen) & Dengan Interaksi (Dependen)

#Model log linier independen (tanpa interaksi)
model_indep <- loglm(~ Kecamatan + Risiko, data = tbl_stunting)
summary(model_indep)

## Formula:
## ~Kecamatan + Risiko
## attr(,"variables")
## list(Kecamatan, Risiko)
## attr(,"factors")
##           Kecamatan Risiko
## Kecamatan         1      0
## Risiko            0      1
## attr(,"term.labels")
## [1] "Kecamatan" "Risiko"   
## 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 3918.986  5        0
## Pearson          3654.237  5        0

#Model log linier saturated (dengan interaksi)
model_sat <- loglm(~ Kecamatan * Risiko, data = tbl_stunting)
summary(model_indep)

## Formula:
## ~Kecamatan + Risiko
## attr(,"variables")
## list(Kecamatan, Risiko)
## attr(,"factors")
##           Kecamatan Risiko
## Kecamatan         1      0
## Risiko            0      1
## attr(,"term.labels")
## [1] "Kecamatan" "Risiko"   
## 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 3918.986  5        0
## Pearson          3654.237  5        0

23.5 Menghitung Odds Ratio

# Menghitung Odds Ratio per kecamatan:
# OR = (resiko / tidak_resiko) / (total resiko / total tidak_resiko) sebagai referensi keseluruhan
total_resiko <- sum(resiko)
total_tidak_resiko <- sum(tidak_resiko)
overall_odds <- total_resiko / total_tidak_resiko

df <- data.frame(kecamatan, resiko, tidak_resiko, stringsAsFactors = FALSE)
df$odds <- df$resiko / df$tidak_resiko
df$odds_ratio <- df$odds / overall_odds

df[, c("kecamatan", "odds", "odds_ratio")]

Interpretasi :

Kecamatan	Odds	Interpretasi Odds
Pontianak_Selatan	1.69	Setiap 1,69 kasus berisiko ada 1 kasus tidak berisiko. Risiko relatif sedang.
Pontianak_Timur	4.11	Setiap 4,11 kasus berisiko ada 1 kasus tidak berisiko. Risiko cukup tinggi.
Pontianak_Barat	3.66	Setiap 3,66 kasus berisiko ada 1 kasus tidak berisiko. Risiko cukup tinggi.
Pontianak_Utara	11.26	Setiap 11,26 kasus berisiko ada 1 kasus tidak berisiko. Risiko sangat tinggi.
Pontianak_Kota	2.38	Setiap 2,38 kasus berisiko ada 1 kasus tidak berisiko. Risiko sedang.
Pontianak_Tenggara	2.00	Setiap 2 kasus berisiko ada 1 kasus tidak berisiko. Risiko sedang.

Untuk menghitung persentase kenaikan odd ratio terhadap referensi (Kota Pontianak), rumusnya adalah:

(OR−1)×100%

Kecamatan	Odds Ratio	Interpretasi Odds Ratio
Pontianak_Selatan	0.49	Risiko 51% lebih rendah dari rata-rata Kota Pontianak.
Pontianak_Timur	1.19	Risiko 19% lebih tinggi dari rata-rata Kota Pontianak.
Pontianak_Barat	1.06	Risiko 6% lebih tinggi dari rata-rata Kota Pontianak.
Pontianak_Utara	3.25	Risiko 225% lebih tinggi dari rata-rata Kota Pontianak.
Pontianak_Kota	0.69	Risiko 31% lebih rendah dari rata-rata Kota Pontianak.
Pontianak_Tenggara	0.58	Risiko 42% lebih rendah dari rata-rata Kota Pontianak.

23.6 Fitting Model Log-Linear dengan R

kecamatan <- c("Pontianak_Selatan", "Pontianak_Timur", "Pontianak_Barat",
"Pontianak_Utara", "Pontianak_Kota", "Pontianak_Tenggara")

resiko <- c(5238, 9740, 9886, 14349, 9195, 3137)
tidak_resiko <- c(3107, 2367, 2699, 1274, 3865, 1567)

#Membuat tabel kontingensi
data_matrix <- matrix(c(resiko, tidak_resiko), nrow = length(kecamatan), byrow = FALSE)
rownames(data_matrix) <- kecamatan
colnames(data_matrix) <- c("Resiko", "Tidak_Resiko")

data_matrix

##                    Resiko Tidak_Resiko
## Pontianak_Selatan    5238         3107
## Pontianak_Timur      9740         2367
## Pontianak_Barat      9886         2699
## Pontianak_Utara     14349         1274
## Pontianak_Kota       9195         3865
## Pontianak_Tenggara   3137         1567

df <- data.frame(
Kecamatan = rep(kecamatan, each = 2),
Risiko = rep(c("Resiko", "Tidak_Resiko"), times = length(kecamatan)),
Frekuensi = as.vector(t(data_matrix))
)

23.7 Model tanpa interaksi: ~ Kecamatan + Risiko

#Model tanpa interaksi: ~ Kecamatan + Risiko
model_indep <- glm(Frekuensi ~ Kecamatan + Risiko, family = poisson(link = "log"), data = df)

summary(model_indep)

## 
## Call:
## glm(formula = Frekuensi ~ Kecamatan + Risiko, family = poisson(link = "log"), 
##     data = df)
## 
## Coefficients:
##                              Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)                  9.186658   0.009155 1003.511  < 2e-16 ***
## KecamatanPontianak_Kota      0.037048   0.012491    2.966  0.00302 ** 
## KecamatanPontianak_Selatan  -0.410843   0.014117  -29.103  < 2e-16 ***
## KecamatanPontianak_Tenggara -0.984092   0.017089  -57.585  < 2e-16 ***
## KecamatanPontianak_Timur    -0.038722   0.012730   -3.042  0.00235 ** 
## KecamatanPontianak_Utara     0.216239   0.011978   18.053  < 2e-16 ***
## RisikoTidak_Resiko          -1.242504   0.009306 -133.512  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 33040  on 11  degrees of freedom
## Residual deviance:  3919  on  5  degrees of freedom
## AIC: 4055.2
## 
## Number of Fisher Scoring iterations: 4

🔹 Intercept Estimate = 9.1867

Interpretasi: log dari frekuensi keluarga berisiko di Pontianak Barat adalah 9.1867

Dalam bentuk eksponensial: exp(9.1867) ≈ 9,825 → estimasi jumlah keluarga berisiko di Pontianak Barat

Kecamatan	Estimate	exp(Estimate)	Interpretasi
Pontianak_Kota	0.0370	1.038	Frekuensi 3.8% lebih tinggi dari Pontianak Barat
Pontianak_Selatan	-0.4108	0.663	Frekuensi 33.7% lebih rendah dari Pontianak Barat
Pontianak_Tenggara	-0.9841	0.374	Frekuensi 62.6% lebih rendah dari Pontianak Barat
Pontianak_Timur	-0.0387	0.962	Frekuensi 3.8% lebih rendah dari Pontianak Barat
Pontianak_Utara	0.2162	1.241	Frekuensi 24.1% lebih tinggi dari Pontianak Barat

Semua koefisien kecamatan signifikan (p < 0.01), menunjukkan bahwa ada perbedaan signifikan antar kecamatan dalam jumlah keluarga berisiko.

🔹 Risiko (Kategori)

Risiko	Estimate	exp(Estimate)	Interpretasi
Tidak_Resiko	-1.2425	0.289	Frekuensi untuk kategori Tidak Berisiko adalah 71.1% lebih rendah

⚖️ Goodness of Fit Null Deviance: 33,040 on 11 df

Residual Deviance: 3,919 on 5 df

AIC: 4055.2

Penurunan deviance menunjukkan bahwa model ini jauh lebih baik dibanding model tanpa prediktor (null model), meskipun deviasi residual masih agak besar → kemungkinan masih ada variabel penting lain atau perlu mempertimbangkan interaksi.

23.8 Model dengan interaksi: ~ Kecamatan * Risiko

#Model dengan interaksi: ~ Kecamatan * Risiko
model_sat <- glm(Frekuensi ~ Kecamatan * Risiko, family = poisson(link = "log"), data = df)

summary(model_sat)

## 
## Call:
## glm(formula = Frekuensi ~ Kecamatan * Risiko, family = poisson(link = "log"), 
##     data = df)
## 
## Coefficients:
##                                                Estimate Std. Error z value
## (Intercept)                                     9.19887    0.01006 914.629
## KecamatanPontianak_Kota                        -0.07246    0.01449  -5.001
## KecamatanPontianak_Selatan                     -0.63518    0.01709 -37.167
## KecamatanPontianak_Tenggara                    -1.14785    0.02049 -56.014
## KecamatanPontianak_Timur                       -0.01488    0.01428  -1.042
## KecamatanPontianak_Utara                        0.37256    0.01307  28.503
## RisikoTidak_Resiko                             -1.29824    0.02172 -59.778
## KecamatanPontianak_Kota:RisikoTidak_Resiko      0.43154    0.02897  14.897
## KecamatanPontianak_Selatan:RisikoTidak_Resiko   0.77596    0.03138  24.731
## KecamatanPontianak_Tenggara:RisikoTidak_Resiko  0.60413    0.03780  15.984
## KecamatanPontianak_Timur:RisikoTidak_Resiko    -0.11638    0.03157  -3.686
## KecamatanPontianak_Utara:RisikoTidak_Resiko    -1.12328    0.03642 -30.844
##                                                Pr(>|z|)    
## (Intercept)                                     < 2e-16 ***
## KecamatanPontianak_Kota                        5.69e-07 ***
## KecamatanPontianak_Selatan                      < 2e-16 ***
## KecamatanPontianak_Tenggara                     < 2e-16 ***
## KecamatanPontianak_Timur                       0.297339    
## KecamatanPontianak_Utara                        < 2e-16 ***
## RisikoTidak_Resiko                              < 2e-16 ***
## KecamatanPontianak_Kota:RisikoTidak_Resiko      < 2e-16 ***
## KecamatanPontianak_Selatan:RisikoTidak_Resiko   < 2e-16 ***
## KecamatanPontianak_Tenggara:RisikoTidak_Resiko  < 2e-16 ***
## KecamatanPontianak_Timur:RisikoTidak_Resiko    0.000228 ***
## KecamatanPontianak_Utara:RisikoTidak_Resiko     < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3.3040e+04  on 11  degrees of freedom
## Residual deviance: 1.1844e-12  on  0  degrees of freedom
## AIC: 146.23
## 
## Number of Fisher Scoring iterations: 2

🔹 (Intercept) = 9.1989

Log dari frekuensi keluarga berisiko di Pontianak Barat

exp(9.1989) ≈ 9,899 → estimasi jumlah keluarga berisiko di Pontianak Barat

Kecamatan	Estimate	exp(Estimate)	Interpretasi
Pontianak Kota	-0.0725	0.930	Frekuensi 7% lebih rendah dari Pontianak Barat (kategori berisiko)
Pontianak Selatan	-0.6352	0.530	Frekuensi 47% lebih rendah dari Pontianak Barat
Pontianak Tenggara	-1.1479	0.317	Frekuensi 68% lebih rendah
Pontianak Timur	-0.0149	0.985	Hampir sama
Pontianak Utara	0.3726	1.451	Frekuensi 45% lebih tinggi dibanding Pontianak Barat

🔹 Efek Utama Risiko

Risiko	Estimate	exp(Estimate)	Interpretasi
Tidak Berisiko	-1.2982	0.273	Di Pontianak Barat, jumlah keluarga tidak berisiko 72.7% lebih rendah

🔹 Interaksi Kecamatan × Risiko

Interaksi	Estimate	exp(Estimate)	Interpretasi
Kota × Tidak_Resiko	0.4315	1.539	Efek tambahan: frekuensi 53.9% lebih tinggi dari Pontianak Barat
Selatan × Tidak_Resiko	0.7760	2.173	Efek tambahan: frekuensi 117% lebih tinggi
Tenggara × Tidak_Resiko	0.6041	1.829	Efek tambahan: frekuensi 82.9% lebih tinggi
Timur × Tidak_Resiko	-0.1164	0.890	Efek tambahan: frekuensi 11% lebih rendah
Utara × Tidak_Resiko	-1.1233	0.325	Efek tambahan: frekuensi 67.5% lebih rendah

⚖️ Goodness of Fit Null Deviance: 33,040 (df = 11)

Residual Deviance: ~0 (df = 0)

AIC: 146.23 (sangat kecil)

🔎 Interpretasi:

Model fit sangat baik terhadap data (residual deviasi nyaris nol)

Menunjukkan bahwa interaksi sangat penting dan harus disertakan dalam model

24. Referensi

• Dobson, A. J. (1983). Introduction to statistical modelling. Springer US.

• Breiman, L. (2001). Random Forests. Machine Learning

• Agresti, A. (2007). An Introduction to Categorical Data Analysis (2nd ed.). Wiley.

• Greenacre, M. (2007). Correspondence Analysis in Practice. CRC Press.

• Agresti, A. (2013). Categorical Data Analysis (3rd ed.). Wiley.

• Nugraha, J. (2016). Pengantar Analisis Data Kategorik: Metode dan Aplikasi Menggunakan Program R. Deepublish.

• Agresti, A. (2018). An Introduction to Categorical Data Analysis (2nd ed.). Wiley

• Rudas, T. (2018). Lectures on Categorical Data Analysis. Springer.

• Astuti, A. B., Effendi, A., Astutik, S., & Sumarminingsih, E. (2020). Analisis Data Kategorik Menggunakan R: Teori dan Aplikasinya di Berbagai Bidang. UB Press.

• Mangiafico, S. S. (2025). rcompanion: Functions to Support Extension Education Program Evaluation. Rutgers Cooperative Extension, New Brunswick, New Jersey. Version 2.5.0.

• CRAN Project. (2025). Package rcompanion.

• Hosmer, D.W., Lemeshow, S., & Sturdivant, R.X. (2013). Applied Logistic Regression (3rd ed.). Wiley.

• Agresti, A. (2018). Statistical Methods for the Social Sciences (5th ed.). Pearson.

• James, G., Witten, D., Hastie, T., & Tibshirani, R. (2013). An Introduction to Statistical Learning. Springer.

• Saito, T., & Rehmsmeier, M. (2015). The Precision-Recall Plot Is More Informative than the ROC Plot When Evaluating Binary Classifiers on Imbalanced Datasets. PLOS ONE.

• Davis, J., & Goadrich, M. (2006). The Relationship Between Precision-Recall and ROC Curves. ICML.

• Google Developers. (2025). Machine Learning Crash Course: ROC and AUC.

• Staf ULM. (2017). Menghitung Kinerja Algoritma Klasifikasi: Pilih ROC Curve atau Precision-Recall Curve?

• MySkill Blog. (2024). Memahami Receiver Operating Characteristic (ROC) Curve dalam Data Analysis.

• Halo Ryan Blog. (2024). ROC-AUC: Pengertian, Fungsi, dan Cara Menggunakannya dalam Machine Learning.