Kata Pengantar
1 Pendahuluan
2 Metode dalam Analisis Kategori
3 Distribusi Probabilitas dalam Data Kategori
4 Desain Sampling dalam Analisis Data Kategori
- 4.1 Prostective Sampling
  - 4.1.1 Eksperimen
  - 4.1.2 Studi Kohort
- 4.2 Retrospective Sampling
  - 4.2.1 Studi Kasus-Kontrol
  - 4.2.2 Studi Kohort Retrospektif
5 Tabel Kontingensi 2 x 2
- 5.1 Distribusi Peluang dalam Tabel Kontingensi 2x2
- 5.2 Ukuran Asosiasi dalam Kategori 2x2
6 Inferensi Tabel Kontingensi Dua Arah
7 Tabel Kontingensi Tiga Arah
8 Generalized Linear Model (GLM)
9 Inferensi GLM
10 Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio Regresi logistik
11 Pemilihan Model Regresi Logistik dan Evaluasi
12 Regresi Logistik Multinomial
13 Regresi Logistik Ordinal
14 Log Linear Model
15 Model Log Linear 2 Arah
16 Model Log-Linear Tiga Arah
Referensi

Kata Pengantar

Puji dan syukur penulis panjatkan ke hadirat Tuhan Yang Maha Esa karena atas limpahan rahmat dan karunia-Nya, penulis dapat menyelesaikan penulisan eBook ini yang berjudul “Analisis Data Kategori”. eBook ini merupakan hasil dari proses pembelajaran dan pendalaman materi selama mengikuti perkuliahan Analisis Data Kategori yang diampu oleh Dr. I Gede Nyoman Mindra Jaya, dosen Program Studi Statistika, Universitas Padjadjaran.

Selama proses perkuliahan, penulis memperoleh banyak wawasan dan pemahaman mendalam mengenai prinsip-prinsip dasar hingga lanjutan dalam analisis data kategorik. Materi yang dikaji mencakup teori probabilitas diskrit, tabel kontingensi, model log-linear, regresi logistik biner, multinomial, dan ordinal, hingga berbagai uji asosiasi yang relevan untuk data-data kategori.

Adapun tujuan utama dari penulisan eBook ini adalah untuk merangkum, menyusun ulang, dan menyajikan seluruh materi yang telah dipelajari dalam format yang lebih sistematis, aplikatif, dan mudah dipahami. Penulis berharap eBook ini dapat menjadi sumber pembelajaran yang bermanfaat, khususnya bagi mahasiswa, peneliti, maupun praktisi yang ingin mendalami analisis data kategori. Pembahasan dalam buku ini dilengkapi dengan dasar teori, rumus matematis, contoh kasus nyata, serta interpretasi hasil analisis menggunakan perangkat lunak statistik.

Penulis menyadari bahwa karya ini masih jauh dari sempurna, baik dari segi isi maupun penyajian. Oleh karena itu, saran dan kritik yang membangun sangat penulis harapkan demi perbaikan dan penyempurnaan pada edisi-edisi berikutnya.

Akhir kata, penulis menyampaikan terima kasih yang sebesar-besarnya kepada Dr. I Gede Nyoman Mindra Jaya atas ilmu, bimbingan, serta inspirasi akademik yang telah diberikan selama proses perkuliahan. Semoga eBook ini dapat memberikan manfaat dan kontribusi positif dalam pengembangan ilmu statistik, khususnya di bidang analisis data kategori.

1 Pendahuluan

Dalam statistik, data kategori merujuk pada data yang terbagi ke dalam kelompok atau kategori tertentu, bukan berupa angka dengan nilai kontinu. Jenis data ini sering digunakan di berbagai bidang untuk menganalisis pola, hubungan, dan tren yang tidak bisa diukur langsung dengan skala numerik. Contoh umum data kategori dapat berupa status mengidap penyakit (sakit/tidak sakit), golongan darah (A/B/AB/O), dan metrik kepuasan (puas/tidak puas).

Berbeda dengan data kontinu, analisis kategori perlu pendekatan khusus seperti tabel kontingensi, uji chi-square, analisis regresi logistik, dan metode estimasi berbasis model probabilitas. Analisis data kategori juga mulai banyak digunakan pada machine learning dan big data analytics.

1.1 Tujuan Analisis Data Kategori

Analisis data kategori mempunyai berbagai tujuan yang membantu dalam pengambilan keputusan berbasis data. Beberapa tujuannya yaitu :

Mengidentifikasi Pola dan Tren

Misalnya pada sektor ekonomi, analisis ini dapat membantu memaham tren konsumsi yang sedang ada di masyarakat
Menganalisis Hubungan Antarvariabel

Misalnya dalam penelitian olahraga, analisis data kategori dapat membantu memahami hubungan efek latihan tertentu dan performa atlet
Membantu dalam Pengambilan Keputusan

Misalnya pada perusahaan dapat digunakan untuk menentukan kebijakan pengeluaran produk baru
Mengembangkan Model Prediktif

contohnya digunakan pada regresi logistik untuk memprediksi kejadian biner (misalnya memprediksi kemungkinan seseorang gagal membayar pinjaman bank)

1.2 Ruang Lingkup dan Jenis Data

data kategori berbeda dari data kuantitatif, perbedaan utamanya data kategorik tidak dapat diukur dengan skala numerik kontinu. Dibandingkan data kuantitatif yang berbasis angka, data kategori lebih menekankan pada klasifikasi datanya.

data yang digunakan pada data kategori adalah :

Nominal : data yang tidak memiliki urutan/hierarki (contoh : golongan darah, asal negara)
Ordinal : data ini memiliki hierarki atau urutannya (contoh : tingkat pendidikan, tingkat prioritas)

selain itu juga ada data biner (dua kategori) atau multikategori (lebih dari dua).

1.3 Manfaat Luas Analisis Data Kategori

beberapa contoh analisis kategori dalam berbagai bidang :

pendidikan : membantu memahami efektivitas pembelajaran terhadap nilai siswa
kriminalitas : membantu memetakan pola kejahatan yang sering terjadi di suatu tempat
kesehatan : studi mengenai penyakit, memahami penyebab dan membantu cara pengobatan
olahraga : memahami riawayat cedera atlet dan mencegah penyebab

2 Metode dalam Analisis Kategori

analisis data kategori memiliki berbagai macam metode yang dapat digunakan, berikut ini beberapa metode yang umum dipakai :

Tabel Kontingensi dan Uji Chi-Square

tujuan dari metode ini adalah untuk melihat hubungan antar dua variabel kategori

contoh : apakah rajin berolahraga berpengaruh terhadap kondisi tubuh
Regresi Logistik

digunakan untuk memprediksi suatu kejadian berdasarkan data kategori, metode ini juga cocok digunakan dalam prediksi risiko, klasifikasi, dan pengambilan keputusan

contoh : memprediksi calon presiden pilihan seseorang berdasarkan pekerjaan, preferensi politik, dan status ekonomi
Correspondence Analysis (CA)

metode ini digunakan untuk memetakan hubungan kompleks pada data multikategori

contoh : analisis jumlah pendapatan terhadap kebiasaan belanja
Decision tree dan Random Forest

metode berbasis machine learning ini berguna untuk klasifikasi dan prediksi berdasarkan data kategori

contoh : mengklasifikasikan risiko siswa putus sekolah berdasarkan penghasilan orang tua dan jenis sekolah

3 Distribusi Probabilitas dalam Data Kategori

Dalam analisis data kategori, distribusi probabilitas digunakan untuk memodelkan kemungkinan terjadinya berbagai hasil dari variabel kategori. Berikut adalah beberapa distribusi probabilitas yang umum digunakan.

3.1 Distribusi Bernoulli

Distribusi Bernoulli digunakan ketika percobaan hanya memiliki dua hasil: sukses (1) atau gagal (0).

Fungsi Probabilitas:

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

\(X\): variabel acak (0 atau 1)
\(p\): probabilitas sukses

Contoh:

misalnya, kita ingin mensimulasikan hasil ujian (Lulus = 1, Gagal = 0) dengan probabilitas lulus = 0.6.

set.seed(789) 
bernoulli_sample <- rbinom(n = 15, size = 1, prob = 0.6) 
bernoulli_sample

##  [1] 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1

3.2 Distribusi Binomial

Distribusi Binomial digunakan untuk n percobaan independen Bernoulli.

Fungsi Probabilitas:

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

\(X\): Jumlah keberhasilan dalam \(n\) percobaan
\(n\): Jumlah percobaan
\(k\): Jumlah sukses
\(p\): Probabilitas sukses

Contoh:

Jumlah keberhasilan dalam 10 kali lemparan koin

set.seed(789)
binomial_sample <- rbinom(n = 10, size = 5, prob = 0.5)
binomial_sample

##  [1] 3 1 0 3 2 0 3 1 2 2

3.3 Distribusi Multinomial

Distribusi Multinomial merupakan generalisasi dari binomial untuk kasus lebih dari dua kategori.

Fungsi Probabilitas:

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

Contoh:

Simulasi 25 suara pemilih terhadap 4 kandidat: A, B, C, dan D dengan proporsi 0.25, 0.35, 0.20, dan 0.20.

set.seed(789)
votes <- rmultinom(n = 1, size = 25, prob = c(0.25, 0.35, 0.20, 0.20))
rownames(votes) <- c("Kandidat A", "Kandidat B", "Kandidat C", "Kandidat D")
votes

##            [,1]
## Kandidat A    7
## Kandidat B    6
## Kandidat C    2
## Kandidat D   10

3.4 Distribusi Poisson

Distribusi Poisson digunakan untuk menghitung jumlah kejadian dalam interval waktu atau ruang tertentu.

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

Contoh:

Simulasi jumlah kendaraan melintasi jalan tol dalam 10 menit (λ = 5):

set.seed(789)
trafik <- rpois(n = 10, lambda = 5)
trafik

##  [1] 6 2 1 5 5 1 5 3 4 4

4 Desain Sampling dalam Analisis Data Kategori

Dalam analisis data kategori, desain sampling sangat menentukan validitas dan reliabilitas hasil penelitian. Pemilihan desain yang tepat bergantung pada tujuan penelitian dan jenis data yang dikumpulkan. Umumnya, desain ini diklasifikasikan menjadi dua pendekatan utama:

Prospective sampling
Retrospective sampling

Setiap pendekatan memiliki karakteristik dan metode pengambilan sampel yang berbeda, serta biasa digunakan dalam studi eksperimental maupun observasional seperti eksperimen, studi kohort, dan studi kasus-kontrol.

4.1 Prostective Sampling

Prospective sampling merupakan metode pengambilan data secara ke depan, di mana subjek diamati dalam jangka waktu tertentu untuk melihat perkembangan variabel yang diteliti.

4.1.1 Eksperimen

Eksperimen melibatkan alokasi acak subjek ke kelompok perlakuan dan kontrol. Beberapa teknik sampling yang digunakan antara lain:

Simple Random Sampling (SRS): Semua individu memiliki peluang yang sama untuk dipilih.
Stratified Sampling: Populasi dibagi menjadi strata, lalu diambil sampel acak dari setiap strata.
Cluster Sampling: Populasi dibagi ke dalam klaster, dan beberapa klaster dipilih secara acak.

4.1.2 Studi Kohort

Studi kohort adalah penelitian observasional di mana individu dengan karakteristik tertentu diikuti dari waktu ke waktu. Teknik sampling yang digunakan antara lain:

Census Sampling: Seluruh populasi dijadikan sampel.
Systematic Sampling: Sampel dipilih berdasarkan interval tertentu.
Matched Sampling: Subjek dipasangkan berdasarkan karakteristik serupa.

4.2 Retrospective Sampling

Retrospective sampling menggunakan data historis atau peristiwa masa lalu untuk menganalisis hubungan antara faktor risiko dan outcome.

4.2.1 Studi Kasus-Kontrol

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

Purposive Sampling: Pemilihan berdasarkan karakteristik relevan.
Snowball Sampling: Responden awal membantu merekrut responden lain.
Incidence Density Sampling: Mempertimbangkan waktu kejadian kasus dalam pemilihan sampel.

4.2.2 Studi Kohort Retrospektif

Studi ini menggunakan data historis untuk mengelompokkan individu berdasarkan paparan, kemudian menelusuri outcome-nya. Teknik sampling meliputi:

Convenience Sampling: Berdasarkan data yang tersedia.
Quota Sampling: Sampel dipilih berdasarkan proporsi tertentu.
Case-Based Sampling: Berdasarkan karakteristik kasus.

Berikut adalah perbandingan berbagai jenis desain sampling yang umum digunakan dalam analisis data kategori:

Perbandingan Desain Sampling dalam Analisis Data Kategori
Jenis.Studi	Pendekatan	Metode.Sampling	Keuntungan	Kelemahan
Eksperimen	Prospektif	SRS, Stratified, Cluster	Kontrol tinggi terhadap variabel, dapat analisis sebab-akibat	Biaya tinggi, etika dan validitas perlu diperhatikan
Studi Kohort	Prospektif	Census, Systematic, Matched	Mengamati perkembangan kejadian jangka panjang	Waktu lama, potensi kehilangan partisipan
Studi Kasus-Kontrol	Retrospektif	Purposive, Snowball, Incidence Density	Cepat dan efisien untuk kasus langka	Sulit kontrol variabel pengganggu, bias recall
Studi Kohort Retrospektif	Retrospektif	Convenience, Quota, Case-Based	Memanfaatkan data historis, biaya rendah	Bergantung pada kualitas data historis

Pemilihan desain sampling yang sesuai merupakan langkah krusial dalam analisis data kategori. Setiap pendekatan memiliki kelebihan dan kekurangan yang perlu dipertimbangkan sesuai dengan tujuan penelitian dan sumber daya yang tersedia.

5 Tabel Kontingensi 2 x 2

Tabel kontingensi 2×2 merupakan bentuk paling dasar dari tabel kontingensi yang biasa digunakan dalam analisis statistik. Fungsinya adalah untuk mengevaluasi adanya hubungan antara dua variabel yang bersifat kategorik. Dalam praktiknya, tabel ini sering digunakan untuk mengetahui apakah terdapat hubungan atau asosiasi antara dua hal, misalnya antara metode pengobatan dengan tingkat kesembuhan pasien, atau antara kebiasaan merokok dengan kemungkinan terkena kanker paru-paru. Karena kesederhanaannya, tabel 2×2 menjadi alat yang umum digunakan dalam berbagai studi statistik untuk menganalisis keterkaitan antar variabel kategori.

Tabel kontingensi 2×2 memiliki struktur sebagai berikut:

	Kategori 1 (+)	Kategori 2 (−)	Total
Grup 1	\(n_{11}\)	\(n_{12}\)	\(n_{1.}\)
Grup 2	\(n_{21}\)	\(n_{22}\)	\(n_{2.}\)
Total	\(n_{.1}\)	\(n_{.2}\)	\(n\)

Notasi dan Penjelasan

\(n_{11}\) = Jumlah observasi dalam Grup 1 dan Kategori 1
\(n_{12}\) = Jumlah observasi dalam Grup 1 dan Kategori 2
\(n_{21}\) = Jumlah observasi dalam Grup 2 dan Kategori 1
\(n_{22}\) = Jumlah observasi dalam Grup 2 dan Kategori 2
\(n_{1.} = n_{11} + n_{12}\) = Total observasi di Grup 1
\(n_{2.} = n_{21} + n_{22}\) = Total observasi di Grup 2
\(n_{.1} = n_{11} + n_{21}\) = Total observasi dalam Kategori 1
\(n_{.2} = n_{12} + n_{22}\) = Total observasi dalam Kategori 2
\(n = n_{11} + n_{12} + n_{21} + n_{22}\) = Total seluruh observasi

contoh tabel hubungan jenis kelamin dengan jenis bacaan:

	Fiksi Ilmiah	Novel	Total
Laki-Laki	38	25	63
Perempuan	20	15	35
Total	58	40	98

dengan r:

# Data Observasi
data <- matrix(c(38, 25, 20, 15), nrow = 2, byrow = TRUE)
colnames(data) <- c("Fiksi Ilmiah", "Novel")
rownames(data) <- c("Laki-laki", "Perempuan")
n <- sum(data)
data

##           Fiksi Ilmiah Novel
## Laki-laki           38    25
## Perempuan           20    15

5.1 Distribusi Peluang dalam Tabel Kontingensi 2x2

5.1.1 Peluang Bersama

Peluang bersama (\(P_{ij}\)) menunjukkan probabilitas kedua variabel berada dalam suatu sel tabel kontingensi secara bersamaan.

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

Hitung peluang bersama secara manual:

\(P(\text{Laki-laki, Fiksi Ilmiah}) = \frac{38}{98} = 0.39\) \(P(\text{Laki-laki, Novel}) = \frac{25}{98} = 0.26\) \(P(\text{Perempuan, Fiksi Ilmiah}) = \frac{20}{98} = 0.20\) \(P(\text{Perempuan, Novel}) = \frac{15}{98} = 0.15\)

dengan r:

# Peluang Bersama
P_joint <- data / n
list(Peluang_Bersama = P_joint)

## $Peluang_Bersama
##           Fiksi Ilmiah     Novel
## Laki-laki    0.3877551 0.2551020
## Perempuan    0.2040816 0.1530612

5.1.2 Peluang Marginal

Peluang marginal adalah peluang terjadinya suatu peristiwa tanpa memperhatikan peristiwa lainnya

Rumus	Keterangan
\(P(A_i) = \frac{n_{i.}}{n}\)	Peluang marginal baris
\(P(B_j) = \frac{n_{.j}}{n}\)	Peluang marginal kolom

Hitung peluang marginal secara manual:

\(P(\text{Laki-laki}) = \frac{63}{98} = 0.64\) \(P(\text{Perempuan}) = \frac{35}{98} = 0.36\) \(P(\text{Fiksi Ilmiah}) = \frac{58}{98} = 0.59\) \(P(\text{Novel}) = \frac{40}{98} = 0.41\)

dengan r:

# Peluang Marginal
P_marginal_rows <- rowSums(data) / n
P_marginal_cols <- colSums(data) / n
list(Peluang_Marginal_Baris = P_marginal_rows, Peluang_Marginal_Kolom = P_marginal_cols)

## $Peluang_Marginal_Baris
## Laki-laki Perempuan 
## 0.6428571 0.3571429 
## 
## $Peluang_Marginal_Kolom
## Fiksi Ilmiah        Novel 
##    0.5918367    0.4081633

5.1.3 Peluang Bersyarat

Peluang bersyarat adalah peluang terjadinya suatu peristiwa dengan syarat bahwa peristiwa lain sudah terjadi.

\[ P(B_j|A_i)=\frac{n_{ij}}{n_{i.}} \] Hitung peluang bersyarat secara manual:

\(P(\text{Fiksi Ilmiah|Laki-laki}) = \frac{38}{63} = 0.60\) \(P(\text{Fiksi Ilmiah|Perempuan}) = \frac{20}{35} = 0.40\) \(P(\text{Novel|Laki-laki}) = \frac{25}{63} = 0.57\) \(P(\text{Novel|Perempuan}) = \frac{15}{35} = 0.43\)

dengan r:

# Peluang Bersyarat
P_conditional <- data / rowSums(data)
list(Peluang_Bersyarat = P_conditional)

## $Peluang_Bersyarat
##           Fiksi Ilmiah     Novel
## Laki-laki    0.6031746 0.3968254
## Perempuan    0.5714286 0.4285714

Interpretasi

Peluang bersama menggambarkan kemungkinan terjadinya dua kejadian sekaligus dalam satu sel pada tabel kontingensi.
Peluang marginal merepresentasikan probabilitas dari satu kejadian saja, tanpa mempertimbangkan kejadian lainnya.
Peluang bersyarat menunjukkan bagaimana peluang suatu kejadian berubah ketika diketahui informasi atau kondisi dari variabel lain.

Jika \(P(\text{Fiksi Ilmiah} \mid \text{Laki-laki}) > P(\text{Fiksi Ilmiah} \mid \text{Perempuan})\) maka ini menunjukkan bahwa Laki-laki lebih memilih bacaan fiksi ilmiah dibanding perempuan

5.2 Ukuran Asosiasi dalam Kategori 2x2

Dalam analisis statistik, tabel kontingensi 2 × 2 digunakan untuk mengevaluasi apakah terdapat keterkaitan antara dua variabel kategori. Tabel ini menyajikan frekuensi gabungan dari dua variabel yang masing-masing memiliki dua kategori, sehingga memudahkan kita untuk mengamati dan mengukur adanya hubungan statistik di antara keduanya.

Risk Diference (RD): Mengukur selisih risiko antara dua kelompok.
Relative Risk (RR): Membandingkan besarnya risiko kejadian di satu kelompok terhadap kelompok lainnya.
Odds Ratio (OR): Membandingkan peluang kejadian antara dua kelompok.
Uji Chi-Square dan Uji Fisher: Digunakan untuk menguji apakah hubungan antara dua variabel tersebut bersifat signifikan secara statistik.

5.2.1 Risk Difference

Risk Difference (RD) atau Perbedaan Risiko adalah ukuran yang digunakan untuk mengetahui seberapa besar perbedaan risiko suatu kejadian terjadi antara dua kelompok.

Biasanya digunakan dalam bidang kesehatan, terutama untuk membandingkan risiko suatu penyakit atau kondisi antara kelompok yang terpapar suatu faktor (misalnya merokok, obat, atau vaksin) dan kelompok yang tidak terpapar.

RD dihitung dengan cara mengurangkan proporsi kejadian di kelompok pertama dengan proporsi kejadian di kelompok kedua. Hasilnya bisa positif, negatif, atau nol

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

Jika RD > 0, maka risiko kejadian lebih tinggi di Grup 1 dibandingkan Grup 2.
Jika RD < 0, maka risiko kejadian lebih rendahdi Grup 1 dibandingkan Grup 2.
Jika RD = 0, maka tidak ada perbedaan risiko antara dua kelompok.

RD <- function(n11, n12, n21, n22) {
(n11 / (n11 + n12)) - (n21 / (n21 + n22))
}
RD(38, 25, 20, 15)

## [1] 0.03174603

5.2.2 Relative Risk

Relative Risk (RR) atau Risiko Relatif adalah ukuran statistik yang digunakan untuk membandingkan risiko terjadinya suatu peristiwa (misalnya penyakit, keberhasilan pengobatan, dll) antara dua kelompok, yaitu kelompok yang terpapar suatu faktor dan kelompok yang tidak terpapar. RR menunjukkan “Seberapa besar kemungkinan kejadian terjadi pada kelompok yang terpapar dibandingkan dengan yang tidak terpapar?”

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

Jika RR > 1, maka kejadian lebih sering terjadi di Grup 1 dibandingkan Grup 2.
Jika RR < 1, maka kejadian lebih jarang terjadi di Grup 1 dibandingkan Grup 2.
Jika RR = 1, maka tidak ada perbedaan risiko antara dua kelompok.

RR <- function(n11, n12, n21, n22) {
(n11 / (n11 + n12)) / (n21 / (n21 + n22))
}
RR(38, 25, 20, 15)

## [1] 1.055556

5.2.3 Odds Ratio

Odds Ratio (OR) atau Rasio Odds adalah ukuran statistik yang digunakan untuk membandingkan peluang terjadinya suatu kejadian antara dua kelompok, yaitu kelompok yang terpapar dan yang tidak terpapar terhadap suatu faktor.

OR dihitung dengan membandingkan perbandingan antara kejadian dan tidak kejadian di satu kelompok dengan perbandingan yang sama di kelompok lainnya. Dengan kata lain, OR melihat seberapa besar kemungkinan suatu kejadian terjadi pada kelompok terpapar dibandingkan dengan kelompok tidak terpapar. OR sering digunakan dalam studi case-control tetapi juga dapat digunakan dalam studi kohort dan eksperimental.

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

Jika OR > 1, maka peluang kejadian lebih besar di Grup 1 dibandingkan Grup 2.
Jika OR < 1, maka peluang kejadian lebih kecil di Grup 1 dibandingkan Grup 2.
Jika OR = 1, maka tidak ada perbedaan peluang kejadian antara dua kelompok.

OR <- function(n11, n12, n21, n22) {
(n11 * n22) / (n12 * n21)
}
OR(38, 25, 20, 15)

## [1] 1.14

Perbandingan antara RD, RR, dan OR

Ukuran	Definisi	Desain Sampling yang Cocok	Interpretasi
RD	Perbedaan peluang terjadinya suatu peristiwa antara dua grup	Cocok untuk studi kohort dan eksperimen acak	Menggambarkan selisih risiko secara langsung antara dua kelompok
RR	Rasio perbandingan probabilitas kejadian pada dua kelompok	Umumnya digunakan pada studi kohort atau klinis	Menunjukkan seberapa besar kemungkinan suatu kejadian di satu kelompok
OR	Rasio antara peluang kejadian pada satu grup dibandingkan grup lainnya	Relevan untuk studi kasus-kontrol dan observasional	Mengukur tingkat kekuatan hubungan antara dua variabel kategori

Kesimpulan

RD digunakan untuk memahami dampak absolut dari suatu faktor risiko atau intervensi.
RR lebih cocok untuk studi kohort atau eksperimen karena mengukur kemungkinan relatif.
OR seringdigunakan dalam studi kasus-kontrol karena dapat memperkirakan risiko relatif dalam desain penelitian ini.

6 Inferensi Tabel Kontingensi Dua Arah

Inferensi statistik pada tabel kontingensi dua arah bertujuan untuk membuat kesimpulan mengenai hubungan antara dua variabel kategorikal berdasarkan data sampel. Dalam konteks ini, tabel kontingensi menyajikan distribusi frekuensi dari dua variabel kategorikal dalam format matriks.

Inferensi dibagi menjadi dua jenis utama:

Estimasi parameter populasi
Pengujian hipotesis

6.1 Estimasi

6.1.1 Estimasi Titik

stimasi titik bertujuan untuk memberikan satu nilai yang dianggap sebagai perkiraan terbaik dari parameter populasi. Misalnya, proporsi suatu kategori dihitung sebagai:

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

Dimana:

\(\hat{p}\) adalah estimasi proporsi
\(x\) adalah jumlah kasus pada kategori tertentu
\(n\) adalah total observasi

6.1.2 Estimasi Interval

Estimasi interval memberikan rentang nilai yang diyakini mencakup parameter populasi dengan tingkat kepercayaan tertentu. Rumus untuk estimasi interval proporsi adalah:

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

Dimana:

\(\hat{p}\): estimasi titik proporsi
\(n\): ukuran sampel
\(Z_{\alpha/2}\): nilai kritis distribusi normal standar untuk tingkat kepercayaan tertentu (misalnya 1.96 untuk 95%)

6.2 Uji Hipotesis

6.2.1 Uji Proporsi

Uji proporsi dua sampel digunakan untuk menguji apakah dua proporsi populasi berbeda secara signifikan berdasarkan data sampel.

Struktur Tabel Kontingensi 2x2:

	Kejadian (+)	Tidak Kejadian (-)	Total
Grup 1	\(n_{11}\)	\(n_{12}\)	\(n_{1.}\)
Grup 2	\(n_{21}\)	\(n_{22}\)	\(n_{2.}\)
Total	\(n_{.1}\)	\(n_{.2}\)	\(n\)

Hipotesis:

H₀: \(p_1 = p_2\) (tidak ada perbedaan proporsi antara dua kelompok)
H₁: \(p_1 \neq p_2\) (terdapat perbedaan proporsi)

Proporsi masing-masing kelompok:

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

Proporsi gabungan:

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

Statistik uji Z:

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

Jika \(|Z|\) > nilai kritis (misalnya 1.96 untuk \(\alpha = 0.05\)), maka hipotesis nol ditolak.

Misalkan: - \(n_{11} = 50\), \(n_{12} = 30\) → \(n_{1.} = 80\) - \(n_{21} = 30\), \(n_{22} = 50\) → \(n_{2.} = 80\)

Maka:

\[ \hat{p}_1 = \frac{50}{80} = 0.625, \quad \hat{p}_2 = \frac{30}{80} = 0.375 \]

\[ \hat{p} = \frac{50 + 30}{160} = 0.5 \]

\[ Z = \frac{0.625 - 0.375}{\sqrt{0.5(1 - 0.5)(\frac{1}{80} + \frac{1}{80})}} = \frac{0.25}{\sqrt{0.5 \cdot 0.5 \cdot 0.025}} = \frac{0.25}{0.0791} \approx 3.16 \]

Implementasi di R

# Data
n11 <- 50; n12 <- 30
n21 <- 30; n22 <- 50
n1_ <- n11 + n12
n2_ <- n21 + n22

# Estimasi proporsi
p1 <- n11 / n1_
p2 <- n21 / n2_
p_pooled <- (n11 + n21) / (n1_ + n2_)
se <- sqrt(p_pooled * (1 - p_pooled) * (1/n1_ + 1/n2_))
z <- (p1 - p2) / se

list(Z_statistic = z)

## $Z_statistic
## [1] 3.162278

Jika nilai Z > 1.96 (untuk \(\alpha = 0.05\)), maka kita menolak hipotesis nol. Ini menunjukkan terdapat perbedaan proporsi yang signifikan secara statistik antara dua kelompok.

6.2.2 Uji Asosiasi

Uji asosiasi dalam konteks tabel kontingensi 2x2 bertujuan untuk mengukur hubungan antara dua variabel kategorikal. Tiga ukuran asosiasi yang umum digunakan adalah Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR).

Hipotesis Umum

H₀ (Hipotesis Nol): Tidak ada hubungan (asosiasi) antara dua variabel.
H₁ (Hipotesis Alternatif): Terdapat asosiasi antara dua variabel.

6.2.2.1 Risk Difference (RD)

Risk Difference mengukur perbedaan absolut antara probabilitas kejadian pada dua kelompok.

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

Standard Error (SE):

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

Statistik uji Z:

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

6.2.2.2 Relative Risk (RR)

Relative Risk menunjukkan rasio risiko antara dua kelompok:

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

Standard Error untuk log(RR):

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

Statistik uji Z:

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

6.2.2.3 Odds Ratio (OR)

Odds Ratio mengukur perbandingan odds antara dua kelompok:

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

Contoh Perhitungan Manual

Misal:

	Fiksi Ilmiah	Novel	Total
Laki-Laki	38	25	63
Perempuan	20	15	35
Total	58	40	98

Maka:

\(\hat{p}_1 = 0.603\), \(\hat{p}_2 = 0.571\)

Hitung RD, RR, dan OR beserta nilai Z masing-masing.

implementasi pada R

n11 <- 38; n12 <- 25; n21 <- 20; n22 <- 15
n1. <- n11 + n12; n2. <- n21 + n22
# Risk Difference
p1<-(n11/n1.)
p2<-(n21/n2.)
rd <- p1 - p2
se_rd <- sqrt((p1 * (1 - p1) / n1.) + p2*((1 - p2) / n2.))
z_rd <- rd / se_rd
# Relative Risk
rr <- (n11/n1.) / (n21/n2.)
se_ln_rr <- sqrt((1/n11) - (1/n1.) + (1/n21) - (1/n2.))
z_rr <- log(rr) / se_ln_rr
# Odds Ratio
or <- (n11 * n22) / (n12 * n21)
se_ln_or <- sqrt((1/n11) + (1/n12) + (1/n21) + (1/n22))
z_or <- log(or) / se_ln_or
# Hasil
list(RD = rd, SE_RD = se_rd, Z_RD = z_rd, 
     RR = rr, SE_Ln_RR = se_ln_rr, Z_RR = z_rr,
     OR = or, SE_Ln_OR = se_ln_or, Z_OR = z_or)

## $RD
## [1] 0.03174603
## 
## $SE_RD
## [1] 0.1039056
## 
## $Z_RD
## [1] 0.3055277
## 
## $RR
## [1] 1.055556
## 
## $SE_Ln_RR
## [1] 0.1785255
## 
## $Z_RR
## [1] 0.3028544
## 
## $OR
## [1] 1.14
## 
## $SE_Ln_OR
## [1] 0.4277645
## 
## $Z_OR
## [1] 0.3063094

6.2.3 Uji Independensi

Uji independensi bertujuan untuk mengetahui apakah terdapat hubungan statistik antara dua variabel kategori dalam suatu populasi berdasarkan data dalam tabel kontingensi. Jika dua variabel bersifat independen, maka distribusi satu variabel tidak bergantung pada nilai dari variabel lainnya.

6.2.3.1 Uji Chi-Square

Uji Chi-Square merupakan metode paling umum untuk menguji independensi dalam tabel kontingensi. Uji ini digunakan untuk digunakan untuk menguji apakah ada hubungan antara dua variabel kategorikal

Rumus:

\[ \chi^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \]

Dimana:

\(O_{ij}\) adalah frekuensi yang diamati
\(E_{ij}\) adalah frekuensi yang diharapkan, dihitung dengan:

\[ E_{ij} = \frac{(n_{i.})(n_{.j})}{n} \]

Derajat kebebasan (df):

\[ (I - 1)(J - 1) \]

implementasi di r:

# Contoh Data
data <- matrix(c(12, 27, 18, 30), nrow = 2, byrow = TRUE)
dimnames(data) <- list("Jenis Kelamin" = c("Laki-laki", "Perempuan"), "Jenis Bacaan" = c("Fiksi Ilmiah", "Novel"))
print(data)

##              Jenis Bacaan
## Jenis Kelamin Fiksi Ilmiah Novel
##     Laki-laki           12    27
##     Perempuan           18    30

# Uji Chi-Square
chisq_test <- chisq.test(data)
print(chisq_test)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data
## X-squared = 0.18498, df = 1, p-value = 0.6671

6.2.3.2 Partisi Chi-Square

Partisi Chi-Square digunakan untuk mengidentifikasi kategori mana dalam tabel yang berkontribusi terhadap hasil uji signifikan. Langkah-langkahnya:

Hitung χ² keseluruhan.
Bagi tabel menjadi beberapa sub-tabel 2x2.
Hitung χ² masing-masing sub-tabel.
Evaluasi kontribusi tiap bagian terhadap total χ².

6.2.3.3 Uji Likelihood Ratio (G²) Uji rasio kemungkinan (Likelihood Ratio Test) menawarkan alternatif terhadap uji χ² dengan rumus:

\[ G^2 = 2 \sum\sum O_{ij} \ln \left(\frac{O_{ij}}{E_{ij}}\right) \] - \(O_{ij}\): frekuensi observasi - \(E_{ij}\): frekuensi yang diharapkan

Distribusi dari \(G^2\) mengikuti distribusi chi-square dengan df = \((I - 1)(J - 1)\).

# Data Observasi
data_matrix <- matrix(c(53, 40, 30, 57), nrow = 2, byrow = TRUE)
colnames(data_matrix) <- c("Fiksi Ilmiah", "Novel")
rownames(data_matrix) <- c("Politikus", "Penari")

# Hitung Frekuensi Ekspektasi
data_expected <- chisq.test(data_matrix)$expected
# Hitung Statistik G²
G2 <- 2 * sum(data_matrix * log(data_matrix / data_expected))
# Nilai kritis chi-square untuk df = 1 dan alpha = 0.05
critical_value <- qchisq(0.95, df = 1)
# Hasil
list(G2 = G2, Critical_Value = critical_value, Decision = ifelse(G2 > critical_value, "Reject H0", "Fail"))

## $G2
## [1] 9.252464
## 
## $Critical_Value
## [1] 3.841459
## 
## $Decision
## [1] "Reject H0"

6.2.3.4 Uji Exact Fisher Uji Fisher Exact merupakan metode statistik yang digunakan untuk menguji independensi antara dua variabel kategori dalam tabel kontingensi 2x2, terutama ketika ukuran sampel kecil dan asumsi uji chi-square tidak dapat terpenuhi.

Keunggulan: - Cocok untuk data dengan frekuensi kecil. - Tidak memerlukan asumsi distribusi normal.

6.2.3.4.1 Distribusi HIpergeometrik

Distribusi hipergeometrik menggambarkan probabilitas mengambil bola putih dalam pengambilan acak tanpa pengembalian dari kumpulan bola yang terdiri dari bola putih dan bola hitam.

Fungsi Probabilitas

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

( N ): total objek dalam populasi
( K ): jumlah objek dalam kategori tertentu (misalnya, bola putih atau sukses)
( n ): jumlah sampel yang diambil
( x ): jumlah objek dalam kategori tertentu yang diamati dalam sampel

6.3 Analisis Residual dalam Tabel Kontingensi

Analisis residual pada tabel kontingensi dilakukan untuk mengidentifikasi sel mana yang menyimpang jauh dari nilai harapan. Ketidaksesuaian ini dapat mengindikasikan pola tertentu dalam data atau ketidaktepatan model.

Residual dihitung sebagai selisih antara frekuensi observasi (\(n_{ij}\)) dan frekuensi harapan (\(E_{ij}\)), lalu dinormalisasi agar dapat dibandingkan antar sel.

6.3.1 Jenis Residual

Pearson Residual

Ini adalah bentuk residual yang paling umum digunakan dalam uji chi-square.

\[ r_{ij}^{(P)} = \frac{n_{ij} - E_{ij}}{\sqrt{E_{ij}}} \]

Pearson residual mencerminkan kontribusi relatif dari masing-masing sel terhadap nilai total chi-square.

Standardized (Adjusted) Residual

Residual ini telah disesuaikan untuk mempertimbangkan proporsi baris dan kolom, serta digunakan untuk interpretasi yang lebih presisi.

\[ r_{ij}^{(S)} = \frac{n_{ij} - E_{ij}}{\sqrt{E_{ij}(1 - \hat{p}_i)(1 - \hat{p}_j)}} \]

Dimana: - \(\hat{p}_i = \frac{n_{i.}}{n}\) adalah proporsi total baris ke-\(i\) - \(\hat{p}_j = \frac{n_{.j}}{n}\) adalah proporsi total kolom ke-\(j\) - \(n\) adalah total keseluruhan sampel

Standardized residual digunakan untuk mengukur seberapa jauh observasi menyimpang dari harapan dengan memperhitungkan variabilitas tambahan.

6.3.2 Deteksi Outlier dalam Analisis Data Kategori

Dalam analisis data kategori, outlier merujuk pada sel dalam tabel kontingensi yang memiliki residual sangat tinggi, baik dalam arah positif maupun negatif. Hal ini mengindikasikan bahwa kategori tersebut memiliki frekuensi pengamatan yang secara signifikan lebih besar atau lebih kecil dibandingkan dengan frekuensi yang diharapkan jika kategori-kategori tersebut saling bebas (independen).

# 1. Observasi data (frekuensi aktual)
obs <- matrix(c(63, 35, 39, 48), nrow = 2, byrow = TRUE)
colnames(obs) <- c("Laki-laki", "Perempuan")
rownames(obs) <- c("Politikus", "Penari")
obs

##           Laki-laki Perempuan
## Politikus        63        35
## Penari           39        48

# 2. Total baris, kolom, dan grand total
row_totals <- rowSums(obs)
col_totals <- colSums(obs)
n <- sum(obs)

# 3. Hitung frekuensi harapan (expected)
expected <- outer(row_totals, col_totals) / n
expected

##           Laki-laki Perempuan
## Politikus  54.03243  43.96757
## Penari     47.96757  39.03243

# 4. Hitung residual Pearson manual
residuals_pearson <- (obs - expected) / sqrt(expected)
residuals_pearson

##           Laki-laki Perempuan
## Politikus  1.219965 -1.352410
## Penari    -1.294794  1.435363

# Threshold umum
outlier_flag <- abs(residuals_pearson) > 2
# Tampilkan residual dan apakah termasuk outlier
data.frame(
 Cell = c("Politikus x Laki-laki", "Politikus x Perempuan", "Penari x Laki-laki", "Penari x Perempuan"),
 Residual = as.vector(residuals_pearson),
 Outlier = as.vector(outlier_flag))

hasil perhitungan menunjukkan tidak adanya residual yang menjadi outlier pada data

7 Tabel Kontingensi Tiga Arah

Tabel kontingensi tiga arah merupakan pengembangan dari tabel dua arah yang digunakan untuk menganalisis hubungan antara dua variabel kategorik dengan mempertimbangkan pengaruh dari variabel ketiga. Tabel ini memiliki struktur tiga dimensi (X × Y × Z) dan umum digunakan dalam studi observasional untuk mengontrol variabel perancu atau melakukan analisis berdasarkan strata tertentu.

7.1 Peluang dan Tabel

7.1.1 Parsial

Tabel parsial menunjukkan hubungan antara \(X\) dan \(Y\) pada level tetap dari \(Z\), misalnya \(Z = k\). Ini memungkinkan kita mengevaluasi hubungan yang dikendalikan terhadap \(Z\).

Formulasi:

\[ n_{ij|k} = n_{ijk} \]

Contoh struktur tabel parsial untuk \(Z = 1\):

	Y = 1	Y = 2	Total
X = 1	\(n_{111}\)	\(n_{121}\)	\(n_{1+1}\)
X = 2	\(n_{211}\)	\(n_{221}\)	\(n_{2+1}\)
Total	\(n_{+11}\)	\(n_{+21}\)	\(n_{++1}\)

Menurut Agresti (2007), tabel parsial diperlukan untuk mengevaluasi asosiasi bersyarat dan menghindari kesimpulan keliru dari tabel agregat.

7.1.2 Marginal

Tabel marginal pada kontingensi tiga arah dibentuk dengan menjumlahkan nilai dari variabel ketiga (\(Z\)), menghasilkan distribusi dua dimensi antara \(X\) dan \(Y\) tanpa mempertimbangkan pengaruh \(Z\).

Rumus:

\[ n_{ij+} = \sum_k n_{ijk} \]

Struktur Tabel Marginal \(X \times Y\)

	Y = 1	Y = 2	Total
X = 1	\(n_{11+}\)	\(n_{12+}\)	\(n_{1+}\)
X = 2	\(n_{21+}\)	\(n_{22+}\)	\(n_{2+}\)
Total	\(n_{+1+}\)	\(n_{+2+}\)	\(n_{++}\)

Marginalisasi ini digunakan untuk mengetahui hubungan secara keseluruhan antara \(X\) dan \(Y\) tanpa mempertimbangkan kondisi lain.

Namun, perlu diperhatikan bahwa teknik ini rentan terhadap efek Simpson’s Paradox, yaitu kondisi ketika tren atau asosiasi yang tampak pada data agregat berbalik arah ketika data dipisah ke dalam subkelompok.

Peluang marginal didapat dari penjumlahan peluang bersama terhadap satu variabel.

Misalnya:

\[ P(X = i) = \sum_{j,k} P(X = i, Y = j, Z = k) \]

Atau:

\[ P(X = i, Y = j) = \sum_{k} P(X = i, Y = j, Z = k) \]

Peluang marginal berguna untuk melihat distribusi umum variabel tunggal, tetapi tidak menunjukkan hubungan antar variabel.

7.1.3 Bersyarat

Peluang bersyarat menunjukkan probabilitas dari satu variabel, dengan syarat nilai tertentu dari dua variabel lainnya. Dalam konteks tabel kontingensi tiga arah, formulanya adalah:

\[ P(Y = j \mid X = i, Z = k) = \frac{n_{ijk}}{n_{i+k}} \]

Menurut Agresti (2007), peluang bersyarat merupakan fondasi dalam analisis hubungan serta estimasi parameter pada model log-linear dan regresi logistik bersyarat.

Tabel bersyarat menyajikan peluang bersyarat antara \(X\) dan \(Y\) pada setiap nilai tetap dari variabel \(Z\). Nilai-nilai dalam tabel ini berbentuk proporsi, bukan frekuensi absolut.

\[ P(Y = j \mid X = i, Z = k) = \frac{n_{ijk}}{n_{i+k}} \]

Struktur tabel bersyarat untuk \(Z = k\):

	Y = 1	Y = 2	Total
X = 1	\(P(1,1 \mid Z = k)\)	\(P(1,2 \mid Z = k)\)	\(P(1,+ \mid Z = k)\)
X = 2	\(P(2,1 \mid Z = k)\)	\(P(2,2 \mid Z = k)\)	\(P(2,+ \mid Z = k)\)
Total	\(P(+1 \mid Z = k)\)	\(P(+2 \mid Z = k)\)	\(P(+,+ \mid Z = k)\)

Tabel bersyarat ini memungkinkan kita untuk menganalisis pengaruh langsung antara dua variabel (X dan Y) setelah mengontrol efek variabel ketiga (Z). Oleh karena itu, ia menjadi dasar dari analisis multivariat dalam konteks data kategorik.

7.2 Ukuran Asosiasi

Ukuran asosiasi dalam tabel kontingensi digunakan untuk mengukur kekuatan hubungan antara dua variabel kategori. Tiga ukuran asosiasi yang paling umum digunakan adalah:

1. Risk Difference (RD)

\[ BP = P(Y \mid X_1, Z) - P(Y \mid X_2, Z) \]

2. Risiko Relatif (RR)

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

3. Odds Ratio (OR)

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

7.3 Conditional Independence

Independensi bersyarat terjadi ketika dua variabel kategori tidak saling bergantung secara statistik, dengan syarat bahwa nilai variabel ketiga sudah ditentukan atau dikendalikan. Dalam kasus ini, yang ingin dianalisis adalah apakah terdapat hubungan antara jenis kelamin dan jenis bacaan, setelah diklasifikasikan berdasarkan profesi.

Untuk mengevaluasi apakah independensi bersyarat terpenuhi, kita dapat:

Membandingkan ukuran asosiasi, seperti Odds Ratio (OR), pada setiap tabel parsial yang diperoleh dari memfiksasi variabel ketiga.
Bila nilai OR dari semua tabel parsial tersebut relatif konsisten, maka kita bisa mengasumsikan adanya independensi bersyarat.
Sebaliknya, jika terdapat perbedaan yang signifikan antar OR, maka hal ini menunjukkan bahwa variabel ketiga (misalnya profesi) memengaruhi hubungan antara dua variabel lainnya.

7.5 Marginal Y dan X

data <- array(
c(38, 20, 25, 15, 12, 18, 27, 30),
dim = c(2, 2, 2),
dimnames = list(
        JenisKelamin = c("Laki-laki", "Perempuan"),
        JenisBacaan = c("Fiksi Ilmiah", "Novel"),
        Profesi = c("Politikus", "Penari")
 )
)
data

## , , Profesi = Politikus
## 
##             JenisBacaan
## JenisKelamin Fiksi Ilmiah Novel
##    Laki-laki           38    25
##    Perempuan           20    15
## 
## , , Profesi = Penari
## 
##             JenisBacaan
## JenisKelamin Fiksi Ilmiah Novel
##    Laki-laki           12    27
##    Perempuan           18    30

Politikus

a <- data["Laki-laki", "Fiksi Ilmiah", "Politikus"]
b <- data["Perempuan", "Fiksi Ilmiah", "Politikus"]
c <- data["Laki-laki", "Novel", "Politikus"]
d <- data["Perempuan", "Novel", "Politikus"]
p1 <- a / (a + b)
p2 <- c / (c + d)
RD_Pol <- p1 - p2
RR_Pol <- p1 / p2
OR_Pol <- (a * d) / (b * c)
data_RD_RR_OR_Pol <- data.frame(RD = RD_Pol, RR = RR_Pol, OR = OR_Pol)
kable(data_RD_RR_OR_Pol, format = "latex", booktabs = TRUE)

a <- data["Laki-laki", "Fiksi Ilmiah", "Penari"]
b <- data["Perempuan", "Fiksi Ilmiah", "Penari"]
c <- data["Laki-laki", "Novel", "Penari"]
d <- data["Perempuan", "Novel", "Penari"]
p1 <- a / (a + b)
p2 <- c / (c + d)
RD_Pen <- p1 - p2
RR_Pen <- p1 / p2
OR_Pen <- (a * d) / (b * c)
data_RD_RR_OR_Pen <- data.frame(RD = RD_Pen, RR = RR_Pen, OR = OR_Pen)
kable(data_RD_RR_OR_Pen, format = "latex", booktabs = TRUE)

kesimpulan

profesi memengaruhi hubungan antara jenis kelamin dan jenis bacaan. bacaan fiksi ilmiah lebih banyak dibaca oleh politikus dibanding penari.

7.6 Inferensi Tabel kontingensi 3 arah

8 Generalized Linear Model (GLM)

Generalized Linear Model (GLM) merupakan pengembangan dari model regresi linear klasik yang dirancang untuk mengakomodasi berbagai tipe distribusi data. Model ini sangat bermanfaat ketika data yang dianalisis tidak memenuhi asumsi normalitas, seperti data kategori dan data count.

Komponen Utama GLM

GLM terdiri dari tiga komponen pokok:

Distribusi: Berasal dari keluarga distribusi eksponensial, seperti binomial, Poisson, dan normal.
Fungsi Link: Menghubungkan rata-rata (\(\mu\)) dengan kombinasi linear dari prediktor (\(\eta = X\beta\)).
Fungsi Varian: Menentukan bentuk varian, dinyatakan sebagai \(\text{Var}(Y) = \phi V(\mu)\).

8.1 Exponential Family

Distribusi yang digunakan dalam GLM harus termasuk dalam bentuk umum berikut:

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

Contoh Distribusi yang termasuk Exponential Family

Normal: \(\mu = E(Y), \quad \text{Var}(Y) = \sigma^2\)
Binomial: \(\mu = np, \quad \text{Var}(Y) = np(1-p)\)
Poisson: \(\mu = \lambda, \quad \text{Var}(Y) = \lambda\)
Gamma: \(\text{Var}(Y) = \mu^2\)

Model GLM mempermudah analisis data yang memiliki karakteristik berbeda-beda dengan tetap menggunakan pendekatan regresi yang fleksibel.

8.2 Model Regresi Logistik

Model regresi logistik digunakan ketika data respons bersifat biner (contohnya: sukses/gagal). Model ini menggunakan fungsi link logit:

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

Sehingga probabilitas kejadian dapat dituliskan sebagai:

\[ P(Y = 1 \mid X) = \frac{1}{1 + \exp(-X\beta)} \]

Komponen Model:

Distribusi: Binomial
Fungsi Link: Logit
Fungsi Varian: \(\mu (1 - \mu)\)

Contoh Implementasi di R

# Membuat data frame contoh
logit_data <- data.frame(
  usia = c(30, 33, 35, 30, 59, 35, 36, 39, 41, 34),
  pinjaman = c(0, 1, 0, 1, 0, 0, 0, 0, 0, 1)
)

# Membangun model regresi logistik
model_logit <- glm(pinjaman ~ usia, family = binomial(link = "logit"), data = logit_data)

# Menampilkan ringkasan model
summary(model_logit)

## 
## Call:
## glm(formula = pinjaman ~ usia, family = binomial(link = "logit"), 
##     data = logit_data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  13.7681    10.8155   1.273    0.203
## usia         -0.4250     0.3217  -1.321    0.187
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 12.2173  on 9  degrees of freedom
## Residual deviance:  8.8139  on 8  degrees of freedom
## AIC: 12.814
## 
## Number of Fisher Scoring iterations: 6

8.3 Model Regresi Poisson

Model ini digunakan untuk data yang berbentuk hitungan (count), seperti jumlah kejadian dalam interval waktu tertentu. Model ini menggunakan fungsi link log:

\[ \log(\mu) = X \beta \]

Sehingga:

\[ \mu = \exp(X \beta) \]

Komponen Model:

Distribusi: Poisson
Fungsi Link: Log
Fungsi Varian: \(\mu\)

Contoh Implementasi di R

# Membuat data frame baru
poisson_data <- data.frame(
  waktu = c(1, 2, 3, 4, 5, 6),
  jumlah = c(2, 3, 4, 6, 9, 12)
)

# Membangun model regresi Poisson
model_poisson <- glm(jumlah ~ waktu, family = poisson(link = "log"), data = poisson_data)

# Menampilkan ringkasan model
summary(model_poisson)

## 
## Call:
## glm(formula = jumlah ~ waktu, family = poisson(link = "log"), 
##     data = poisson_data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.3541     0.5159   0.686    0.492   
## waktu         0.3590     0.1092   3.289    0.001 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 12.136851  on 5  degrees of freedom
## Residual deviance:  0.038077  on 4  degrees of freedom
## AIC: 24.957
## 
## Number of Fisher Scoring iterations: 3

9 Inferensi GLM

9.1 Mencari Ekspektasi dan Varians dalam GLM

Dalam GLM, nilai ekspektasi (rataan) dan varians dari variabel respons ditentukan oleh fungsi link dan jenis distribusi yang digunakan. Secara umum, rumusnya sebagai berikut:

Ekspektasi:

\[ E(Y) = \mu = g^{-1}(X\beta) \]

Varians:

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

Dengan keterangan:

\(g^{-1}\) adalah fungsi invers dari fungsi link
\(\phi\) adalah parameter dispersi (contoh: \(\phi = 1\) untuk distribusi Binomial dan Poisson)
\(V(\mu)\) adalah fungsi varian, tergantung pada distribusi yang digunakan:
- Binomial: \(V(\mu) = \mu(1 - \mu)\)
- Poisson: \(V(\mu) = \mu\)
- Normal: \(V(\mu) = 1\)
- Gamma: \(V(\mu) = \mu^2\)

9.2 Metode Penaksiran Parameter

Penaksiran parameter dalam Generalized Linear Model (GLM) umumnya menggunakan pendekatan Maximum Likelihood Estimation (MLE). Mengingat struktur non-linear dari GLM, solusi analitik sulit diperoleh. Oleh karena itu, pendekatan numerik seperti Iteratively Reweighted Least Squares (IRLS) digunakan untuk memperoleh estimasi parameter.

9.2.1 Fungsi Likelihood

Fungsi likelihood dalam GLM berdasarkan distribusi keluarga eksponensial didefinisikan sebagai:

\[ L(\beta) = \prod_{i=1}^n f(y_i; \mu_i(\beta), \phi) \]

Log-likelihood dari fungsi tersebut adalah:

\[ \ell(\beta) = \sum_{i=1}^n \left[ \frac{y_i \theta_i - b(\theta_i)}{a(\phi)} + c(y_i, \phi) \right] \]

9.2.2 Iteratively Reweighted Least Squares (IRLS)

IRLS adalah algoritma yang digunakan untuk menyelesaikan estimasi parameter dalam GLM. Prosedurnya meliputi langkah-langkah berikut:

Tentukan tebakan awal untuk parameter \(\beta^{(0)}\)
Hitung nilai prediksi \(\mu^{(t)} = g^{-1}(X\beta^{(t)})\)
Bentuk matriks bobot \(W\) serta vektor respons kerja \(z\)
Perbarui parameter menggunakan regresi tertimbang:

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

Ulangi langkah-langkah tersebut hingga diperoleh konvergensi.

Contoh Implementasi di R

# Contoh data untuk distribusi Poisson
data_pois <- data.frame(
  x = c(1, 2, 3, 4, 5),
  y = c(1, 3, 4, 7, 10)
)

# Penerapan model GLM dengan link log (Poisson)
model_poisson <- glm(y ~ x, family = poisson(link = "log"), data = data_pois)
summary(model_poisson)

## 
## Call:
## glm(formula = y ~ x, family = poisson(link = "log"), data = data_pois)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -0.06153    0.66238  -0.093  0.92599   
## x            0.48287    0.16275   2.967  0.00301 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 10.50458  on 4  degrees of freedom
## Residual deviance:  0.37675  on 3  degrees of freedom
## AIC: 20.599
## 
## Number of Fisher Scoring iterations: 4

# Contoh data untuk distribusi Binomial (logit)
data_logit <- data.frame(
  x = c(0.2, 0.4, 0.6, 0.8, 1.0),
  y = c(0, 0, 1, 1, 1)
)

# Penerapan model GLM dengan link logit (Binomial)
model_log <- glm(y ~ x, family = binomial(link = "logit"), data = data_logit)

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

summary(model_log)

## 
## Call:
## glm(formula = y ~ x, family = binomial(link = "logit"), data = data_logit)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -115.0   218656.4  -0.001        1
## x              229.9   421323.8   0.001        1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 6.7301e+00  on 4  degrees of freedom
## Residual deviance: 4.1812e-10  on 3  degrees of freedom
## AIC: 4
## 
## Number of Fisher Scoring iterations: 24

Interpretasi Output

Estimate: Merupakan nilai koefisien estimasi \(\hat{\beta}\)
Std. Error: Kesalahan baku dari estimasi koefisien
z value: Statistik uji untuk hipotesis nol \(H_0: \beta = 0\)
Pr(>|z|): Nilai-p untuk menguji signifikansi koefisien

Metode ini penting untuk mengidentifikasi kontribusi masing-masing prediktor terhadap model serta menilai signifikansi hubungan antara variabel bebas dengan respons.

9.3 Diagnostik Model GLM

Setelah model GLM terbentuk, langkah penting berikutnya adalah melakukan diagnostik untuk memastikan bahwa model sesuai dan asumsi yang digunakan terpenuhi.

9.3.1 Deviance dan Null Deviance

Null Deviance: Deviance dari model yang hanya menggunakan intercept.
Residual Deviance: Deviance dari model setelah menambahkan prediktor.

Penurunan nilai deviance mengindikasikan bahwa model mengalami peningkatan kualitas.

9.3.2 AIC (Akaike Information Criterion)

AIC digunakan untuk membandingkan beberapa model: - Nilai AIC yang lebih kecil menunjukkan model yang lebih baik.

9.4 Diagnostik Residual

Residual Pearson:

\[ r_i = \frac{y_i - \hat{\mu}_i}{\sqrt{V(\hat{\mu}_i)}} \]

Deviance Residual: Mengukur kontribusi deviance dari masing-masing observasi.

9.4.1 Leverage dan Influential Points

Leverage mengukur seberapa besar pengaruh suatu observasi terhadap prediksi.
Cook’s Distance digunakan untuk mengidentifikasi observasi yang memiliki pengaruh besar terhadap keseluruhan model.

library(car)

## Loading required package: carData

# Contoh data baru untuk analisis diagnostik
set.seed(123)
data_diag <- data.frame(
  x = 1:10,
  y = rpois(10, lambda = exp(0.5 + 0.3 * (1:10)))
)

model_diag <- glm(y ~ x, family = poisson(link = "log"), data = data_diag)
cooks.distance(model_diag)

##          1          2          3          4          5          6          7 
## 0.07355730 0.03797404 0.03572944 0.10659530 0.22563295 0.24280889 0.01206427 
##          8          9         10 
## 0.01811583 0.31099573 0.17720238

9.4.2 Goodness-of-Fit

Uji goodness-of-fit dapat digunakan untuk mengevaluasi kesesuaian model terhadap data.

# Uji goodness-of-fit menggunakan chi-square
with(model_diag, 1 - pchisq(deviance, df.residual))

## [1] 0.2907491

9.4.3 Visualisasi Diagnostik

# Visualisasi diagnostik model
par(mfrow = c(2, 2))
plot(model_diag)

Kesimpulan

Diagnostik model sangat penting untuk memastikan bahwa model yang dibangun tidak hanya cocok secara statistik, namun juga valid dalam interpretasi. Analisis residual, deviance, leverage, dan Cook’s distance merupakan alat utama dalam mengidentifikasi potensi masalah dalam model GLM.

9.5 Detail Metode Estimasi dan Inferensi Regresi Logistik

Regresi logistik adalah bentuk khusus dari GLM yang digunakan ketika variabel respons bersifat biner (0/1). Estimasi parameter dilakukan melalui pendekatan Maximum Likelihood Estimation (MLE), dan inferensi didasarkan pada distribusi normal asimtotik dari koefisien estimasi.

Fungsi Likelihood untuk Regresi Logistik

\[ \ell(\beta) = \sum_{i=1}^n \left[ y_i \log(p_i) + (1 - y_i) \log(1 - p_i) \right], \quad \text{dengan } p_i = \frac{1}{1 + e^{-x_i^T\beta}} \]

Estimasi Parameter

Estimasi parameter \(\beta\) dilakukan dengan memaksimumkan fungsi log-likelihood di atas melalui metode numerik, umumnya menggunakan algoritma IRLS.

Inferensi: Confidence Interval dan Uji Hipotesis

Setelah diperoleh estimasi \(\hat{\beta}\), dilakukan inferensi sebagai berikut:

Asumsi Normal Asimtotik: \(\hat{\beta} \sim N(\beta, \text{Var}(\hat{\beta}))\)
Interval Kepercayaan:

\[ \hat{\beta}_j \pm Z_{1 - \alpha/2} \cdot SE(\hat{\beta}_j) \]

Uji Wald untuk hipotesis nol \(H_0: \beta_j = 0\)

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

Interpretasi Koefisien

Koefisien \(\beta_j\) merepresentasikan log odds ratio, sedangkan nilai \(e^{\beta_j}\) memberikan interpretasi sebagai odds ratio, yaitu perubahan odds untuk setiap peningkatan satu unit pada variabel prediktor.

Implementasi di R

# Contoh data baru untuk regresi logistik
set.seed(42)
data_logit2 <- data.frame(
  x = c(0.1, 0.3, 0.5, 0.7, 0.9, 1.1),
  y = c(0, 0, 1, 1, 1, 1)
)

model_logit2 <- glm(y ~ x, family = binomial(link = "logit"), data = data_logit2)

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

# Estimasi koefisien dan uji Wald
summary(model_logit2)

## 
## Call:
## glm(formula = y ~ x, family = binomial(link = "logit"), data = data_logit2)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -93.68  218291.05       0        1
## x              233.96  517295.76       0        1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 7.6382e+00  on 5  degrees of freedom
## Residual deviance: 2.7778e-10  on 4  degrees of freedom
## AIC: 4
## 
## Number of Fisher Scoring iterations: 25

# Interval kepercayaan berbasis profil likelihood
confint(model_logit2)

## Waiting for profiling to be done...

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

##                  2.5 %    97.5 %
## (Intercept) -150026.80  40794.08
## x            -60082.47 368039.50

# Odds ratio
exp(coef(model_logit2))

##   (Intercept)             x 
##  2.059192e-41 4.035105e+101

# Interval kepercayaan untuk odds ratio
exp(confint(model_logit2))

## Waiting for profiling to be done...

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

##             2.5 % 97.5 %
## (Intercept)     0    Inf
## x               0    Inf

Kesimpulan

Dalam regresi logistik, pendekatan MLE digunakan untuk memperoleh estimasi parameter yang efisien. Inferensi statistik dilakukan melalui uji Wald dan interval kepercayaan berdasarkan asumsi distribusi normal asimtotik. Koefisien dalam regresi logistik diinterpretasikan dalam konteks log odds dan odds ratio.

9.6 Detail Metode Estimasi dan Inferensi Regresi Poisson

Regresi Poisson digunakan untuk memodelkan data hitungan (count), seperti jumlah kejadian dalam suatu interval waktu. Model ini menggunakan distribusi Poisson dan fungsi link logaritmik.

Fungsi Likelihood untuk Regresi Poisson

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

Estimasi Parameter

Estimasi parameter \(\beta\) dilakukan dengan memaksimumkan fungsi log-likelihood. Prosedur ini umumnya dilakukan menggunakan metode numerik IRLS (Iteratively Reweighted Least Squares).

Confidence Interval dan Uji Hipotesis

Asumsi Normal Asimtotik: \(\hat{\beta} \sim N(\beta, \text{Var}(\hat{\beta}))\)
Interval Kepercayaan:

\[ \hat{\beta}_j \pm Z_{1 - \alpha/2} \cdot SE(\hat{\beta}_j) \]

Uji Wald untuk hipotesis nol \(H_0: \beta_j = 0\)

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

Interpretasi Koefisien

Koefisien \(\beta_j\) merepresentasikan perubahan pada log-count. Nilai \(e^{\beta_j}\) mengindikasikan rate ratio, yaitu perubahan expected count untuk setiap peningkatan satu unit pada prediktor.

Implementasi di R

# Contoh data baru untuk regresi Poisson
set.seed(7)
data_pois2 <- data.frame(
  x = 1:6,
  y = rpois(6, lambda = exp(0.8 + 0.4 * 1:6))
)

model_pois2 <- glm(y ~ x, family = poisson(link = "log"), data = data_pois2)

# Estimasi koefisien dan uji Wald
summary(model_pois2)

## 
## Call:
## glm(formula = y ~ x, family = poisson(link = "log"), data = data_pois2)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.93845    0.39311   2.387    0.017 *  
## x            0.34253    0.08375   4.090 4.32e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 27.0829  on 5  degrees of freedom
## Residual deviance:  8.5551  on 4  degrees of freedom
## AIC: 36.194
## 
## Number of Fisher Scoring iterations: 5

# Interval kepercayaan untuk koefisien
confint(model_pois2)

## Waiting for profiling to be done...

##                 2.5 %   97.5 %
## (Intercept) 0.1151943 1.661870
## x           0.1828311 0.512404

# Rate ratio
exp(coef(model_pois2))

## (Intercept)           x 
##    2.556024    1.408500

# Interval kepercayaan untuk rate ratio
exp(confint(model_pois2))

## Waiting for profiling to be done...

##                2.5 %   97.5 %
## (Intercept) 1.122091 5.269154
## x           1.200612 1.669299

10 Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio Regresi logistik

digunakan untuk memodelkan hubungan antara variabel dependen kategorik biner dengan satu atau lebih variabel independen, baik itu bersifat nominal, ordinal, maupun rasio. Bab ini akan menjelaskan cara menyimulasikan data dengan kombinasi prediktor tersebut, serta bagaimana mengevaluasi pengaruh masing-masing prediktor terhadap probabilitas kejadian suatu peristiwa.

10.1 Simulasi Data

untuk tujuan ilustrasi, kita akan mensimulasikan sebuah data fiktif tentang status penyakit jantung (penyakit = 0: tidak sakit, 1: sakit), yang dipengaruhi oleh:

jenis_kelamin (Nominal: “Laki-laki”, “Perempuan”)
tingkat_stres (Ordinal: “Rendah”, “Sedang”, “Tinggi”)
usia (Rasio, dalam tahun) Kesimpulan

Regresi Poisson merupakan metode yang efektif untuk menganalisis data hitungan. Estimasi parameter dilakukan dengan MLE melalui IRLS, sedangkan inferensi didasarkan pada distribusi normal asimtotik. Interpretasi koefisien umumnya difokuskan pada perubahan expected count atau rate ratio sebagai hasil dari perubahan nilai prediktor.

set.seed(123)

n <- 300

# Prediktor nominal
jenis_kelamin <- sample(c("Laki-laki", "Perempuan"), size = n, replace = TRUE)

# Prediktor ordinal
tingkat_stres <- sample(c("Rendah", "Sedang", "Tinggi"), size = n, replace = TRUE, prob = c(0.3, 0.5, 0.2))

# Prediktor rasio
usia <- round(rnorm(n, mean = 45, sd = 10), 1)

# Konversi ke nilai numerik untuk model
jk_num <- ifelse(jenis_kelamin == "Laki-laki", 1, 0)
stres_num <- as.numeric(factor(tingkat_stres, levels = c("Rendah", "Sedang", "Tinggi")))

# Fungsi logit: log(p / (1 - p)) = β0 + β1*jk + β2*stres + β3*usia
linpred <- -3 + 0.8 * jk_num + 0.6 * stres_num + 0.05 * usia
prob <- 1 / (1 + exp(-linpred))

# Variabel dependen biner
penyakit <- rbinom(n, size = 1, prob = prob)

# Buat data frame
data_simulasi <- data.frame(
  penyakit = factor(penyakit, labels = c("Tidak", "Ya")),
  jenis_kelamin,
  tingkat_stres,
  usia
)

# Tampilkan data simulasi
head(data_simulasi)

library(ggplot2)

ggplot(data_simulasi, aes(x = usia, fill = penyakit)) +
  geom_histogram(binwidth = 5, position = "dodge") +
  facet_wrap(~ jenis_kelamin) +
  labs(title = "Distribusi Usia Berdasarkan Jenis Kelamin dan Status Penyakit",
       x = "Usia", y = "Frekuensi") +
  theme_minimal()

10.2 Eksplorasi Data

Sebelum melakukan analisis lebih lanjut menggunakan regresi logistik, langkah penting yang perlu dilakukan adalah eksplorasi data. Tujuan dari eksplorasi ini adalah untuk memahami struktur data, distribusi variabel, serta potensi hubungan awal antara prediktor dengan variabel respons.

10.2.1 Statistik Deskriptif

# Paket yang digunakan 
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

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

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

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

# Ringkasan variabel numerik

summary(data_simulasi$usia)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   16.90   38.98   45.45   45.09   51.85   70.70

# Distribusi kategori

table(data_simulasi$jenis_kelamin)

## 
## Laki-laki Perempuan 
##       154       146

table(data_simulasi$tingkat_stres)

## 
## Rendah Sedang Tinggi 
##     80    155     65

table(data_simulasi$penyakit)

## 
## Tidak    Ya 
##    97   203

Interpretasi:

Rata-rata usia responden adalah sekitar nilai mean yang disimulasikan (45 tahun), dengan variasi sesuai standar deviasi.
Jenis kelamin dan tingkat stres tersebar cukup merata, meskipun distribusi stres mengikuti probabilitas yang telah ditentukan (lebih banyak yang ‘Sedang’).
Distribusi status penyakit (sakit atau tidak) menunjukkan bahwa sebagian besar responden tidak menderita penyakit jantung, namun proporsi antara kategori cukup seimbang untuk analisis klasifikasi.

10.2.2 Eksplorasi Hubungan Antara Variabel

Usia terhadap Status Penyakit

boxplot(usia ~ penyakit, data = data_simulasi,
        main = "Boxplot Usia berdasarkan Status Penyakit",
        xlab = "Status Penyakit", ylab = "Usia", col = c("lightblue", "salmon"))

Interpretasi:

Dari boxplot terlihat bahwa rata-rata usia responden yang menderita penyakit jantung cenderung lebih tinggi dibandingkan yang tidak sakit.
Sebaran usia juga lebih lebar pada kelompok penderita penyakit, menunjukkan bahwa usia bisa menjadi prediktor penting.

Jenis Kelamin terhadap Status Penyakit

table_gender_penyakit <- table(data_simulasi$jenis_kelamin, data_simulasi$penyakit)
prop.table(table_gender_penyakit, margin = 1)

##            
##                 Tidak        Ya
##   Laki-laki 0.2532468 0.7467532
##   Perempuan 0.3972603 0.6027397

Interpretasi:

Tabel proporsi menunjukkan bahwa laki-laki memiliki proporsi lebih tinggi untuk menderita penyakit dibandingkan perempuan.
Hal ini mengindikasikan bahwa jenis kelamin (nominal) mungkin memiliki pengaruh terhadap risiko penyakit jantung.

Tingkat Stres terhadap Status Penyakit

table_stres_penyakit <- table(data_simulasi$tingkat_stres, data_simulasi$penyakit)
prop.table(table_stres_penyakit, margin = 1)

##         
##              Tidak        Ya
##   Rendah 0.4125000 0.5875000
##   Sedang 0.3225806 0.6774194
##   Tinggi 0.2153846 0.7846154

Interpretasi:

Responden dengan tingkat stres ‘Tinggi’ cenderung memiliki proporsi penyakit lebih besar dibandingkan kelompok dengan stres ‘Rendah’ atau ‘Sedang’.
Pola ini mendukung asumsi bahwa tingkat stres (ordinal) berpotensi meningkatkan risiko penyakit jantung.

10.2.3 Visualisasi Data Kategorik

library(ggplot2)

ggplot(data_simulasi, aes(x = tingkat_stres, fill = penyakit)) +
  geom_bar(position = "fill") +
  facet_wrap(~ jenis_kelamin) +
  labs(title = "Proporsi Penyakit berdasarkan Tingkat Stres dan Jenis Kelamin",
       x = "Tingkat Stres", y = "Proporsi", fill = "Penyakit") +
  theme_minimal()

Interpretasi:

Visualisasi ini menampilkan bahwa, baik pada laki-laki maupun perempuan, proporsi penderita penyakit meningkat seiring peningkatan tingkat stres.
Ini mendukung eksplorasi sebelumnya dan menunjukkan adanya interaksi potensial antara jenis kelamin dan stres terhadap penyakit.

10.3 Perlakuan Variabel Ordinal

Variabel ordinal memiliki urutan tetapi tidak memiliki jarak yang tetap antar kategorinya. Dalam analisis regresi logistik, variabel ordinal dapat diperlakukan dengan dua pendekatan umum, yaitu:

Mengubahnya menjadi variabel dummy (perlakuan sebagai nominal)
Mengubahnya menjadi skala numerik peringkat (perlakuan sebagai rasio)

Eksperimen terhadap kedua pendekatan ini penting untuk melihat mana yang lebih baik merepresentasikan pengaruh variabel terhadap probabilitas kejadian.

10.3.1 Treat sebagai Nominal (Dummy)

Dalam pendekatan ini, setiap kategori dari variabel ordinal akan diubah menjadi variabel dummy terpisah. Ini berarti tidak ada asumsi tentang urutan antar kategori.

# Ubah tingkat stres menjadi variabel dummy
library(dplyr)
data_dummy <- data_simulasi %>%
  mutate(stres_sedang = ifelse(tingkat_stres == "Sedang", 1, 0),
         stres_tinggi = ifelse(tingkat_stres == "Tinggi", 1, 0))

# Regresi logistik dengan dummy variabel
model_dummy <- glm(penyakit ~ jenis_kelamin + stres_sedang + stres_tinggi + usia,
                   data = data_dummy, family = binomial)
summary(model_dummy)

## 
## Call:
## glm(formula = penyakit ~ jenis_kelamin + stres_sedang + stres_tinggi + 
##     usia, family = binomial, data = data_dummy)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            -2.07663    0.66687  -3.114  0.00185 ** 
## jenis_kelaminPerempuan -0.59140    0.26359  -2.244  0.02486 *  
## stres_sedang            0.31934    0.30053   1.063  0.28798    
## stres_tinggi            0.88675    0.39315   2.256  0.02410 *  
## usia                    0.06311    0.01415   4.461 8.16e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 377.61  on 299  degrees of freedom
## Residual deviance: 342.18  on 295  degrees of freedom
## AIC: 352.18
## 
## Number of Fisher Scoring iterations: 4

Interpretasi: - Kategori “Rendah” menjadi baseline/reference. - Koefisien untuk stres_sedang dan stres_tinggi menunjukkan log odds peningkatan risiko dibanding kategori “Rendah”. - Jika koefisien signifikan dan positif, artinya semakin tinggi stres, semakin tinggi pula probabilitas sakit.

10.3.2 Treat sebagai Rasio (Numeric Rank)

Jika kita yakin bahwa peningkatan tingkat stres merepresentasikan kenaikan tingkat risiko yang proporsional, maka kita dapat mengubah variabel ordinal menjadi angka peringkat.

# Konversi ordinal ke numeric
data_ordinal <- data_simulasi %>%
  mutate(stres_rank = as.numeric(factor(tingkat_stres, levels = c("Rendah", "Sedang", "Tinggi"))))

# Regresi logistik dengan stres sebagai numerik
model_ordinal <- glm(penyakit ~ jenis_kelamin + stres_rank + usia,
                     data = data_ordinal, family = binomial)
summary(model_ordinal)

## 
## Call:
## glm(formula = penyakit ~ jenis_kelamin + stres_rank + usia, family = binomial, 
##     data = data_ordinal)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            -2.53485    0.73387  -3.454 0.000552 ***
## jenis_kelaminPerempuan -0.59567    0.26342  -2.261 0.023743 *  
## stres_rank              0.42655    0.19069   2.237 0.025290 *  
## usia                    0.06275    0.01411   4.448 8.68e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 377.61  on 299  degrees of freedom
## Residual deviance: 342.39  on 296  degrees of freedom
## AIC: 350.39
## 
## Number of Fisher Scoring iterations: 3

Interpretasi: - Koefisien stres_rank menunjukkan pengaruh log-linear dari peningkatan tingkat stres. - Jika signifikan, ini berarti setiap kenaikan satu tingkat stres menaikkan log odds risiko penyakit. - Pendekatan ini lebih sederhana dan efisien jika asumsi linearitas antar level valid.

10.3.3 Perbandingan Pendekatan

# Bandingkan AIC kedua model
AIC(model_dummy, model_ordinal)

Interpretasi: - Nilai AIC yang lebih rendah menunjukkan model yang lebih baik dalam menjelaskan data. - Jika perbedaan kecil, model ordinal lebih disukai karena lebih sederhana. - Jika model dummy jauh lebih baik, maka hubungan antar level tidak linear dan lebih cocok ditangani secara terpisah.

Kesimpulan

Perlakuan terhadap variabel ordinal sangat mempengaruhi interpretasi model. Perlakuan dummy cocok jika tidak ada asumsi urutan logis antar level, sedangkan konversi ke numerik cocok jika peningkatan level diharapkan memberikan peningkatan efek secara bertahap dan proporsional. Pemilihan metode terbaik dapat dilihat dari evaluasi model seperti AIC atau goodness-of-fit secara keseluruhan.

11 Pemilihan Model Regresi Logistik dan Evaluasi

11.1 Membangun Model Regresi Logistik: Pendekatan Konfirmatori dan Eksploratori

Dalam penerapan regresi logistik, proses pemilihan model sangat penting karena akan menentukan sejauh mana model tersebut mampu memprediksi probabilitas dari suatu kejadian biner. Secara umum, terdapat dua pendekatan utama yang sering digunakan untuk membangun model, yaitu pendekatan Konfirmatori (Confirmatory) dan Eksploratori (Exploratory).

11.1.1 Pendekatan Konfirmator

Pendekatan konfirmatori digunakan ketika peneliti memiliki dasar teori atau dugaan awal (hipotesis) mengenai bagaimana hubungan antar variabel prediktor dengan variabel respon terbentuk. Pendekatan ini bersifat top-down, karena analisis dilakukan untuk menguji kebenaran struktur model yang telah dirancang sebelumnya berdasarkan referensi ilmiah atau studi sebelumnya.

Ciri-ciri pendekatan konfirmatori: - Model dibangun berdasarkan teori yang sudah ada atau didukung oleh penelitian sebelumnya. - Tujuan utamanya adalah untuk menguji apakah efek dari variabel prediktor signifikan terhadap respon, bukan untuk mencari model terbaik secara statistik. - Biasanya, peneliti langsung membangun model lengkap yang memuat semua variabel yang dianggap relevan, kemudian dilakukan pengujian terhadap pengaruh penambahan atau penghapusan komponen model (seperti interaksi). - Salah satu metode pengujian yang umum adalah Likelihood Ratio Test (LRT), yang membandingkan model lengkap dengan model yang telah dikurangi efek tertentu.

Contoh penerapan: Apabila teori menyebutkan bahwa variabel x1 dan x2 berpengaruh terhadap keputusan konsumen dalam membeli produk, maka model logistik akan langsung dibangun dengan memasukkan kedua variabel tersebut. Selanjutnya, dilakukan pengujian terhadap kontribusi x2 dalam model untuk melihat apakah variabel tersebut memiliki pengaruh signifikan.

11.1.2 Pendekatan Eksploratori

Pendekatan eksploratori digunakan dalam situasi ketika peneliti belum memiliki teori yang kuat atau ingin menelusuri kemungkinan pola hubungan antar variabel dari data yang tersedia. Tujuan utamanya adalah membangun model prediksi terbaik berdasarkan performa statistik.

Karakteristik pendekatan eksploratori: - Pembangunan model dilakukan secara bertahap, dengan fokus utama pada pencarian kombinasi variabel prediktor yang memberikan hasil terbaik. - Pemilihan variabel dalam model dilakukan dengan mempertimbangkan indikator statistik seperti Akaike Information Criterion (AIC), nilai deviance, atau log-likelihood. - Proses seleksi model dapat dilakukan dengan: - Forward Selection: Memulai dari model kosong, kemudian menambahkan satu per satu variabel jika terbukti signifikan. - Backward Elimination: Memulai dari model lengkap, kemudian secara bertahap menghapus variabel yang tidak signifikan. - Stepwise Selection: Kombinasi dari kedua metode di atas, di mana variabel dapat ditambahkan maupun dihapus secara iteratif berdasarkan evaluasi statistik.

Tujuan akhir: Menghasilkan model yang parsimonious, yaitu model yang cukup sederhana namun tetap memberikan performa prediksi yang tinggi.

11.1.3 Konsep AIC dan BIC

Dalam analisis regresi logistik, memilih model yang tepat sangat penting untuk memastikan bahwa model tidak terlalu sederhana (underfitting) atau terlalu kompleks (overfitting). Dua kriteria yang umum digunakan untuk pemilihan model adalah Akaike Information Criterion (AIC) dan Bayesian Information Criterion (BIC).

AIC mengukur kualitas relatif dari model statistik untuk kumpulan data tertentu. AIC menghargai model dengan goodness-of-fit yang baik dan penalti untuk kompleksitas model. Rumus AIC adalah:

\[ \text{AIC} = 2k - 2\ln(\hat{L}) \]

di mana:
- \(k\) adalah jumlah parameter dalam model,
- \(\hat{L}\) adalah nilai maksimum dari likelihood fungsi model.
BIC mirip dengan AIC tetapi memberikan penalti yang lebih besar untuk kompleksitas model, terutama saat ukuran sampel besar. Rumus BIC adalah:

\[ \text{BIC} = \ln(n)k - 2\ln(\hat{L}) \]

di mana:
- \(n\) adalah jumlah observasi dalam data.

Kedua kriteria ini membantu dalam memilih model yang seimbang antara fit dan kompleksitas. Model dengan nilai AIC atau BIC yang lebih rendah dianggap lebih baik dalam konteks relatif antar model.

Simulasi Data Untuk keperluan ilustrasi, kita akan membuat simulasi data baru dengan tiga variabel prediktor: - umur sebagai variabel numerik kontinu (diasumsikan berdistribusi normal), - status_merokok sebagai variabel biner (0 = tidak merokok, 1 = merokok), - aktivitas_fisik sebagai variabel numerik kontinu yang merepresentasikan intensitas aktivitas (misalnya skor skala 0–10).

Kemudian, kita bangun probabilitas kejadian (misalnya, mengalami_penyakit) berdasarkan fungsi logit dari ketiga variabel tersebut.

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

## 
## Attaching package: 'MASS'

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

library(caret)

## Warning: package 'caret' was built under R version 4.4.3

## Loading required package: lattice

library(pROC)

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

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

## 
## Attaching package: 'pROC'

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

library(DescTools)

## 
## Attaching package: 'DescTools'

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

## The following object is masked from 'package:car':
## 
##     Recode

set.seed(2025)

n <- 200
umur <- rnorm(n, mean = 40, sd = 12)
status_merokok <- rbinom(n, size = 1, prob = 0.4)
aktivitas_fisik <- rnorm(n, mean = 5, sd = 2)

# Fungsi linier prediktor
lin_pred <- -1 + 0.04 * umur + 1.1 * status_merokok - 0.6 * aktivitas_fisik

# Menghitung probabilitas
p <- 1 / (1 + exp(-lin_pred))

# Menentukan respons biner
mengalami_penyakit <- rbinom(n, size = 1, prob = p)

# Membuat data frame
simulasi_data <- data.frame(
  mengalami_penyakit = as.factor(mengalami_penyakit), umur, status_merokok, aktivitas_fisik
)

# Menampilkan data awal
head(simulasi_data)

Catatan: - Koefisien positif untuk status_merokok mengindikasikan bahwa perokok memiliki risiko lebih tinggi terhadap penyakit. - Koefisien negatif pada aktivitas_fisik mencerminkan bahwa aktivitas fisik tinggi dapat menurunkan risiko penyakit. - umur memberikan pengaruh positif moderat terhadap kemungkinan mengalami penyakit.

Data ini dapat digunakan untuk melanjutkan analisis regresi logistik dan membandingkan model sesuai pendekatan konfirmatori atau eksploratori.

Pemilihan Model

model_full <- glm(mengalami_penyakit ~ umur + status_merokok + aktivitas_fisik, data = simulasi_data, family = binomial)
summary(model_full)

## 
## Call:
## glm(formula = mengalami_penyakit ~ umur + status_merokok + aktivitas_fisik, 
##     family = binomial, data = simulasi_data)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -1.40439    0.89716  -1.565  0.11750    
## umur             0.04861    0.01792   2.713  0.00667 ** 
## status_merokok   0.27670    0.40445   0.684  0.49389    
## aktivitas_fisik -0.49032    0.11543  -4.248 2.16e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 194.49  on 199  degrees of freedom
## Residual deviance: 164.04  on 196  degrees of freedom
## AIC: 172.04
## 
## Number of Fisher Scoring iterations: 5

11.2 Metode Stepwise:Forward, Backward, dan Kedua Arah

null_model <- glm(mengalami_penyakit ~ 1, data = simulasi_data, family = binomial)
step_forward <- step(null_model, direction = "forward", scope = formula(model_full), trace = FALSE)
step_backward <- step(model_full, direction = "backward", trace = FALSE)
step_both <- step(null_model, direction = "both", scope = formula(model_full), trace = FALSE)
AIC(model_full, step_forward, step_backward, step_both)

11.3 Evaluasi Model: ROC dan AUC

pred_prob <- predict(step_both, type = "response")
roc_obj <- roc(simulasi_data$mengalami_penyakit, pred_prob)

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

## Setting direction: controls < cases

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

auc(roc_obj)

## Area under the curve: 0.755

11.4 Pseudo R-Squared

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

##   CoxSnell Nagelkerke   McFadden 
##  0.1392308  0.2239002  0.1541771

11.5 Tabel Klasifikasi

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

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 161  32
##          1   1   6
##                                           
##                Accuracy : 0.835           
##                  95% CI : (0.7762, 0.8836)
##     No Information Rate : 0.81            
##     P-Value [Acc > NIR] : 0.2105          
##                                           
##                   Kappa : 0.2206          
##                                           
##  Mcnemar's Test P-Value : 1.767e-07       
##                                           
##             Sensitivity : 0.1579          
##             Specificity : 0.9938          
##          Pos Pred Value : 0.8571          
##          Neg Pred Value : 0.8342          
##              Prevalence : 0.1900          
##          Detection Rate : 0.0300          
##    Detection Prevalence : 0.0350          
##       Balanced Accuracy : 0.5759          
##                                           
##        'Positive' Class : 1               
##

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

## Sensitivity Specificity 
##   0.1578947   0.9938272

11.6 Metode Perbandingan Model dalam Regresi Logistik

Dalam regresi logistik, kita sering membandingkan beberapa model untuk memilih model terbaik yang merepresentasikan data secara optimal. Salah satu metode perbandingan model yang umum digunakan adalah Likelihood-Ratio Test (LRT).

Rumus Likelihood-Ratio Test:

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

\(L_0\): Likelihood dari model sederhana (misalnya model null).
\(L_1\): Likelihood dari model kompleks (misalnya model dengan lebih banyak prediktor).
\(G^2\) mengikuti distribusi chi-square dengan derajat bebas sebesar selisih jumlah parameter antara kedua model.

# Simulasi data
set.seed(123)
n <- 200
x1 <- rnorm(n, mean = 0, sd = 1)
x2 <- rbinom(n, 1, 0.5)
linpred <- -0.5 + 1 * x1 + 0.8 * x2
p <- 1 / (1 + exp(-linpred))
y <- rbinom(n, 1, p)
data <- data.frame(y, x1, x2)

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

# Model penuh
model_full <- glm(y ~ x1 + x2, family = binomial, data = data)

# Likelihood-Ratio Test
anova(model_null, model_full, test = "Chisq")

Interpretasi: Hasil anova menunjukkan apakah model penuh secara signifikan lebih baik dari model null. Jika nilai p < 0.05, maka model penuh memberikan peningkatan signifikan dalam penjelasan variabilitas respons.

11.7 Likelihood-Ratio Test

Uji LRT adalah teknik penting untuk memutuskan apakah penambahan variabel ke dalam model secara signifikan meningkatkan kualitas model.

Lanjutkan perhitungan \(G^2\) secara eksplisit:

G2 <- -2 * (logLik(model_null) - logLik(model_full))
G2

## 'log Lik.' 57.02612 (df=1)

Interpretasi: Nilai \(G^2\) menunjukkan peningkatan log-likelihood saat menambahkan prediktor. Semakin besar \(G^2\), semakin baik model menjelaskan data.

11.8 Prinsip Parsimony

Prinsip parsimony menyatakan bahwa dari sekian banyak model yang mungkin, kita memilih model yang paling sederhana yang tetap menjelaskan data dengan baik.

Hindari model yang terlalu kompleks (overfitting).
Prioritaskan model dengan prediktor yang benar-benar berkontribusi signifikan.

summary(model_full)

## 
## Call:
## glm(formula = y ~ x1 + x2, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.2333     0.2337  -0.998    0.318    
## x1            1.4250     0.2317   6.150 7.77e-10 ***
## x2            0.2867     0.3291   0.871    0.384    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 276.76  on 199  degrees of freedom
## Residual deviance: 219.73  on 197  degrees of freedom
## AIC: 225.73
## 
## Number of Fisher Scoring iterations: 4

Interpretasi: Evaluasi nilai p untuk setiap koefisien. Jika p > 0.05, pertimbangkan untuk menghapus prediktor tersebut demi penyederhanaan model.

11.9 Evaluasi Tabel Klasifikasi dan Akurasi Model

Setelah memilih model terbaik, kita perlu mengevaluasi kinerjanya melalui:

Tabel klasifikasi
Akurasi
Sensitivitas & Spesifisitas

# Probabilitas prediksi
prob_pred <- predict(model_full, type = "response")
class_pred <- ifelse(prob_pred > 0.5, 1, 0)

table(Actual = data$y, Predicted = class_pred)

##       Predicted
## Actual  0  1
##      0 79 26
##      1 35 60

# Akurasi
mean(class_pred == data$y)

## [1] 0.695

Interpretasi: - Tabel klasifikasi menunjukkan distribusi prediksi benar/salah. - Akurasi menunjukkan proporsi prediksi yang benar.

11.9.1 Sensitivitas dan Spesifisitas

Sensitivitas (True Positive Rate): Kemampuan model untuk mendeteksi kelas positif dengan benar.
Spesifisitas (True Negative Rate): Kemampuan model untuk mendeteksi kelas negatif dengan benar.

TP <- sum(data$y == 1 & class_pred == 1)
TN <- sum(data$y == 0 & class_pred == 0)
FP <- sum(data$y == 0 & class_pred == 1)
FN <- sum(data$y == 1 & class_pred == 0)

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

list(Sensitivitas = sensitivitas, Spesifisitas = spesifisitas)

## $Sensitivitas
## [1] 0.6315789
## 
## $Spesifisitas
## [1] 0.752381

Interpretasi: - Sensitivitas tinggi berarti model baik dalam menangkap kelas positif. - Spesifisitas tinggi berarti model baik dalam menghindari false positive.

Kesimpulan: Pada bagian ini, kita telah: - Menggunakan LRT untuk membandingkan model - Menerapkan prinsip parsimony untuk memilih model terbaik - Mengevaluasi kinerja model melalui akurasi, sensitivitas, dan spesifisitas

Langkah-langkah ini penting dalam memilih model logistik yang tidak hanya fit terhadap data, tetapi juga generalizable.

11.10 Penjelasan Kurva ROC (Receiver Operating Characteristic)

Kurva ROC digunakan untuk mengevaluasi kemampuan model dalam membedakan antara dua kelas.

Sumbu x: False Positive Rate (1 - spesifisitas)
Sumbu y: True Positive Rate (sensitivitas)
AUC (Area Under Curve): Luas di bawah kurva ROC, menggambarkan akurasi model secara keseluruhan.

library(pROC)
roc_obj <- roc(data$y, prob_pred)

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

## Setting direction: controls < cases

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

auc(roc_obj)

## Area under the curve: 0.7843

Interpretasi: - Semakin dekat AUC ke 1, semakin baik performa model. - AUC = 0.5 menunjukkan model tidak baik

11.11 Precision-Recall Curve (PR Curve)

Kurva PR adalah alternatif ROC, khususnya saat data tidak seimbang (banyak kelas 0 atau 1).

Precision: Proporsi prediksi positif yang benar-benar positif
Recall: Sama dengan sensitivitas

library(PRROC)

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

## Loading required package: rlang

pr <- pr.curve(scores.class0 = prob_pred[data$y == 1],
               scores.class1 = prob_pred[data$y == 0],
               curve = TRUE)
plot(pr)

Interpretasi: - PR curve lebih informatif daripada ROC jika proporsi kelas positif sangat kecil. - Area yang lebih besar menunjukkan kinerja model yang lebih baik.

11.12 Pseudo R-squared pada Regresi Logistik

Pseudo R-squared digunakan untuk mengukur seberapa baik model logistik menjelaskan variasi data, mirip dengan \(R^2\) pada regresi linear.

Beberapa jenis pseudo R-squared:

McFadden’s R²:

\[ R^2_{McF} = 1 - \frac{\log L_{model}}{\log L_{null}} \]

Cox & Snell R²:

\[ R^2_{CS} = 1 - \left(\frac{L_{null}}{L_{model}}\right)^{\frac{2}{n}} \]

Nagelkerke R²:

Merupakan normalisasi dari Cox & Snell agar maksimum 1.

library(pscl)

## Warning: package 'pscl' was built under R version 4.4.3

## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.

pR2(model_full)

## fitting null model for pseudo-r2

##         llh     llhNull          G2    McFadden        r2ML        r2CU 
## -109.866271 -138.379332   57.026122    0.206050    0.248084    0.331055

Interpretasi: - Nilai yang lebih tinggi menunjukkan model lebih baik. - Tidak bisa dibandingkan langsung dengan \(R^2\) pada regresi linear, tapi bermanfaat untuk membandingkan antar model logistik.

Kesimpulan: Bagian ini memperluas evaluasi model logistik melalui: - Visualisasi ROC & PR curve - Perhitungan AUC dan precision-recall - Evaluasi menggunakan pseudo R-squared

Evaluasi ini memberi wawasan menyeluruh tentang seberapa baik model klasifikasi bekerja dan seberapa akurat prediksinya pada data nyata.

12 Regresi Logistik Multinomial

12.1 Distribusi Multinomial

Distribusi multinomial adalah perluasan dari distribusi binomial, yang digunakan ketika setiap percobaan tidak hanya menghasilkan dua hasil (sukses/gagal), tetapi bisa menghasilkan salah satu dari \(K\) kategori. Misalnya, melempar sebuah dadu dengan \(K\) sisi, setiap sisi memiliki probabilitas muncul tertentu.

Notasi dan Definisi

Misalkan terdapat \(K\) kategori, masing-masing dengan probabilitas:

\[ P = [P_1, P_2, \dots, P_K]^T \]

dengan syarat:

\[ \sum_{k=1}^{K} P_k = 1 \]

Setelah dilakukan \(n\) percobaan, misalkan hasilnya:

\[ x = [x_1, x_2, \dots, x_K]^T \]

dengan jumlah total:

\[ \sum_{k=1}^{K} x_k = n \]

Maka vektor acak \(x\) mengikuti distribusi multinomial:

\[ x \sim \text{Mult}(x \mid n, P) \]

Fungsi Probabilitas

Fungsi probabilitas dari distribusi multinomial adalah:

\[ P(x) = \text{Mult}(x \mid n, P) = \frac{n!}{x_1! x_2! \cdots x_K!} \prod_{k=1}^{K} P_k^{x_k} \]

Rumus ini menyatakan probabilitas mendapatkan kombinasi hasil \(x_1, x_2, \dots, x_K\) dari \(n\) percobaan dengan probabilitas kategori masing-masing \(P_k\).

Sifat Statistik

Nilai harapan:

\[ \mathbb{E}[x] = nP = [nP_1, nP_2, \dots, nP_K]^T \]

Varians tiap komponen:

\[ \text{Var}(x_k) = n P_k (1 - P_k), \quad k = 1, 2, \dots, K \]

Kovarians antar komponen (untuk \(i \ne j\)):

\[ \text{Cov}(x_i, x_j) = -n P_i P_j \]

Kasus Khusus: Distribusi Kategorikal

Jika hanya dilakukan satu percobaan (\(n = 1\)), maka distribusi multinomial menjadi distribusi kategorikal, yaitu:

\[ x \sim \text{Categorical}(P) \]

Distribusi ini adalah bentuk umum dari distribusi Bernoulli, ketika terdapat lebih dari dua kategori.

Contoh Kasus

Dilakukan sebuah survey kepada 20 fans F1 mereka disuruh memilih antara 4 tim top di F1 : Ferrari, McLaren, Red Bull, dan Mercedes Hasil survei:

Ferrari: 9 orang
McLaren: 4 orang
Red Bull: 4 orang
Mercedes: 3 orang Probabilitas teoretik preferensi:
\(p_F = 0{,}45\)
\(p_Mc = 0{,}2\)
\(p_R = 0{,}2\)
\(p_M = 0{,}15\)

Pertanyaannya:

Berapa peluang bahwa dalam 20 orang akan ada 9 yang memilih Ferrari, 4 memilih McLaren, 4 memilih Red Bull dan 3 memilih Mercedes?

Rumus Distribusi Multinomial

Distribusi peluang multinomial:

\[ P(X_1 = x_1, \dots, X_K = x_K) = \frac{n!}{x_1! x_2! \cdots x_K!} p_1^{x_1} p_2^{x_2} \cdots p_K^{x_K} \]

dengan:

\(n = 20\), \(x_1 = 9\), \(x_2 = 4\), \(x_3 = 4\), \(x_4 = 3\)
\(p_1 = 0{,}45\), \(p_2 = 0{,}2\), \(p_3 = 0{,}2\), \(p_4 = 0{,}15\)

perhitungan manual di R

# Jumlah percobaan
n <- 20

#Hasil percobaan
x <- c(9,4,4,3)
# Probabilitas masing-masing kategori (jumlah harus = 1)
p <- c(0.45, 0.2, 0.2, 0.15)

# Hitung komponen-koefisien
faktorial_total <- factorial(n)
faktorial_x <- prod(factorial(x))
koefisien <- faktorial_total / faktorial_x
# Hitung peluang
peluang <- koefisien * prod(p^x)
peluang

## [1] 0.01268277

Probabilitas hal tersebut dapat terjadi ialah 0.01268277. Distribusi multinomial digunakan untuk menghitung peluang dalam percobaan dengan beberapa kategori hasil. Rumus dasarnya merupakan generalisasi dari binomial untuk lebih dari dua kategori. Distribusi multinomial memodelkan percobaan dengan hasil kategori jamak. Ini merupakan generalisasi dari distribusi binomial dan mencakup distribusi Bernoulli dan kategorikal sebagai kasus khusus.

12.2 Regresi Logistik Multinomial

Model regresi logistik multinomial merupakan perluasan dari regresi logistik biner yang digunakan ketika variabel respon bersifat kategorik dengan lebih dari dua kategori (kategori nominal), dan kita ingin memodelkan pengaruh satu atau lebih variabel prediktor terhadap peluang masing-masing kategori respon.

Model ini berguna ketika kita ingin mengetahui faktor-faktor yang memengaruhi kemungkinan suatu pengamatan termasuk ke dalam salah satu dari \(K\) kategori yang tidak berurutan. Misalnya, dalam konteks preferensi produk, pilihan jurusan, atau hasil diagnosis penyakit dengan beberapa kemungkinan diagnosis yang berbeda.

Untuk \(K\) kategori, model membandingkan setiap kategori \(j = 1, 2, \dots, K - 1\) terhadap kategori referensi (biasanya kategori ke-\(K\)), dengan bentuk model sebagai berikut:

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

untuk \(j = 1, 2, \dots, K - 1\).

Setiap persamaan logit merepresentasikan log odds dari kategori \(j\) terhadap kategori referensi \(K\) sebagai fungsi linear dari prediktor \(x_1, x_2, \dots, x_p\). Dengan demikian, model ini memiliki \((K-1)(p+1)\) parameter, mencakup satu intercept dan \(p\) koefisien regresi untuk setiap logit.

Estimasi parameter biasanya dilakukan dengan metode maximum likelihood estimation (MLE), dan interpretasi koefisien mengikuti prinsip yang serupa dengan regresi logistik biner, namun dengan mempertimbangkan bahwa interpretasi berlaku terhadap kategori referensi.

12.3 Contoh Kasus

seorang peneliti ingin mengetahui pengaruh jenis kelamin terhadap pemilihan tempat makan favorit diantara 3 tempat yaitu : Burger King, McDonald, dan KFC. Ia melakukan survey terhadap 200 orang

# Gunakan kode berikut secara lengkap dan terstruktur
library(nnet)

set.seed(123)
n <- 200

# Gender vector (200 obs)
Gender <- sample(c("Male", "Female"), size = n, replace = TRUE)

# Food preference conditional on gender
food <- sapply(Gender, function(gen) {
  if (gen == "Male") {
    sample(c("Burger King", "McDonald", "KFC"), size = 1, prob = c(0.3, 0.5, 0.2))
  } else {
    sample(c("Burger King", "McDonald", "KFC"), size = 1, prob = c(0.1, 0.7, 0.2))
  }
})

# Gabungkan jadi data frame
food_data <- data.frame(
  food = factor(food),
  Gender = factor(Gender)
)

# Pastikan level referensi diatur (opsional)
food_data$food <- relevel(food_data$food, ref = "Burger King")
head(food_data)

# Fit model multinomial logistic regression
model_mlr <- multinom(food ~ Gender, data = food_data, trace = FALSE)

# Ringkasan hasil
summary(model_mlr)

## Call:
## multinom(formula = food ~ Gender, data = food_data, trace = FALSE)
## 
## Coefficients:
##          (Intercept) GenderMale
## KFC        0.5978643  -1.259268
## McDonald   1.7917398  -1.200394
## 
## Std. Errors:
##          (Intercept) GenderMale
## KFC        0.3753753  0.4854524
## McDonald   0.3256687  0.3951899
## 
## Residual Deviance: 364.1613 
## AIC: 372.1613

# Z dan p-value
library(knitr)
z_values <- summary(model_mlr)$coefficients / summary(model_mlr)$standard.errors
p_values <- 2 * (1 - pnorm(abs(z_values)))

coef_table <- cbind(summary(model_mlr)$coefficients,
                    "z value" = round(z_values, 2),
                    "p value" = round(p_values, 4))

knitr::kable(coef_table, caption = "Koefisien Model Multinomial Logistic Regression")

Koefisien Model Multinomial Logistic Regression
	(Intercept)	GenderMale	(Intercept)	GenderMale	(Intercept)	GenderMale
KFC	0.5978643	-1.259268	1.59	-2.59	0.1112	0.0095
McDonald	1.7917398	-1.200394	5.50	-3.04	0.0000	0.0024

Interpretasi : - Gender memiliki pengaruh terhadap preferensi tempat makan - peluang Male memilih KFC dibanding Burger King adalah 72% lebih rendah daripada Female. - peluang Male memilih McDonald dibanding Burger King adalah 70% lebih rendah daripada Female.

prediksi

# Buat data baru berdasarkan level Gender
new_data <- expand.grid(
  Gender = levels(food_data$Gender)
)

# Prediksi probabilitas (hasil berupa matriks)
pred_probs <- predict(model_mlr, newdata = new_data, type = "probs")

# Gabungkan dengan data baru
final_data <- cbind(new_data, pred_probs)
print(final_data)

##   Gender Burger King       KFC  McDonald
## 1 Female   0.1134029 0.2061928 0.6804043
## 2   Male   0.3009741 0.1553406 0.5436852

# Model null: hanya intercept
model_null <- multinom(food ~ 1, data = food_data, trace = FALSE)

# Log-likelihood
ll_null <- logLik(model_null)
ll_full <- logLik(model_mlr)

# G2 (deviance)
G2 <- -2 * (as.numeric(ll_null) - as.numeric(ll_full))

# Derajat kebebasan
df <- attr(ll_full, "df") - attr(ll_null, "df")

# p-value
pval <- 1 - pchisq(G2, df)

# McFadden's pseudo R2
pseudo_r2 <- 1 - (as.numeric(ll_full) / as.numeric(ll_null))

# Tampilkan hasil
cat("Nilai G² (Deviance):", round(G2, 4), "\n")

## Nilai G² (Deviance): 11.007

cat("Derajat Kebebasan:", df, "\n")

## Derajat Kebebasan: 2

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

## p-value: 0.0041

cat("McFadden's Pseudo R²:", round(pseudo_r2, 4), "\n")

## McFadden's Pseudo R²: 0.0293

13 Regresi Logistik Ordinal

Regresi logistik ordinal digunakan ketika variabel respon \(Y\) bersifat ordinal (memiliki urutan), misalnya tingkat kepuasan: Rendah, Sedang, Tinggi.

Model ini berbeda dengan:

Regresi logistik biner: hanya 2 kategori
Regresi logistik multinomial: kategori > 2 tetapi tidak berurutan

13.1 Konsep Cumulative Logit Model

Model yang digunakan adalah Cumulative Logit Model dengan asumsi proportional odds:

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

\(\alpha_j\): intercept khusus untuk kategori ke-\(j\)
\(\beta\): koefisien regresi (sama untuk semua kategori kumulatif)

Untuk \(c\) kategori, terdapat \((c - 1)\) model logit kumulatif.

13.2 Interpretasi Koefisien

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

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

Odds ratio:

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

13.3 Contoh

misalkan ada data tingkat kepuasan pasien (1: Buruk, 2: Sedang, 3: Baik) terhadap waktu tunggu sebelum dilayani dokter (dalam menit).

set.seed(2025)
n <- 200

# Waktu tunggu antara 5–60 menit
wait_time <- round(runif(n, 5, 60))

# Skor acak tergantung waktu tunggu (semakin lama, semakin tidak puas)
# Asumsikan threshold 25 dan 40 menit sebagai batas cut-off
latent_score <- 10 - 0.1 * wait_time + rnorm(n)

# Konversi ke ordinal satisfaction
satisfaction <- cut(latent_score,
                    breaks = c(-Inf, 6, 8, Inf),
                    labels = c("Buruk", "Sedang", "Baik"),
                    ordered_result = TRUE)

# Buat data frame
df <- data.frame(satisfaction, wait_time)

# Lihat data awal
head(df)

estimasi model

model_ordinal <- polr(satisfaction ~ wait_time, data = df, Hess = TRUE)

summary(model_ordinal)

## Call:
## polr(formula = satisfaction ~ wait_time, data = df, Hess = TRUE)
## 
## Coefficients:
##             Value Std. Error t value
## wait_time -0.2045    0.02133  -9.586
## 
## Intercepts:
##              Value   Std. Error t value
## Buruk|Sedang -8.3552  0.9033    -9.2494
## Sedang|Baik  -4.0999  0.5111    -8.0214
## 
## Residual Deviance: 212.1853 
## AIC: 218.1853

uji signifikansi dan P-value

# Hitung z dan p-value
coefs <- coef(summary(model_ordinal))
p_values <- pnorm(abs(coefs[, "t value"]), lower.tail = FALSE) * 2
cbind(coefs, "p value" = round(p_values, 4)) %>%
  knitr::kable(caption = "Koefisien dan p-value dari Model Regresi Logistik Ordinal")

Koefisien dan p-value dari Model Regresi Logistik Ordinal
	Value	Std. Error	t value
wait_time	-0.2044537	0.0213290	-9.585694
Buruk\|Sedang	-8.3552105	0.9033290	-9.249355
Sedang\|Baik	-4.0999423	0.5111286	-8.021352

Interpretasi

Koefisien wait_time = -0.2045 - Nilai koefisien negatif menunjukkan bahwa semakin lama waktu tunggu, semakin kecil kemungkinan pasien memberikan nilai kepuasan yang lebih tinggi. - Artinya, waktu tunggu yang lebih panjang menurunkan kepuasan pasien secara signifikan.

Odds Ratio: \[ \text{OR} = e^{-0.2045} \approx 0.815 \] Interpretasi: setiap kenaikan 1 menit waktu tunggu, odds untuk berada di kategori kepuasan yang lebih tinggi menurun sekitar 18.5%.

Cut Points (Intercepts) - Buruk|Sedang = -8.3552 → Batas antara kategori Buruk dan (Sedang atau Baik). - Sedang|Baik = -4.0999 → Batas antara kategori (Buruk atau Sedang) dan Baik.

Cut point adalah bagian dari model logit kumulatif dan digunakan untuk mengestimasi peluang, namun tidak diinterpretasikan secara langsung seperti koefisien prediktor.

Signifikansi Statistik

Koefisien	Estimate	Std. Error	t value
`wait_time`	-0.2045	0.0213	-9.59
`Buruk\|Sedang`	-8.3552	0.9033	-9.25
`Sedang\|Baik`	-4.0999	0.5111	-8.02

Semua nilai p sangat kecil (p < 0.001), menunjukkan bahwa seluruh parameter model signifikan secara statistik.
Efek waktu tunggu pada kepuasan sangat nyata dan tidak terjadi secara kebetulan.

Kesimpulan

Waktu tunggu pasien berpengaruh signifikan terhadap tingkat kepuasan.
Pasien yang menunggu lebih lama cenderung melaporkan kepuasan yang lebih rendah.
Model ini cocok untuk data ordinal karena mempertimbangkan urutan antar kategori.
Regresi logistik ordinal dengan pendekatan cumulative logit memberikan pemodelan yang efektif untuk hubungan antara variabel ordinal dan prediktor kontinu seperti wait_time.

14 Log Linear Model

Model log-linier adalah bentuk khusus dari Generalized Linear Model (GLM) yang digunakan pada frekuensi sel dalam tabel kontingensi, dan mengasumsikan distribusi Poisson.

Tidak seperti regresi logistik, model log-linier tidak menetapkan variabel mana yang dependen dan independen, karena seluruh variabel diperlakukan secara simetris. Model ini lebih cocok digunakan ketika tujuan analisis adalah untuk memahami struktur asosiasi atau independensi antar variabel, bukan untuk prediksi.

Struktur model log-linier ditentukan berdasarkan: - Efek utama dari masing-masing variabel - Interaksi antar variabel

Sebagai contoh, dalam tabel kontingensi tiga arah (misalnya: jenis kelamin, status merokok, dan penyakit paru-paru), model log-linier dapat membantu mengidentifikasi apakah: - Interaksi dua variabel cukup menjelaskan struktur data, - Atau diperlukan interaksi tiga arah untuk menangkap keseluruhan asosiasi.

Penyesuaian model dapat dilakukan menggunakan likelihood ratio test, yang digunakan untuk membandingkan model sederhana dengan model yang lebih kompleks. Pendekatan ini membantu dalam memilih model yang paling sesuai berdasarkan keseimbangan antara kesederhanaan dan kecocokan terhadap data.

Model Log-Linier

Model log-linier memodelkan logaritma dari ekspektasi frekuensi sel dalam tabel kontingensi:

\[ \log(\mu_{ij}) = \lambda + \lambda^A_i + \lambda^B_j + \lambda^{AB}_{ij} \]

Cocok untuk menyelidiki asosiasi dan independensi antar variabel.
Tidak membedakan antara variabel respon dan penjelas.
Umumnya digunakan pada tabel multi-dimensi (2 arah, 3 arah, dst).

Model Regresi Logistik

Model regresi logistik biner:

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

Digunakan jika ada variabel dependen kategorik (biasanya biner).
Bertujuan untuk memprediksi probabilitas suatu outcome.
Umumnya digunakan dalam studi observasional atau eksperimental.

14.1 Tabel Kontingensi dan Loglinier

Tabel kontingensi menyajikan frekuensi dari kombinasi kategori antar dua atau lebih variabel. Misalnya:

# Contoh tabel 2x2
matrix(c(45, 15, 60, 80), nrow = 2,
       dimnames = list(Vaksin = c("Sudah", "Belum"),
                       Sakit = c("Ya", "Tidak")))

##        Sakit
## Vaksin  Ya Tidak
##   Sudah 45    60
##   Belum 15    80

Model log-linier untuk tabel I x J dapat dituliskan sebagai:

\[ \log(\mu_{ij}) = \lambda + \lambda^V_i + \lambda^S_j + \lambda^{VS}_{ij} \]

Dengan: - \(\mu_{ij}\): ekspektasi frekuensi untuk sel baris ke-\(i\) dan kolom ke-\(j\) - \(\lambda\): intercept (efek umum) - \(\lambda^V_i\): efek utama dari kategori vaksinasi - \(\lambda^S_j\): efek utama dari status sakit - \(\lambda^{VS}_{ij}\): interaksi antara vaksinasi dan sakit

14.2 Model Saturated

Model saturated atau model penuh menyertakan seluruh efek utama dan interaksi:

Cocok sempurna terhadap data
Tidak mengasumsikan independensi antar variabel Contoh formulasi untuk tabel 2x2:

library(MASS)
data1 <- matrix(c(45, 15, 60, 80), nrow = 2)
dimnames(data1) = list(Vaksin = c("Sudah", "Belum"),Sakit = c("Ya", "Tidak"))
ftable(data1)

##        Sakit Ya Tidak
## Vaksin               
## Sudah        45    60
## Belum        15    80

Model Saturated

model_saturated <- loglm(~ Vaksin * Sakit, data = data1)
summary(model_saturated)

## Formula:
## ~Vaksin * Sakit
## attr(,"variables")
## list(Vaksin, Sakit)
## attr(,"factors")
##        Vaksin Sakit Vaksin:Sakit
## Vaksin      1     0            1
## Sakit       0     1            1
## attr(,"term.labels")
## [1] "Vaksin"       "Sakit"        "Vaksin:Sakit"
## 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

14.3 Model Indepen

Model independen mengasumsikan bahwa tidak ada interaksi antara variabel dalam tabel kontingensi.

Model log-linier yang digunakan untuk menyatakan independensi dua variabel \(X\) dan \(Y\) dalam tabel dua arah (I x J) adalah:

\[ \log(\mu_{ij}) = \mu + \lambda^T_i + \lambda^R_j \]

model_indep <- loglm(~ Vaksin + Sakit, data = data1)
summary(model_indep)

## Formula:
## ~Vaksin + Sakit
## attr(,"variables")
## list(Vaksin, Sakit)
## attr(,"factors")
##        Vaksin Sakit
## Vaksin      1     0
## Sakit       0     1
## attr(,"term.labels")
## [1] "Vaksin" "Sakit" 
## 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 18.06418  1 2.135823e-05
## Pearson          17.40064  1 3.027233e-05

14.4 Odds Ratio

Odds Ratio untuk Tabel 2x2

Odds Ratio (OR) adalah ukuran asosiasi antara dua variabel kategorik dalam tabel kontingensi 2x2. Rumus umum untuk OR adalah:

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

Dengan: - \(n_{11}\) = frekuensi pada baris pertama kolom pertama - \(n_{12}\) = baris pertama kolom kedua - \(n_{21}\) = baris kedua kolom pertama - \(n_{22}\) = baris kedua kolom kedua

Interpretasi Nilai OR - OR = 1: Tidak ada asosiasi antara variabel baris dan kolom (independen) - OR > 1: Terdapat asosiasi positif (kemungkinan kejadian lebih besar di kelompok pertama) - OR < 1: Terdapat asosiasi negatif (kemungkinan kejadian lebih kecil di kelompok pertama)

14.5 estimasi parameter

Dalam model saturated:

Estimasi dilakukan dengan pembatasan seperti sum-to-zero
Estimasi parameter dilakukan dengan iterative proportional fitting (IPF)

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

## [1] 1.386294

14.6 Perbandingan Model

Perbandingan antar model dilakukan dengan menggunakan statistik deviance (G²) atau likelihood ratio test.

anova(model_indep, model_saturated)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Vaksin + Sakit 
## Model 2:
##  ~Vaksin * Sakit 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   18.06418  1                                    
## Model 2    0.00000  0   18.06418         1          2e-05
## Saturated  0.00000  0    0.00000         0          1e+00

15 Model Log Linear 2 Arah

Model log-linear merupakan pendekatan statistik yang digunakan untuk mengevaluasi hubungan antara dua atau lebih variabel kategorik yang disajikan dalam bentuk tabel kontingensi. Model ini bekerja dengan mengasumsikan bahwa logaritma dari nilai harapan frekuensi pada setiap sel dalam tabel (\(\mu_{ij}\)) dapat dijelaskan sebagai hasil penjumlahan dari efek masing-masing variabel, serta (jika diperlukan) interaksi antara variabel-variabel tersebut.

Perbedaan Log-Linear dan Regresi Logistik - Model log-linear digunakan untuk memodelkan frekuensi (count) pada tabel kontingensi dan menguji asosiasi antar variabel kategorik, tanpa menganggap ada variabel respon dan prediktor.

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

15.1 Estimasi Parameter

Sistem Persamaan Model Log-Linear

Persamaan model log-linear untuk tabel 2x2 dapat dituliskan sebagai:

\[ \begin{aligned} \log(\mu_{11}) &= \lambda + \lambda^A_1 + \lambda^B_1 + \lambda^{AB}_{11} \\ \log(\mu_{12}) &= \lambda + \lambda^A_1 + \lambda^B_2 + \lambda^{AB}_{12} \\ \log(\mu_{21}) &= \lambda + \lambda^A_2 + \lambda^B_1 + \lambda^{AB}_{21} \\ \log(\mu_{22}) &= \lambda + \lambda^A_2 + \lambda^B_2 + \lambda^{AB}_{22} \end{aligned} \]

Constraint Sum-to-Zero

Untuk memastikan estimasi parameter yang unik, digunakan kendala jumlah nol (sum-to-zero constraint):

\[ \lambda^A_1 + \lambda^A_2 = 0 \] \[ \lambda^B_1 + \lambda^B_2 = 0 \] \[ \lambda^{AB}_{11} + \lambda^{AB}_{12} + \lambda^{AB}_{21} + \lambda^{AB}_{22} = 0 \]

Rumus Estimasi Parameter dengan Sum-to-Zero Constraint

Parameter-parameter model dapat diestimasi menggunakan rumus berikut:

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

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

\[ \lambda^{AB}_{12} = \frac{1}{4}[\log \mu_{12} - \log \mu_{11} - \log \mu_{22} + \log \mu_{21}] \]

15.2 Analisis Data Tabel Kontingensi 2x2

seorang peneliti ingin melihat efektivitas sebuah kursus terhadap keberhasilan ujian masuk universitas

Ikut pelatihan|Lulus : 65
Ikut pelatihan|Tidak Lulus : 15
Tidak pelatihan|Lulus : 40
Tidak pelatihan|Tidak Lulus : 30

Misalkan:

A1 = Ikut Pelatihan, A2 = Tidak Ikut
B1 = Lulus, B2 = Tidak Lulus

Data observasi tabel kontingensi:

	Lulus (B1)	Tidak Lulus (B2)
Ikut Pelatihan (A1)	\(n_{11} = 65\)	\(n_{12} = 15\)
Tidak Ikut (A2)	\(n_{21} = 40\)	\(n_{22} = 30\)

Model log-linear untuk data ini mengikuti sistem persamaan:

COnstraint sum-to-zero:

\[ \lambda^A_1 + \lambda^A_2 = 0, \] \[ \quad \lambda^B_1 + \lambda^B_2 = 0, \] \[ \quad \lambda^{AB}_{11} + \lambda^{AB}_{12} + \lambda^{AB}_{21} + \lambda^{AB}_{22} = 0 \] ### 15.2.1 dengan cara manual Estimasi Parameter Langkah pengerjaan: 1. Hitung rata-rata log frekuensi \[ \lambda = \frac{1}{4} [\log(65) + \log(15) + \log(40) + \log(30)] \]

lamda = (1/4)*(log(65)+log(15)+log(40)+log(30));lamda

## [1] 3.493129

\[ \lambda = 3.493129 \]

# Matriks frekuensi dari data pelatihan
n <- matrix(c(65, 15, 40, 30), nrow = 2, byrow = TRUE)

# Hitung rata-rata log frekuensi sel
log_mean <- mean(log(c(n)))
log_mean

## [1] 3.493129

Efek Utama A (ikut pelatihan) \[ \lambda^{A}_1 = \frac{1}{2} [\log(65) + \log(15) - \log(40) + \log(30)] \]

lambda_A1 <- 0.5 * ((log(65) + log(15)) - (log(40) + log(30)))
lambda_A2 <- -lambda_A1
lambda_A1; lambda_A2

## [1] -0.1038197

## [1] 0.1038197

efek utama B (kelulusan) \[ \lambda^B_1 = \frac{1}{2}\left[(\log 65 + \log 40) - (\log 15 + \log 30)\right] \]

lambda_B1 <- 0.5 * ((log(65) + log(40)) - (log(15) + log(30)))
lambda_B2 <- -lambda_B1
lambda_B1; lambda_B2

## [1] 0.8770096

## [1] -0.8770096

efek interaksi \[ \lambda^{AB}_{12} = \frac{1}{4} [\log(15) - \log(65) - \log(30) + \log(40)] \]

lambda_AB12 <- 0.25 * (log(15) - log(65) - log(30) + log(40))
lambda_AB11 <- -lambda_AB12
lambda_AB21 <- -lambda_AB12
lambda_AB22 <- lambda_AB12

lambda_AB11; lambda_AB12; lambda_AB21; lambda_AB22

## [1] 0.2946637

## [1] -0.2946637

## [1] 0.2946637

## [1] -0.2946637

Hasil : \[ \lambda = 3.493129 \] \[ \lambda^{A}_1 = -0.1038197 \] \[ \lambda^{A}_2 = 0.1038197 \] \[ \lambda^{B}_1 = 0.8770096 \] \[ \lambda^{B}_2 = -0.8770096 \] \[ \lambda^{AB}_{11} = 0.2946637 \] \[ \lambda^{AB}_{12} = -0.0.2946637 \] \[ \lambda^{AB}_{21} = 0.2946637 \] \[ \lambda^{AB}_{22} = -0.0.2946637 \]

Menghitung Odds Ratio \[ OR = \frac{n_{11}n_{22}}{n_{12}n_{21}}=\frac{65\times30}{40\times15}=3.25 \] Log OR :

log(3.25)

## [1] 1.178655

\[ Log(OR)=Log(3.25)=1.78655 \]

Standard Error \[ SE = \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}} = \sqrt{\frac{1}{65} + \frac{1}{15} + \frac{1}{40} + \frac{1}{30}} \]

\[ = \sqrt{0.0154 + 0.0667 + 0.025 + 0.0333} = \sqrt{0.1404} = 0.3747 \] Confidence Interval 95% untuk log(OR): \[ \log(OR) \pm 1.96 \times SE = 1.1787 \pm 1.96 \times 0.3747 \]

\[ = (1.1787 - 0.7344, 1.1787 + 0.7344) \]

\[ = (0.4443, 1.9131) \] Kembali ke skala OR:

\[ \text{Lower} = e^{0.4443} = 1.559 \quad \text{Upper} = e^{1.9131} = 6.776 \]

Sehingga diperoleh hasil akhir:

\[ OR = 3.25 \quad (95\%\ CI: 1.56 - 6.78) \]

15.2.2 fitting model dengan r

tabel <- matrix(c(65, 15, 40, 30), nrow = 2, byrow = TRUE)
colnames(tabel) <- c("Lulus", "Tidak Lulus")
rownames(tabel) <- c("Ikut", "Tidak")
tabel

##       Lulus Tidak Lulus
## Ikut     65          15
## Tidak    40          30

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

#model tanpa interaksi
fit_no_inter <- glm(Freq ~ Pelatihan + Status, family = poisson, data = data)
summary(fit_no_inter)

## 
## Call:
## glm(formula = Freq ~ Pelatihan + Status, family = poisson, data = data)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         4.0254     0.1239  32.482  < 2e-16 ***
## PelatihanTidak     -0.1335     0.1637  -0.816    0.415    
## StatusTidak Lulus  -0.8473     0.1782  -4.755 1.98e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 35.792  on 3  degrees of freedom
## Residual deviance: 10.440  on 1  degrees of freedom
## AIC: 37.787
## 
## Number of Fisher Scoring iterations: 4

#model dengan interaksi
fit_inter <- glm(Freq ~ Pelatihan * Status, family = poisson, data = data)
summary(fit_inter)

## 
## Call:
## glm(formula = Freq ~ Pelatihan * Status, family = poisson, data = data)
## 
## Coefficients:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                        4.1744     0.1240  33.655  < 2e-16 ***
## PelatihanTidak                    -0.4855     0.2010  -2.416  0.01569 *  
## StatusTidak Lulus                 -1.4663     0.2864  -5.119 3.07e-07 ***
## PelatihanTidak:StatusTidak Lulus   1.1787     0.3747   3.146  0.00166 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  3.5792e+01  on 3  degrees of freedom
## Residual deviance: -8.8818e-15  on 0  degrees of freedom
## AIC: 29.347
## 
## Number of Fisher Scoring iterations: 3

Interpretasi : - Parameter utama (intercept) menunjukkan rata-rata log frekuensi sel.

Efek “Pelatihan” dan “Status” menunjukkan perbedaan log frekuensi antar kategori.
Interaksi signifikan menunjukkan adanya asosiasi antara Pelatihan dan Status.
Nilai log(3.25) =1.787 itu sama dengan efek interaksi ouput R

15.3 Bentuk Model Log-Linear untuk Tabel 2x3

dilakukan survei untuk mengetahui apakah Olahraga dapat membantu meningkatkan tinggi

tabel2x3 <- matrix(c(16, 18, 6, 8, 22, 10), nrow = 2, byrow = TRUE)
colnames(tabel2x3) <- c("Tinggi", "Menengah", "Rendah")
rownames(tabel2x3) <- c("Olahraga", "Tidak Olahraga")
tabel2x3

##                Tinggi Menengah Rendah
## Olahraga           16       18      6
## Tidak Olahraga      8       22     10

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

\[ \log(\mu_{ij}) = \lambda + \lambda^C_i + \lambda^D_j + \lambda^{CD}_{ij} \]

dengan:

\((\mu_{ij})\): ekspektasi frekuensi pada baris ke-\(i\), kolom ke-\(j\)
C: Status Olaharaga \((i = 1 : Olahraga, i = 2 : Tidak~Olahraga)\)
D: Kategori Tinggi \((j = 1 : Tinggi, j = 2 : Menengah, j = 3 : Rendah)\)
Constraint: \(\sum_i \lambda^C_i = 0 , \sum_j \lambda^D_j = 0 , \sum_i \lambda^{CD}_{ij} = 0 , dan \sum_j \lambda^{CD}_{ij} = 0\)

Secara eksplisit:

\[ \log(\mu_{ij}) = \lambda + \lambda^C_1 \text{ (olahraga)}, \lambda^C_2 \text{ (Tidak~Olahraga)} + \lambda^D_1 \text{ (Tinggi)}, \lambda^D_2 \text{ (Menengah)}, \lambda^D_3 \text{ (Rendah)} + \lambda^{CD}_{ij} \text{ (interaksi jika ada)} \]

15.3.1 Fitting Model

data2x3 <- as.data.frame(as.table(tabel2x3))
colnames(data2x3) <- c("Status", "TinggiBadan", "Freq")
data2x3

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

## 
## Call:
## glm(formula = Freq ~ Status + TinggiBadan, family = poisson, 
##     data = data2x3)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           2.485e+00  2.327e-01  10.677   <2e-16 ***
## StatusTidak Olahraga  6.871e-14  2.236e-01   0.000   1.0000    
## TinggiBadanMenengah   5.108e-01  2.582e-01   1.978   0.0479 *  
## TinggiBadanRendah    -4.055e-01  3.227e-01  -1.256   0.2090    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 15.1632  on 5  degrees of freedom
## Residual deviance:  4.1297  on 2  degrees of freedom
## AIC: 38.177
## 
## Number of Fisher Scoring iterations: 4

# Model log-linear dengan interaksi (untuk cek asosiasi)
fit_inter <- glm(Freq ~ Status * TinggiBadan, family = poisson, data = data2x3)
summary(fit_inter)

## 
## Call:
## glm(formula = Freq ~ Status * TinggiBadan, family = poisson, 
##     data = data2x3)
## 
## Coefficients:
##                                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                2.7726     0.2500  11.090   <2e-16
## StatusTidak Olahraga                      -0.6931     0.4330  -1.601   0.1094
## TinggiBadanMenengah                        0.1178     0.3436   0.343   0.7317
## TinggiBadanRendah                         -0.9808     0.4787  -2.049   0.0405
## StatusTidak Olahraga:TinggiBadanMenengah   0.8938     0.5371   1.664   0.0961
## StatusTidak Olahraga:TinggiBadanRendah     1.2040     0.6739   1.787   0.0740
##                                             
## (Intercept)                              ***
## StatusTidak Olahraga                        
## TinggiBadanMenengah                         
## TinggiBadanRendah                        *  
## StatusTidak Olahraga:TinggiBadanMenengah .  
## StatusTidak Olahraga:TinggiBadanRendah   .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1.5163e+01  on 5  degrees of freedom
## Residual deviance: 2.6645e-15  on 0  degrees of freedom
## AIC: 38.048
## 
## Number of Fisher Scoring iterations: 3

Interpretasi - model tanpa interaksi : hasil yang didapat kurang signifikan - Model dengan interaksi : terdapat hasil yang signifikan dan nilai AIC lebih kecil sehingga lebih baik menggunakan model ini - sehingga interaksi antara status berolahraga dan tinggi badan memiliki interaksi

16 Model Log-Linear Tiga Arah

model ini merupakan model yang lebih komplek karena melibatkan tiga variabel kategorik (misal:X,Y,Z), sehingga kemungkinan untuk terjadinya interaksi semakin besar. Berikut beberapa alternatif model log-linear yang umum digunakan

Model Saturated

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk} \]

Model ini memuat semua kemungkinan interaksi, termasuk interaksi tiga arah antara X, Y, dan Z. Model ini sepenuhnya menyesuaikan data dan tidak memiliki derajat bebas.

Model Homogen

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \]

Model ini hanya mencakup interaksi dua arah, tanpa memasukkan interaksi tiga arah. Digunakan untuk mengevaluasi apakah asosiasi tiga arah diperlukan atau cukup dengan asosiasi dua arah.

Model Conditional

Conditional pada X:

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} \] Memuat interaksi X dengan Y dan X dengan Z.

Conditional pada Y:

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk} \] Memuat interaksi Y dengan X dan Y dengan Z.

Conditional pada Z:

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \] Memuat interaksi Z dengan X dan Z dengan Y.

Model Joint Independence

Independensi antara X dan Y:

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} \]

Independensi antara X dan Z:

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} \]

Independensi antara Y dan Z:

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{YZ}_{jk} \]

Model Tanpa Interaksi

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k \]

Model ini hanya mencakup efek utama masing-masing variabel tanpa adanya interaksi. Digunakan sebagai model paling sederhana untuk menguji independensi penuh antar semua variabel.

16.1 Pengujian Interaksi

Dalam analisis model log-linear tiga arah, pengujian interaksi dilakukan untuk mengetahui ada atau tidaknya interaksi antar variabel. Pengujian ini dilakukan secara bertahap, dimulai dari tingkat interaksi tertinggi ke yang lebih rendah. Untuk model log-linear dengan tiga peubah (X, Y, dan Z), tahapan pengujian meliputi:

Pengujian Interaksi Tiga Arah (XYZ):

Bandingkan model saturated dengan model homogenous.

Pengujian Interaksi Dua Arah (XY, XZ, YZ):

Bandingkan model homogenous dengan model conditional.
Bandingkan model conditional dengan model joint independence.
Bandingkan model joint independence dengan model tanpa interaksi.

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

Biasanya, uji dilakukan menggunakan likelihood ratio test (G² test), dan pemilihan model terbaik mempertimbangkan nilai G², derajat bebas (df), dan signifikansi (p-value).

16.2 Analisis Log Linear Tiga Arah

seorang dokter melakukan penelitian untuk melihat pengaruh jenis vaksin dan usia terhadap efektivitasnya dalam menahan infeksi virus

library("epitools")
library("DescTools")
library("lawstat")

## Warning: package 'lawstat' was built under R version 4.4.3

## 
## Attaching package: 'lawstat'

## The following object is masked from 'package:car':
## 
##     levene.test

z.vak <- factor(rep(c("Pfizer", "Moderna", "Sinovac"), each = 4))
x.us   <- factor(rep(c("Dewasa", "Lansia"), each = 2, times = 3))
y.inf <- factor(rep(c("Terinfeksi", "Tidak"), times = 6))

# Frekuensi observasi fiktif berdasarkan kejadian
frekuensi <- c(25, 75, 30, 70, 40, 60, 50, 50, 35, 65, 45, 55)

# Buat data frame
data <- data.frame(
  Vaksin    = z.vak,
  Usia      = x.us,
  Infeksi   = y.inf,
  Frekuensi = frekuensi
)

data

Membentuk tabel kontingensi 3 arah

table3d <- xtabs(Frekuensi ~ Vaksin + Usia + Infeksi, data = data)
ftable(table3d)

##                Infeksi Terinfeksi Tidak
## Vaksin  Usia                           
## Moderna Dewasa                 40    60
##         Lansia                 50    50
## Pfizer  Dewasa                 25    75
##         Lansia                 30    70
## Sinovac Dewasa                 35    65
##         Lansia                 45    55

16.3 Uji Model Interaksi Arah (Saturated VS Homogen)

x.us  <- relevel(x.us, ref = "Lansia")
y.inf  <- relevel(y.inf, ref = "Tidak")
z.vak <- relevel(z.vak, ref = "Sinovac")

16.3.1 Model Saturated

model_saturated <- glm(frekuensi ~ x.us + y.inf + z.vak +
             x.us*y.inf + x.us*z.vak + y.inf*z.vak +
             x.us*y.inf*z.vak,
             family = poisson(link = "log"))
summary(model_saturated)

## 
## Call:
## glm(formula = frekuensi ~ x.us + y.inf + z.vak + x.us * y.inf + 
##     x.us * z.vak + y.inf * z.vak + x.us * y.inf * z.vak, family = poisson(link = "log"))
## 
## Coefficients:
##                                         Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              4.00733    0.13484  29.719   <2e-16
## x.usDewasa                               0.16705    0.18321   0.912   0.3619
## y.infTerinfeksi                         -0.20067    0.20101  -0.998   0.3181
## z.vakModerna                            -0.09531    0.19540  -0.488   0.6257
## z.vakPfizer                              0.24116    0.18019   1.338   0.1808
## x.usDewasa:y.infTerinfeksi              -0.41837    0.29045  -1.440   0.1497
## x.usDewasa:z.vakModerna                  0.01527    0.26502   0.058   0.9541
## x.usDewasa:z.vakPfizer                  -0.09806    0.24736  -0.396   0.6918
## y.infTerinfeksi:z.vakModerna             0.20067    0.28356   0.708   0.4791
## y.infTerinfeksi:z.vakPfizer             -0.64663    0.29669  -2.179   0.0293
## x.usDewasa:y.infTerinfeksi:z.vakModerna  0.01290    0.40746   0.032   0.9747
## x.usDewasa:y.infTerinfeksi:z.vakPfizer   0.16705    0.43048   0.388   0.6980
##                                            
## (Intercept)                             ***
## x.usDewasa                                 
## y.infTerinfeksi                            
## z.vakModerna                               
## z.vakPfizer                                
## x.usDewasa:y.infTerinfeksi                 
## x.usDewasa:z.vakModerna                    
## x.usDewasa:z.vakPfizer                     
## y.infTerinfeksi:z.vakModerna               
## y.infTerinfeksi:z.vakPfizer             *  
## x.usDewasa:y.infTerinfeksi:z.vakModerna    
## x.usDewasa:y.infTerinfeksi:z.vakPfizer     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  5.6788e+01  on 11  degrees of freedom
## Residual deviance: -8.8818e-15  on  0  degrees of freedom
## AIC: 92.436
## 
## Number of Fisher Scoring iterations: 3

exp(model_saturated$coefficients)

##                             (Intercept)                              x.usDewasa 
##                              55.0000000                               1.1818182 
##                         y.infTerinfeksi                            z.vakModerna 
##                               0.8181818                               0.9090909 
##                             z.vakPfizer              x.usDewasa:y.infTerinfeksi 
##                               1.2727273                               0.6581197 
##                 x.usDewasa:z.vakModerna                  x.usDewasa:z.vakPfizer 
##                               1.0153846                               0.9065934 
##            y.infTerinfeksi:z.vakModerna             y.infTerinfeksi:z.vakPfizer 
##                               1.2222222                               0.5238095 
## x.usDewasa:y.infTerinfeksi:z.vakModerna  x.usDewasa:y.infTerinfeksi:z.vakPfizer 
##                               1.0129870                               1.1818182

Model yang digunakan adalah model log-linear saturated dengan semua efek utama, interaksi dua arah, dan interaksi tiga arah. Model ini memodelkan hubungan antara usia (x.us), status terinfeksi (y.inf), dan jenis vaksin (z.vak) terhadap frekuensi responden.

Interpretasi Koefisien - (Intercept) = 4.007 → log(ekspektasi frekuensi) untuk kombinasi referensi (Lansia, Tidak Terinfeksi, Sinovac)

x.usDewasa = 0.167 → odds ratio ≈ 1.18:
- Dewasa memiliki peluang sedikit lebih tinggi dibanding Lansia untuk kombinasi referensi, tetapi tidak signifikan (p = 0.3619)
y.infTerinfeksi = -0.201 → odds ratio ≈ 0.82:
- Peluang terinfeksi sedikit lebih kecil dibanding tidak terinfeksi (tidak signifikan, p = 0.3181)
z.vakndModerna = -0.095 dan z.vakPfizer = 0.241 → dibandingkan Sinovac:
- Moderna ≈ sama peluangnya (OR ≈ 0.91, p = 0.6257)
- Pfizer cenderung meningkatkan peluang (OR ≈ 1.27), tetapi juga tidak signifikan (p = 0.1808)
x.usDewasa:y.infTerinfeksi = -0.418 → OR ≈ 0.66
- Interaksi negatif, artinya efek usia dewasa terhadap infeksi lebih kecil, tapi tidak signifikan (p = 0.1497)
y.infTerinfeksi:z.vakPfizer = -0.647 → OR ≈ 0.52
- Peluang terinfeksi pada vaksin Pfizer lebih kecil dibanding Sinovac saat dikombinasikan (signifikan, p = 0.0293*)
Koefisien interaksi tiga arah tidak signifikan (p > 0.6), menunjukkan bahwa interaksi kompleks antara ketiga variabel tidak berpengaruh secara signifikan terhadap frekuensi sel.

Kesimpulan

Satu-satunya efek yang signifikan secara statistik adalah interaksi antara status infeksi dan vaksin Pfizer (\(p = 0.0293\)), yang menunjukkan bahwa penerima vaksin Pfizer cenderung memiliki peluang lebih kecil untuk terinfeksi.
Tidak ditemukan interaksi yang signifikan pada efek tiga arah, dan sebagian besar koefisien dua arah serta utama tidak signifikan.
Model memberikan deviance residual yang sangat kecil (mendekati nol), menandakan model sangat cocok dengan data (kemungkinan karena model saturated).
AIC rendah (92.436) dan deviance nol mengindikasikan model saturated yang mengakomodasi seluruh variasi data, sehingga tidak ada sisa deviance untuk diuji lebih lanjut.

16.3.2 Model Homogenous

# Homogenous Model
model_homogenous <- glm(frekuensi ~ x.us + y.inf + z.vak +
              x.us*y.inf + x.us*z.vak + y.inf*z.vak,
              family = poisson(link = "log"))
summary(model_homogenous)

## 
## Call:
## glm(formula = frekuensi ~ x.us + y.inf + z.vak + x.us * y.inf + 
##     x.us * z.vak + y.inf * z.vak, family = poisson(link = "log"))
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   4.01874    0.12441  32.303  < 2e-16 ***
## x.usDewasa                    0.14589    0.15757   0.926  0.35453    
## y.infTerinfeksi              -0.22620    0.16681  -1.356  0.17508    
## z.vakModerna                 -0.09690    0.17074  -0.568  0.57036    
## z.vakPfizer                   0.21356    0.16349   1.306  0.19146    
## x.usDewasa:y.infTerinfeksi   -0.36521    0.17151  -2.129  0.03323 *  
## x.usDewasa:z.vakModerna       0.01835    0.20100   0.091  0.92724    
## x.usDewasa:z.vakPfizer       -0.04582    0.20190  -0.227  0.82047    
## y.infTerinfeksi:z.vakModerna  0.20647    0.20340   1.015  0.31007    
## y.infTerinfeksi:z.vakPfizer  -0.56810    0.21508  -2.641  0.00826 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 56.78787  on 11  degrees of freedom
## Residual deviance:  0.18161  on  2  degrees of freedom
## AIC: 88.618
## 
## Number of Fisher Scoring iterations: 3

menguji apakah ada interaksi tiga arah - Hipotesis \(H_0\) : Tidak ada interaksi tiga arah (model homogenous sudah cukup) \(H_1\) : Ada interaksi tiga arah (model saturated diperlukan)

Statistik Uji

# menghitung selisih deviance
Deviance.model <- model_homogenous$deviance - model_saturated$deviance
Deviance.model

## [1] 0.181611

# menghitung derajat bebas
derajat.bebas <- (model_homogenous$df.residual - model_saturated$df.residual)
derajat.bebas

## [1] 2

# menentukan Chi-Square tabel
chi.tabel <- qchisq(1 - 0.05, df = derajat.bebas)
chi.tabel

## [1] 5.991465

Keputusan

Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima H0", "Tolak H0")
Keputusan

## [1] "Terima H0"

• Hipotesis - \(H_0\): \(\lambda^{XYZ}_{ijk} = 0\) (Tidak ada interaksi tiga arah; model yang terbentuk adalah model homogenous) - \(H_1\): \(\lambda^{XYZ}_{ijk} \neq 0\) (Ada interaksi tiga arah; model yang terbentuk adalah model saturated)

• Tingkat Signifikansi - \(\alpha = 5\%\)

• Statistik Uji - \(\Delta\)Deviance = Deviance model homogenous \(-\) Deviance model saturated - \(= 0.181611 - 0 = 0.181611\) - \(df = df_{\text{homogenous}} - df_{\text{saturated}} = 2 - 0 = 2\)

• Daerah Penolakan - Tolak \(H_0\) jika \(\Delta\)Deviance \(> \chi^2_{0.05, df} = \chi^2_{0.05, 2} = 5.991\)

• Keputusan - Karena \(0.181611 < 5.991\), maka tidak tolak \(H_0\)

• Interpretasi - Pada taraf nyata 5%, belum cukup bukti untuk menolak \(H_0\) atau dapat dikatakan bahwa tidak ada interaksi tiga arah antara usia, jenis vaksin, dan status infeksi.

16.4 UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON X)

# Conditional Association on X
model_conditional_X <- glm(frekuensi ~ x.us + y.inf + z.vak +
              x.us*y.inf + x.us*z.vak,
              family = poisson(link = "log"))
summary(model_conditional_X)

## 
## Call:
## glm(formula = frekuensi ~ x.us + y.inf + z.vak + x.us * y.inf + 
##     x.us * z.vak, family = poisson(link = "log"))
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 4.066e+00  1.113e-01  36.543  < 2e-16 ***
## x.usDewasa                  1.335e-01  1.551e-01   0.861  0.38919    
## y.infTerinfeksi            -3.365e-01  1.171e-01  -2.873  0.00406 ** 
## z.vakModerna                2.493e-16  1.414e-01   0.000  1.00000    
## z.vakPfizer                 8.100e-15  1.414e-01   0.000  1.00000    
## x.usDewasa:y.infTerinfeksi -3.567e-01  1.695e-01  -2.105  0.03530 *  
## x.usDewasa:z.vakModerna    -2.220e-16  2.000e-01   0.000  1.00000    
## x.usDewasa:z.vakPfizer     -6.905e-15  2.000e-01   0.000  1.00000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 56.788  on 11  degrees of freedom
## Residual deviance: 14.436  on  4  degrees of freedom
## AIC: 98.872
## 
## Number of Fisher Scoring iterations: 4

• Hipotesis - \(H_0\): \(\lambda^{YZ}_{jk} = 0\) (Tidak ada interaksi antara status terinfeksi (Y) dan jenis vaksin (Z)) - \(H_1\): \(\lambda^{YZ}_{jk} \neq 0\) (Ada interaksi antara status terinfeksi (Y) dan jenis vaksin (Z))

• Tingkat Signifikansi - \(\alpha = 5\%\)

• Statistik Uji - \(\Delta\)Deviance = Deviance model conditional on X \(-\) Deviance model homogenous - \(= 14.436 - 0.18161 = 14.2541\) - \(df = df_{\text{conditional on X}} - df_{\text{homogenous}} = 4 - 2 = 2\)

• Daerah Penolakan - Tolak \(H_0\) jika \(\Delta\)Deviance > \(\chi^2_{0.05,2} = 5.991\)

• Keputusan - Karena \(14.2541 > 5.991\), maka tolak \(H_0\)

• Kesimpulan - Dengan taraf nyata 5%, terdapat bukti yang cukup untuk menolak \(H_0\) atau dapat dikatakan bahwa ada interaksi antara pendapat tentang hukuman mati dan fundamentalisme. Dengan kata lain, model yang terbentuk adalah model dengan parameter \(\lambda^{YZ}_{jk}\).

Pengujian Selisih Deviance (Conditional on X vs Homogenous)

# Deviance of Model
Deviance.model <- model_conditional_X$deviance - model_homogenous$deviance   # model_conditional_X: conditional on X, model_homogenous: homogenous
Deviance.model

## [1] 14.2541

Hitung Derajat Bebas

# Chi Square tabel dengan alpha = 0.05
derajat.bebas <- (4 - 2)
derajat.bebas

## [1] 2

Nilai Chi-Square Tabel

chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel

## [1] 5.991465

Keputusan Uji

Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan

## [1] "Tolak"

16.5 UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON Y)

# Conditional Association on Y
model_conditional_Y <- glm(frekuensi ~ x.us + y.inf + z.vak +
              x.us*y.inf + y.inf*z.vak,
              family = poisson(link = "log"))
summary(model_conditional_Y)

## 
## Call:
## glm(formula = frekuensi ~ x.us + y.inf + z.vak + x.us * y.inf + 
##     y.inf * z.vak, family = poisson(link = "log"))
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   4.02535    0.10668  37.732  < 2e-16 ***
## x.usDewasa                    0.13353    0.10351   1.290  0.19704    
## y.infTerinfeksi              -0.23111    0.16564  -1.395  0.16294    
## z.vakModerna                 -0.08701    0.13200  -0.659  0.50978    
## z.vakPfizer                   0.18924    0.12341   1.533  0.12517    
## x.usDewasa:y.infTerinfeksi   -0.35667    0.16945  -2.105  0.03530 *  
## y.infTerinfeksi:z.vakModerna  0.20479    0.20257   1.011  0.31203    
## y.infTerinfeksi:z.vakPfizer  -0.56394    0.21427  -2.632  0.00849 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 56.78787  on 11  degrees of freedom
## Residual deviance:  0.28753  on  4  degrees of freedom
## AIC: 84.724
## 
## Number of Fisher Scoring iterations: 3

• Hipotesis - \(H_0\): \(\lambda^{XZ}_{jk} = 0\) (Tidak ada interaksi antara usia (X) dan jenis vaksin (Z)) - \(H_1\): \(\lambda^{XZ}_{jk} \neq 0\) (Ada interaksi antara usia (X) dan jenis vaksin (Z))

• Tingkat Signifikansi - \(\alpha = 5\%\)

• Statistik Uji - \(\Delta\)Deviance = Deviance model conditional on Y \(-\) Deviance model homogenous - \(= 0.28753 - 0.18161 = 0.1059176\) - \(df = df_{\text{conditional on Y}} - df_{\text{homogenous}} = 4 - 2 = 2\)

• Daerah Penolakan - Tolak \(H_0\) jika \(\Delta\)Deviance > \(\chi^2_{0.05,2} = 5.991\)

• Keputusan - Karena \(0.1059176 < 5.991\), maka terima \(H_0\)

• Kesimpulan - Dengan taraf nyata 5%, belum cukup bukti untuk menolak \(H_0\) atau dapat dikatakan bahwa tidak ada interaksi antara usia dan jenis vaksin. Dengan kata lain, model yang terbentuk adalah model tanpa parameter \(\lambda^{XZ}_{jk}\).

Pengujian Hipotesis Interaksi X dan Z (Conditional on Y vs Homogenous)

# Deviance of Model
Deviance.model <- model_conditional_Y$deviance - model_homogenous$deviance   # model_conditional_Y: conditional on Y, model_homogenous: homogenous
Deviance.model

## [1] 0.1059176

Hitung Derajat Bebas

derajat.bebas <- (4 - 2)
derajat.bebas

## [1] 2

Nilai Chi-Square Tabel

chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel

## [1] 5.991465

Keputusan Uji

Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan

## [1] "Terima"

16.6 UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON Z)

# Conditional Association on Z
model_conditional_Z <- glm(frekuensi ~ x.us + y.inf + z.vak +
              x.us*z.vak + y.inf*z.vak,
              family = poisson(link = "log"))
summary(model_conditional_Z)

## 
## Call:
## glm(formula = frekuensi ~ x.us + y.inf + z.vak + x.us * z.vak + 
##     y.inf * z.vak, family = poisson(link = "log"))
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   4.094e+00  1.155e-01  35.458  < 2e-16 ***
## x.usDewasa                   -1.854e-16  1.414e-01   0.000  1.00000    
## y.infTerinfeksi              -4.055e-01  1.443e-01  -2.809  0.00497 ** 
## z.vakModerna                 -8.701e-02  1.656e-01  -0.525  0.59929    
## z.vakPfizer                   1.892e-01  1.588e-01   1.191  0.23349    
## x.usDewasa:z.vakModerna       8.400e-17  2.000e-01   0.000  1.00000    
## x.usDewasa:z.vakPfizer       -5.421e-16  2.000e-01   0.000  1.00000    
## y.infTerinfeksi:z.vakModerna  2.048e-01  2.026e-01   1.011  0.31203    
## y.infTerinfeksi:z.vakPfizer  -5.639e-01  2.143e-01  -2.632  0.00849 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 56.788  on 11  degrees of freedom
## Residual deviance:  4.739  on  3  degrees of freedom
## AIC: 91.175
## 
## Number of Fisher Scoring iterations: 4

• Hipotesis - \(H_0\): \(\lambda^{XY}_{jk} = 0\) (Tidak ada interaksi antara usia (X) dan status terinfeksi(Y)) - \(H_1\): \(\lambda^{XY}_{jk} \neq 0\) (Ada interaksi antara usia (X) dan status terinfeksi(Y))

• Tingkat Signifikansi - \(\alpha = 5\%\)

• Statistik Uji - \(\Delta\)Deviance = Deviance model conditional on Z \(-\) Deviance model homogenous - \(= 4.739 - 0.18161 = 4.557343\) - \(df = df_{\text{conditional on Z}} - df_{\text{homogenous}} = 3 - 2 = 1\)

• Daerah Penolakan - Tolak \(H_0\) jika \(\Delta\)Deviance > \(\chi^2_{0.05,2} = 3.841459\)

• Keputusan - Karena \(4.557343 > 3.841459\), maka tolak \(H_0\)

• Kesimpulan - Dengan taraf nyata 5%, terdapat cukup bukti untuk menolak \(H_0\) atau dapat dikatakan bahwa ada interaksi antara usia dan status terinfeksi. Dengan kata lain, model yang terbentuk adalah model dengan parameter \(\lambda^{XY}_{jk}\).

# Deviance of Model
Deviance.model <- model_conditional_Z$deviance - model_homogenous$deviance # model_conditional_Z
Deviance.model

## [1] 4.557343

Hitung Derajat Bebas

derajat.bebas <- (3 - 2)
derajat.bebas

## [1] 1

Nilai Chi-Square Tabel

chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel

## [1] 3.841459

Keputusan Uji

Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan

## [1] "Tolak"

16.7 Pemilihan Model Terbaik

Model	Parameter	Deviance	Jumlah Parameter	df	AIC
Saturated	\(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk}\)	0.00	12	0	92.436
Homogenous	\(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk}\)	0.18161	10	2	88.618
Conditional on X	\(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik}\)	14.436	8	4	98.872
Conditional on Y	\(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk}\)	0.28753	8	4	84.724
Conditional on Z	\(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk}\)	4.739	9	3	91.175

kesimpulan dari hasil diketahui bahwa Conditional on X vs Homogenous dan Conditional on Z vs Homogenous memiliki hasil yang signifikan oleh karena itu dipilih lah model saturated sebagai model terbaik

16.8 Model Terbaik

bestmodel <- model_saturated
summary(bestmodel)

## 
## Call:
## glm(formula = frekuensi ~ x.us + y.inf + z.vak + x.us * y.inf + 
##     x.us * z.vak + y.inf * z.vak + x.us * y.inf * z.vak, family = poisson(link = "log"))
## 
## Coefficients:
##                                         Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              4.00733    0.13484  29.719   <2e-16
## x.usDewasa                               0.16705    0.18321   0.912   0.3619
## y.infTerinfeksi                         -0.20067    0.20101  -0.998   0.3181
## z.vakModerna                            -0.09531    0.19540  -0.488   0.6257
## z.vakPfizer                              0.24116    0.18019   1.338   0.1808
## x.usDewasa:y.infTerinfeksi              -0.41837    0.29045  -1.440   0.1497
## x.usDewasa:z.vakModerna                  0.01527    0.26502   0.058   0.9541
## x.usDewasa:z.vakPfizer                  -0.09806    0.24736  -0.396   0.6918
## y.infTerinfeksi:z.vakModerna             0.20067    0.28356   0.708   0.4791
## y.infTerinfeksi:z.vakPfizer             -0.64663    0.29669  -2.179   0.0293
## x.usDewasa:y.infTerinfeksi:z.vakModerna  0.01290    0.40746   0.032   0.9747
## x.usDewasa:y.infTerinfeksi:z.vakPfizer   0.16705    0.43048   0.388   0.6980
##                                            
## (Intercept)                             ***
## x.usDewasa                                 
## y.infTerinfeksi                            
## z.vakModerna                               
## z.vakPfizer                                
## x.usDewasa:y.infTerinfeksi                 
## x.usDewasa:z.vakModerna                    
## x.usDewasa:z.vakPfizer                     
## y.infTerinfeksi:z.vakModerna               
## y.infTerinfeksi:z.vakPfizer             *  
## x.usDewasa:y.infTerinfeksi:z.vakModerna    
## x.usDewasa:y.infTerinfeksi:z.vakPfizer     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  5.6788e+01  on 11  degrees of freedom
## Residual deviance: -8.8818e-15  on  0  degrees of freedom
## AIC: 92.436
## 
## Number of Fisher Scoring iterations: 3

# Interpretasi koefisien model terbaik
data.frame(
  koef = bestmodel$coefficients,
  exp_koef = exp(bestmodel$coefficients)
)

Interpretasi : -`Koefisien utama seperti usia (Dewasa) dan vaksin Pfizer cenderung meningkatkan frekuensi, meskipun sebagian besar tidak signifikan.

Interaksi dua arah yang paling mencolok adalah:
y.infTerinfeksi:z.vakPfizer = 0.52 → menunjukkan bahwa vaksin Pfizer secara signifikan menurunkan frekuensi kasus infeksi dibandingkan Sinovac.
Interaksi tiga arah memiliki efek kecil dan tidak signifikan (OR mendekati 1).

16.9 Nilai Dugaan Model Terbaik

# Fitted values dari model terbaik
data.frame(
  Vaksin = z.vak,
  Usia = x.us,
  Infeksi = y.inf,
  counts = frekuensi,
  fitted = bestmodel$fitted.values
)

Referensi

Fisher, R. A. (1935). The Design of Experiments. Oliver and Boyd.
Howell, D. C. (2012). Statistical Methods for Psychology (8th ed.). Cengage Learning.
Jaya, I. G. N. M. (2025). E-book Analisis Data Kategori.
Agresti, A. (2013). Categorical Data Analysis. John Wiley & Sons.
Sihotang, S. F. (n.d.). Analisis Model Log Linier Tiga Dimensi untuk Data Kualitatif dengan Metode Forward.
Theodoridis, S. (2020). Machine learning: A bayesian and optimization perspective.
Anderson, C. J. (2018). Log‑linear Models for Contingency Tables
Holt, D. (1979). “Log‑Linear Models for Contingency Table Analysis.” Sociological Methods & Research.
Hilbe, J. M. (2009). Logistic Regression Models. Chapman & Hall/CRC.
Khamis, H. J. (1983). “Log‑linear model analysis of the semi‑symmetric intraclass contingency table.” Communications in Statistics – Theory and Methods.