---

---

Puji syukur ke hadirat Tuhan Yang Maha Esa karena atas limpahan rahmat, taufik, dan hidayah-Nya, penulisan buku berjudul “Analisis Data Kategori: Panduan Teori dan Praktikum dalam R Studio” ini dapat diselesaikan dengan baik dan tepat waktu. Buku ini hadir sebagai bentuk kontribusi dalam memperluas literatur statistik, khususnya dalam bidang analisis data kategori yang memiliki peran penting dalam penelitian sosial, kesehatan, dan ilmu terapan lainnya.

Penulisan buku ini dilandasi oleh kebutuhan akan referensi yang tidak hanya membahas teori secara konseptual, tetapi juga disertai dengan panduan praktis yang aplikatif. Oleh karena itu, buku ini dirancang untuk memberikan pemahaman menyeluruh, mulai dari konsep dasar seperti tabel kontingensi dan uji chi-kuadrat, hingga model lanjutan seperti regresi logistik dan model log-linear. Setiap topik disertai contoh penerapan menggunakan software R Studio, sehingga pembaca dapat langsung menghubungkan teori dengan praktik.

Penulis menyadari bahwa belajar statistik, khususnya topik data kategori, sering kali terasa abstrak dan teknis. Oleh sebab itu, pendekatan yang digunakan dalam buku ini adalah menjelaskan dengan bahasa yang lugas, disertai ilustrasi dan contoh kasus yang relevan agar memudahkan pembaca dalam memahami setiap konsep yang disampaikan.

Ucapan terima kasih penulis sampaikan kepada seluruh pihak yang telah memberikan dukungan selama proses penyusunan buku ini. Terima kasih yang sebesar-besarnya kepada para dosen, rekan mahasiswa, dan semua pihak yang telah memberi masukan, semangat, maupun kritik yang membangun.

Penulis juga menyadari bahwa tidak ada karya yang sempurna. Oleh karena itu, kritik dan saran dari para pembaca sangat penulis harapkan untuk perbaikan dan pengembangan di masa mendatang. Semoga buku ini dapat memberikan manfaat bagi pembaca dalam memahami dan menerapkan analisis data kategori, serta menjadi referensi yang berguna dalam pembelajaran dan penelitian.

Dalam era informasi yang berkembang pesat, kemampuan untuk memahami dan mengolah data menjadi semakin penting. Tidak hanya data numerik yang berbentuk angka-angka kuantitatif, namun juga data yang bersifat kualitatif atau kategorikal. Data kategori merupakan jenis data yang umum ditemukan dalam berbagai bidang seperti kesehatan, pendidikan, pemasaran, politik, dan ilmu sosial. Data ini tidak menggambarkan besaran numerik, melainkan klasifikasi atau kelompok, seperti jenis kelamin, status pekerjaan, atau tingkat kepuasan.

Analisis data kategori menjadi salah satu cabang penting dalam ilmu statistik yang bertujuan untuk memahami karakteristik, distribusi, dan hubungan antara variabel-variabel kategori. Pendekatan ini digunakan untuk merangkum informasi, mendeteksi pola, serta mengevaluasi kecenderungan dalam populasi berdasarkan data kategori yang tersedia. Walaupun tidak bersifat numerik dalam arti konvensional, data kategori memiliki struktur dan informasi yang kaya serta dapat memberikan wawasan yang berarti apabila dianalisis dengan tepat.

Buku ini disusun untuk memberikan pemahaman mendalam mengenai konsep dan aplikasi analisis data kategori. Materi disajikan secara sistematis, dimulai dari pengenalan dasar hingga ke pembahasan analisis yang lebih kompleks, sehingga diharapkan mampu memberikan gambaran yang komprehensif kepada pembaca. Penyajian disertai dengan contoh-contoh kontekstual dan aplikatif guna mendekatkan konsep statistik dengan persoalan nyata yang sering ditemui dalam penelitian.

Diharapkan, buku ini dapat menjadi referensi yang bermanfaat dalam kajian-kajian akademik maupun profesional. Dengan memahami prinsip-prinsip dasar dan logika analisis data kategori, pembaca akan memiliki fondasi yang kuat dalam membuat interpretasi dan pengambilan keputusan yang didasarkan pada bukti data yang sahih.

1 ANALISIS DATA KATEGORI

1.1 Definisi Umum

Merupakan analisis yang fokus terhadap data yang diklasifikasikan ke dalam data kategori, atau yang dapat diklasifikasikan kepada data yang bersifat nominal atau ordinal. Contoh data kategori mencakup jenis kelamin, tingkat pendidikan dan lainnya yang bersifat membedakan.

Analisis dilakukan untuk beberapa hal seperti

1. Mengidentifikasi pola dan tren data

2. menganalisis hubungan antar variabel

3. Mengembangkan model

4. dan lainnya

Analisis data kategori berfokus pada bagaimana pola, frekuensi, dan distribusi dari kategori-kategori tersebut muncul dalam suatu kumpulan data. Tujuannya bukan hanya untuk mengetahui seberapa banyak orang memilih satu kategori tertentu, tetapi juga untuk memahami konteks, hubungan, dan makna dari pilihan-pilihan tersebut. Dengan memahami data kategori, kita bisa memperoleh wawasan tentang sikap, karakteristik, dan preferensi individu atau kelompok dalam sebuah populasi.

Pentingnya analisis data kategori terletak pada kemampuannya menjelaskan hal-hal yang tidak bisa direpresentasikan oleh angka secara langsung. Ia menjadi jembatan antara fenomena sosial, psikologis, atau budaya dengan data yang bisa dianalisis secara sistematis. Oleh karena itu, meskipun sederhana dalam bentuknya, data kategori tetap memiliki kedalaman makna yang besar dan layak dipahami secara cermat.

Analisis ini dilakukan untuk memahami data yang tidak bersifat numerik secara langsung, tetapi juga membawa makna besar dalam bentuk klasifikasi, pilihan dan atribut yang menjadi dasar dalam banyak penelitian empiris modern

1.2 Ruang Lingkup

Analisis data kategori ini merupakan analisis yang dilakukan pada data yang bersifat diskrit, yakni data yang berskala nominal atau ordinal. Data nominal merupakan data yang bersifat sebagai pembeda. Misal jenis kelamin, tinggi-rendah dan lainnya sedangkan data ordinal merupakan data yang selain menjadi pembeda juga memiliki tingkatan seperti tingkat pendidikan (SD,SMP,SMA) dan lainnya.

1.3 Manfaat Analisis Data Kategori

Analisis ini berguna di banyak bidang. Sebagai contoh penggunaan dalam bidang sosial yakni survei mengenai hubungan sosial ekonomi dengan tingkat pendidikan. Dalam Bidang psikologi dapat dilakukan untuk data-data kategori psikologi. Dalam bidang kesehatan seperti efektivitas pengobatan terapi dan banyak manfaat lainnya.

1.4 Metode dan Distribusi dalam Analisis Data Kategori

1.4.1 Metode Analisis

Bnayak metode yang dapat digunakan dalam analisis data kategori tergantung pada tujuan analisisnya. Beberapa contohnya adalah

1. Tabel kontingensi

   Metode ini memperlihatkan distribusi frekuensi gabungan antara dua atau lebih variabel kategori. Dengan ini, kita dapat dengan mudah melihat pola hubungan atau kecenderungan antara variabel.

   Dari tabel kontingensi tersebut, analisis dapat dilanjutkan menggunakan uji chi-square (chi-square test of independence). Uji ini berguna untuk mengetahui apakah terdapat hubungan yang signifikan secara statistik antara dua variabel kategori.

2. Regresi Logistik

   Metode ini digunakan untuk memodelkan hubungan antara satu variabel dependen kategori (biasanya berskala binomial seperti “ya” atau “tidak”) dengan satu atau lebih variabel independen yang dapat berupa numerik maupun kategori. Regresi logistik tidak hanya mampu mengidentifikasi faktor-faktor yang memengaruhi suatu kejadian, tetapi juga dapat memprediksi probabilitas terjadinya kejadian tersebut. Bila variabel dependennya memiliki lebih dari dua kategori, maka digunakan regresi logistik multinomial.

3. Analisis Korespondensi

   Metode ini digunakan untuk memvisualisasikan hubungan antara kategori dalam tabel kontingensi secara grafis dalam ruang dua dimensi. Teknik ini sangat bermanfaat untuk menemukan pola dan asosiasi antara kategori dari dua variabel, terutama dalam analisis survei atau data kualitatif lainnya.

4. Random Forest

   alah satu algoritma yang sangat efektif, khususnya untuk tugas klasifikasi, yaitu ketika tujuan analisis adalah menentukan kategori atau label suatu data berdasarkan sejumlah fitur input. Misalnya, memprediksi apakah seseorang akan membeli produk (kategori: ya atau tidak), mengidentifikasi jenis kelamin berdasarkan preferensi konsumen, atau mengklasifikasikan jenis keluhan pelanggan berdasarkan isi pesan.

1.4.2 Distribusi Probabilitas

Dalam analisis inferensi, perlu dipastikan bahwa asumsi mengenai distribusi peluang terpenuhi. Misal dalam model ANOVA, asumsi normalitas mengambil peran sentral. Oleh karena itu, pada analisis data kategoripun demikian

1. Distribusi Binomial

   Seringkali, data kategorikal berasal dari n percobaan yang independen dan memiliki dua kemungkinan hasil yakni berhasil atau gagal. Percobaan yang dilakukan untuk satu individu disebut sebagai percobaan bernouli. Namun ketika dilakukan n kali percobaan dengan setiap individu dalam n bersifat independen, maka variabel acak Y (respon) akan mengikuti distribusi binomial dengan persamaan

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

   Contoh sebuah kuis memiliki 10 soal dengan jawaban benar atau salah. Seorang siswa tidak siap dengan kuis dadakan tersebut dan menjawab secara acak. Peluang menjawab benar adalah 0,2. Maka peluang menjawab 7 benar adalah

   \[ \begin{aligned} P(X = 7) &= \binom{10}{7} (0.2)^7 (1 - 0.2)^{10 - 7} \\ &= \binom{10}{7} (0.2)^7 (0.8)^3 \\ &= 120 \times 0.000128 \times 0.512 \\ &= 120 \times 0.000065536 \\ &= 0.00786432 \\ &\approx {0.0079} \end{aligned} \]

   Distribusi ini adalah distribusi dasar dalam analisis data kategori.

2. Distribusi Multinomial

   Distribusi ini adalah lanjutan dari distribusi binomial. Ketika hasil yang mungkin terjadi lebih dari dua kemungkinan seperti tinggi, sedang, rendah sehingga persamaannya berubah menjadi

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

   Ini merupakan distribusi yang sering muncul dalam kasus analisis data kategori

3. Distribusi Poisson

   Distribusi ini digunakan untuk menghitung sebuah kejadian yang jarang terjadi dalam interval waktu tertentu. Fungsi probabilitasnya adalah

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

   Ini biasanya terjadi pada kejadian yang tidak diharapkan seperti bencana atau hal lainnya.

1.5 Desain Sampling

Desain sampling memiliki peran yang cukup dalam menentukan validitas dan reliabilitas hasil penelitian. Dalam penelitian terdapat beberapa desain sampling. Pengklasifikasiannya adalah sebagai berikut

1. Prospective Sampling

   Merupakan suatu desain sampling dimana variabel respons Y diukur setelah perlakukan atau pengelompokan dikerjakan terlebih dahulu. Pada desain ini peneliti memiliki kendali utnuk menetapkan siapa yang akan masuk grup mana. Contoh dari sampling ini adalah

   1. Clinical Trial (Studi Eksperimen)

      Subjek secara acak dialokasikan ke dalam kelompok perlakuan dan kontrol. Tekniknya meliputi SRS, Stratisfied dan Cluster sampling.

   2. Cohort Study

      Observasional dimana individu di amati dari waktu ke waktu. tekniknya ada sensus, sistematik, matched sampling dan lainnya.

2. Retrospective Sampling

   merupakan metode yang mana data dikumpulkan dari peristiwa yang telah terjadi sebelumnya. alias variabel Y ditentukan terlebih dahulu, kemudian baru variabel X. Contoh dari sampling ini adalah

   1. Studi kasus-kontrol

      Sekelompok individu dengan kasus A dibandingkan dengan sekelompok kasus yang tidak mengalami kasus A. Contoh tekniknya adlaah puposive, snowball, incidence sampling, dan lainnya.

   2. Studi Kohort Restrospektif

      Menggunakan data historis dalam mengelompokkan individu. Contoh tekniknya adalah convenience sampling, quota sampling dan case-based sampling.

3. Cross Sectional

   Desain ini tidak membatas variabel mana yang menjadi prediktor dan mana yang menjadi respon.

2 Tabel Kontingensi 2 Arah

2.1 Definisi Umum

Tabel kontingensi adalah sebuah cara penyajian data dalam bentuk matriks atau tabel yang digunakan untuk menampilkan hubungan antara dua atau lebih variabel kategori. Dalam tabel ini, baris dan kolom merepresentasikan kategori dari masing-masing variabel, dan setiap sel di dalam tabel menunjukkan jumlah kejadian (frekuensi) yang sesuai dengan kombinasi kategori-kategori tersebut. Tabel kontingensi sering digunakan sebagai langkah awal untuk mengenali pola, hubungan, atau ketidakteraturan dalam data kategori.

Sebagai contoh, sebuah tabel kontingensi bisa menunjukkan hubungan antara jenis kelamin (laki-laki/perempuan) dan preferensi warna (merah/biru/hijau). Setiap sel dalam tabel tersebut akan berisi jumlah orang dalam kombinasi kategori tersebut, misalnya berapa banyak laki-laki yang memilih merah, atau perempuan yang memilih biru. Penyajian seperti ini membantu kita untuk melihat distribusi data secara menyeluruh sekaligus membandingkan antar kategori dengan lebih mudah.

Tabel kontingensi menjadi alat penting dalam analisis data kategori karena ia menyederhanakan informasi yang kompleks ke dalam bentuk yang terstruktur dan mudah dibaca. Lewat tabel ini, peneliti bisa mulai menelusuri apakah ada indikasi hubungan antar variabel, kesamaan distribusi, atau perbedaan yang mencolok. Dengan kata lain, tabel kontingensi tidak hanya berfungsi sebagai alat penyajian data, tapi juga sebagai jendela awal untuk memahami interaksi antara kategori-kategori yang ada.

Pada bab ini akan dibahas mengenai tabel kontingensi 2 arah.

Misalnya variabel X memiliki dua kategori: \(x_1\) dan \(x_2\), serta variabel Y memiliki dua kategori: \(y_1\) dan \(y_2\).

	\(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_{..}\)

Keterangan:

• \(x_i\) adalah kategori ke-\(i\) dari variabel X
• \(y_j\) adalah kategori ke-\(j\) dari variabel Y
• \(n_{ij}\) adalah frekuensi sel baris ke-\(i\) kolom ke-\(j\)
• \(n_{i.} = n_{i1} + n_{i2}\) adalah total baris ke-\(i\)
• \(n_{.j} = n_{1j} + n_{2j}\) adalah total kolom ke-\(j\)
• \(n_{..} = n_{1.} + n_{2.} = n_{.1} + n_{.2}\)

Pada tabel kontingensi, apabila X dan Y merupakan variabel acak dan pemilihan sampelnya juga acak, maka disebutkan bahwa tabel tersebut memiliki distribusi peluang.

Pada tabel kontingensi, distribusi peluang terdiri atas :

1. Joint Probability Peluang bersama menunjukkan peluang dua kejadian terjadi bersamaan.

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

2. Marginal Probability Peluang marjinal didapat dari total baris/kolom dibagi total keseluruhan.

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

3. Conditional Probability Peluang bersyarat adalah peluang satu kejadian terjadi jika kejadian lain sudah terjadi.

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

2.2 Contoh Kasus

Misalkan kita punya data tentang kebiasaan olahraga dan cara bangun tidur pada 100 orang. Pada faktor olahraga terdapat dua taraf yakni kebiasaan sering olahraga dan jarang olahraga sedangkan untuk faktor bangun ada dua taraf yakni langsung bangun atau tunda terlebih dahulu selama 5 menit.

	Langsung Bangun	Tunda-Tunda dulu	Total
Olahraga	40	10	50
Tidak	20	30	50
Total	60	40	100

Langkah 1: Peluang Bersama

Hitung:

• \(P(Olahraga \cap Langsung = \frac{40}{100} = 0.4\)
• \(P(Olahraga \cap Tunda) = \frac{10}{100} = 0.1\)
• \(P(Tidak \cap Langsung) = \frac{20}{100} = 0.2\)
• \(P(Tidak \cap Tunda) = \frac{30}{100} = 0.3\)

Langkah 2: Peluang Marjinal

Hitung:

• \(P(Olahraga) = \frac{50}{100} = 0.5\)
• \(P(Tidak Olahraga) = \frac{50}{100} = 0.5\)
• \(P(Langsung) = \frac{60}{100} = 0.6\)
• \(P(Tunda = \frac{40}{100} = 0.4\)

Langkah 3: Peluang Bersyarat

Hitung:

• \(P(Langsung | Olahraga) = \frac{40/100}{50/100} = \frac{0.4}{0.5} = 0.8\)
• \(P(Tunda| Olahraga) = \frac{10/100}{50/100} = \frac{0.1}{0.5} = 0.2\)
• \(P(Langsung | Tidak Olahraga) = \frac{20/100}{50/100} = \frac{0.2}{0.5} = 0.4\)
• \(P(Tunda| Tidak Olahraga) = \frac{30/100}{50/100} = \frac{0.3}{0.5} = 0.6\)

Pada peluang bersyarat dapat diambil interpretasi yakni \(P(Langsung | Olahraga) = 0.8\)

\(P(Tunda | Tidak Olahraga) = 0.4\)

Menunjukkan bahwa orang yang sering berolahraga lebih cenderung langsung bangun daripada orang yang jarang berolahraga.

2.3 Ukuran Asosiasi Pada Data Kategori 2 x 2

Hubungan antar variabel acak X dan Y yang keduanya bersifat kategori dalam tabel kontingensi dinamakan dengan asosiasi. Asosiasi ini dibentuk dari distribusi peluang bersyaratnya. Ada 3 jenis ukuran asosiasi pada tabel kontingensi 2 arah yakni

1. Beda Peluang (Risk Difference)

   Ukuran absolut yang menggambarkan beda antara probabilitas dari dua kelompok berbeda. Rumusnya sebagai berikut

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

2. Resiko Relatif (Relative Risk)

   Ukuran yang digunakan untuk membandingkan resiko kejadian suatu peristiwa antara dua kelompok, yakni kejadian A dan bukan A.

   Rumus: \[ RR = \frac{ \left( \frac{n_{11}}{n_{1.}} \right) }{ \left( \frac{n_{21}}{n_{2.}} \right) } \]

3. Rasio Odds (Odds ratio)

   Ukuran yang digunakan untuk membandingkan peluang terjadinya kejadian dua kelompok yakni kelompok kejadian A dan bukan A. Bedanya sama resiko relatif, odd ratio ini membandingkan peluangnya sedangkan RR hanya membandingkan resikonya.

   Rumus: \[ OR = \frac{ \left( \frac{n_{11}}{n_{12}} \right) }{ \left( \frac{n_{21}}{n_{22}} \right) } = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} \]

2.3.1 Contoh Perhitungan

Pada soal sebelumnya telah diperoleh beberapa nilai, sehingga dapat dilanjutkan perhitungan ukuran asosiasinya.

1. Beda Peluang

   Beda peluang mengukur selisih risiko absolut antara dua kelompok.

   Rumus: \[ RD = P(Langsung | Olahraga) - P(Tunda | Tidak Olahraga) \]

   Hitung:

   \(RD = 0.8 - 0.4 = 0.4\)

   Orang yang berolahraga memiliki peluang 40% lebih besar untuk langsung bangun dari tidur dibandingkan dengan yang tidak berolahraga.

2. Resiko Relatif

   Mengukur perbandingan risiko antara dua kelompok.

   Rumus: \[ RR = \frac{P(Langsung | Olahraga)}{P(Tunda| Tidak Olahraga)} \]

   Hitung:

   \(RR = \frac{0.8}{0.4} = 2\)

   Orang yang berolahraga memiliki 2 kali lipat kemungkinan untuk dapat langsung bangun dari tidur dibandingkan yang tidak berolahraga.

3. Odd Ratio

   Mengukur perbandingan odds antara dua kelompok.

   Rumus: \[ OR = \frac{(40/10)}{(20/30)} = \frac{4}{0.6667} \approx 6 \]

   Peluang untuk orang dapat langsung bangun dari tidur pada orang yang berolahraga adalah 6 kali lebih tinggi dibandingkan yang tidak berolahraga.

   Contoh Perhitungan R :

   RD <- function(n11, n12, n21, n22) {
     (n11 / (n11 + n12))- (n21 / (n21 + n22))
   }
   RD(40, 10, 20, 30)

   ## [1] 0.4

   RR <- function(n11, n12, n21, n22) {
     (n11 / (n11 + n12)) / (n21 / (n21 + n22))
   }
   RR(40, 10, 20, 30)

   ## [1] 2

   OR <- function(n11, n12, n21, n22) {
     (n11 * n22) / (n12 * n21)
   }
   OR(40, 10, 20, 30)

   ## [1] 6

2.4 Inferensial dalam Tabel Kontingensi Dua Arah

Inferensial dalam konteks tabel kontingensi dua arah mengacu pada proses menarik kesimpulan tentang hubungan antara dua variabel kategori dalam populasi, berdasarkan data sampel yang disajikan dalam tabel tersebut. Dalam hal ini, analisis tidak hanya berhenti pada pengamatan data yang terlihat di tabel, tetapi juga melibatkan pertanyaan lebih lanjut seperti: “Apakah pola ini terjadi karena kebetulan dalam sampel, atau mencerminkan pola nyata di populasi?” Dengan kata lain, pendekatan inferensial bertujuan untuk menguji apakah ada asosiasi yang bermakna antara dua variabel kategori.

Melalui analisis inferensial, kita dapat mengembangkan pemahaman yang lebih dalam mengenai bagaimana dua variabel saling berhubungan, apakah hubungan itu bersifat lemah atau kuat, atau bahkan tidak ada sama sekali. Ini penting karena dalam banyak situasi, perbedaan yang tampak mencolok di tabel belum tentu berarti signifikan secara statistik. Oleh karena itu, pendekatan inferensial membantu membedakan antara pola yang kebetulan muncul dalam sampel dengan pola yang kemungkinan besar berlaku di seluruh populasi.

2.4.1 Estimasi

1. Titik, digunakan untuk menentukan nilai spesifik sebagai perkiraan terbaik parameter.

\[ \hat{p}_1 = \frac{n_{11}}{n_{1.}} \]

2. Interval, digunakan untuk memberikan rentang nilai yang diyakini memiliki parameter populasi

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

2.4.2 Uji Hipotesis

2.4.2.1 Uji Proporsi

Uji ini digunakan untuk membandingkan proporsi kejadian antara dua kelompok dalam tabel kontingensi terutama untuk menentukan apakah terdapat perbedaan signifikan dalam proporsi kejadian.

Untuk menguji hal tersebut, diawali dengan menyatakan hipotesis yakni

H0 : Tidak ada perbedaan antara dua kelompok, yaitu \[ p_1 = p_2 \]

H1 : Terdepat perbedaan proporsi antara dua kelompok.

Untuk formulanya sebagai berikut

1. estimasi proporsi masing masing

   \[ \hat{p}_1 = \frac{n_{11}}{n_{1.}} = \frac{40}{50} = 0.8 \]

   dan

   \[ \hat{p}_2 = \frac{n_{21}}{n_{2.}} = \frac{20}{50} = 0.4 \]

2. estimasi proporsi gabungan

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

3. statistik uji Z

   Statsitik uji Z mengikuti distribusi normal baku (0,1) dan p value dihitung berdasarkan nilai kritis.

\[ Z = \frac{\hat{p}_1 + \hat{p}_2}{\sqrt{\hat{p}(1 - \hat{p})\left( \frac{1}{n_1} + \frac{1}{n_2} \right)}} \\ = \frac{0.8 + 0.4}{\sqrt{0.6(1 - 0.6)\left( \frac{1}{50} + \frac{1}{50} \right)}} \\ = 15.042 \]

4. Evaluasi

   Jika nilai Z lebih dari nilai kritis atau p-value kurang dari alpha yang ditetapkan, maka dapat diyakini bahwa ada perbedaan proporsi dari kedua kelompok tersebut.

Contoh perhitungan menggunakan R

# Matriks data: baris = Olahraga, kolom = Bangun
data <- matrix(c(40, 10, 20, 30), nrow = 2, byrow = TRUE)
dimnames(data) <- list("Olahraga" = c("Ya", "Tidak"),
                       "Sehat" = c("Langsung", "Tunda"))
data

##         Sehat
## Olahraga Langsung Tunda
##    Ya          40    10
##    Tidak       20    30

# Uji Proporsi dengan variabel yang benar
 prop_test <- prop.test(x =c(data[1,1], data[2,1]),
 n =c(sum(data[1,]), sum(data[2,])))
 print(prop_test)

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(data[1, 1], data[2, 1]) out of c(sum(data[1, ]), sum(data[2, ]))
## X-squared = 15.042, df = 1, p-value = 0.0001052
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.2046955 0.5953045
## sample estimates:
## prop 1 prop 2 
##    0.8    0.4

Karena nilai p-value kurang dari 5% sesuai taraf yang ditetapkan, dapat dipercaya bahwa proporsi kejadian antar kebiasaan olahraga dan cara bangun tidur berbeda signifikan.

2.4.2.2 Uji Asosiasi

Uji asosiasi ini bertujuan mengukur hubungan antar dua variabel kategori. Dengan ini disusun hipotesis sebagai berikut

H0 : Tidak ada asosiasi antara dua variabel

H1 : Terdapat asosiasi antara dua variabel

Adapun formulanya sebagai berikut

1. Beda Peluang Rumus RD

   \[ RD = \hat{p}_1 - \hat{p}_2 \]

   Rumus Standar Error RD

   \[ SE(RD) = \sqrt{ \frac{\hat{p}_1 (1 - \hat{p}_1)}{n_1} + \frac{\hat{p}_2 (1 - \hat{p}_2)}{n_2} } \]

   Rumus Statistik Uji Z \[ Z = \frac{RD}{SE(RD)} = \frac{\hat{p}_1 - \hat{p}_2}{\sqrt{ \frac{\hat{p}_1 (1 - \hat{p}_1)}{n_1} + \frac{\hat{p}_2 (1 - \hat{p}_2)}{n_2} }} \]

2. Resiko Relatif

   Rumus :

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

   Standar Error

   \[ SE(\ln RR) = \sqrt{ \frac{1 - \hat{p}_1}{x_1} + \frac{1 - \hat{p}_2}{x_2} } \]

   Statistik Z

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

3. Odd Ratio

   Rumus : \[ OR = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} \]Standard Error dari log(OR): \[ SE(\ln OR) = \sqrt{ \frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} } \]

   Z Statistik: \[ Z = \frac{\ln(OR)}{SE(\ln OR)} \]

Contoh Perhitungan

Pada kasus sebelumnya dijelaskan ada 2 kondisi yakni kebiasaan olahraga dan cara bangun tidur. Untuk perhitungan contoh dilakukan pada interval odd ratio sebagai berikut.
Tadi diperoleh bahwa nilai odd rationya adalah = 6

Maka sesuaikan rumusnya

\(SE(\ln OR) = \sqrt{ \frac{1}{40} + \frac{1}{10} + \frac{1}{20} + \frac{1}{30}} = 0.4564\)

Maka diperoleh nilai Z

\[ Z = \frac{(OR)}{SE(\ln(OR)} = \frac{6}{0.4565} = 3.9255 \]

n11 <- 40 
n12 <- 10 
n21 <- 20 
n22 <- 30
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.4
## 
## $SE_RD
## [1] 0.08944272
## 
## $Z_RD
## [1] 4.472136
## 
## $RR
## [1] 2
## 
## $SE_Ln_RR
## [1] 0.1870829
## 
## $Z_RR
## [1] 3.705028
## 
## $OR
## [1] 6
## 
## $SE_Ln_OR
## [1] 0.4564355
## 
## $Z_OR
## [1] 3.925548

Dengan nilai Z setiap perhitungan dapat dilakukan penilaian terhadap signifikansi masing masing ukuran asosiasi. Jika menggunakan alpha = 5% maka nilai kritisnya adalah 1,96 yang berarti ketiga ukuran asosiasi tersebut memberikan perbedaan yang signifikan karena nilai Z melebihi nilai kritisnya.

2.4.2.3 Uji Independensi

Uji ini digunakan untuk menentukan apakah ada hubungan statistik antara dua variabel kategorikal. Uji yang digunakan salah satunya adlaah uji Chi-Square dengan rumus

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

dimana :

O = nilai observasi dalam tabel kontingensi

E = nilai yang diharapkan dan diperoleh dari rumus

\[ E_{ij} = \frac{R_i \cdot C_j}{N} \]

di mana: - \(R_i\): jumlah total baris ke-\(i\) - \(C_j\): jumlah total kolom ke-\(j\) - \(N\): total seluruh observasi

Sama halnya dengan uji Z, uji chi square juga memiliki ketentuan hipotesis sebagai berikut

H0 : Tidak terdapat hubungan antar kelompok

H1 : Terdapat hubungan signifikan antar kelompok

Dengan ketentuan jika nilai chi-square lebih dari nilai kritis atau p-value kurang dari alpha, maka H0 ditolak.

Contoh Perhitungan : Pada kasus sebelumnya, kita memiliki data, oleh karena itu kita akan melihat apakah ada hubungan antara kebiasaan olahraga dengan cara bangun tidur.

Secara Manual : Misalkan kita memiliki tabel berikut:

	Kejadian Ya	Kejadian Tidak	Total
Terpapar	40	10	50
Tidak Terpapar	20	30	50
Total	50	50	100

Hitung nilai yang diharapkan \(E\):

\[ E_{11} = \frac{50 \times 50}{100} = 25 \]

\[ E_{12} = \frac{50 \times 50}{100} = 25 \]

\[ E_{21} = \frac{50 \times 50}{100} = 25 \]

\[ E_{22} = \frac{50 \times 50}{100} = 25 \]

Hitung nilai \(\chi^2\):

\[ \chi^2 = \frac{(40 - 25)^2}{25} + \frac{(10 - 25)^2}{25} + \frac{(20 - 25)^2}{25} + \frac{(30 - 25)^2}{25} \]

\[ = 15.042 \] Hitung Derajat Kebebasan

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

Jika \(p\)-value < 0.05, maka kita menolak hipotesis nol.

Dengan ini diperoleh bahwa kita menolak H0 yang menandakan kalau ada hubungan antara kebiasaan olahraga dengan cara bangun tidur.

Perhitungan R :

#Load data
data <- matrix(c(40, 10, 20, 30), nrow = 2, byrow = TRUE)
dimnames(data) <- list("Olahraga" = c("Ya", "Tidak"), "Cara Bangun" = c("Ya", "Tidak"))
print(data)

##         Cara Bangun
## Olahraga Ya Tidak
##    Ya    40    10
##    Tidak 20    30

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

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

Dengan perhitungan tersebut diperoleh nilai p-value yang kurang dari alpha (0,05). Oleh karena itu dapat disimpulkan bahwa terdapat hubungan antara kebiasaan olahraga dengan cara bangun tidur

Selain dengan chi-square dapat dihitung juga dengan metode lain seperti

1. Partisi Chi-Square

   metode untuk membagi total nilai chi-square menjadi beberapa bagian sesuai dengan perbandingan yang lebih spesifik antar kategori. Dengan membandingkan setiap bagian, peneliti dapat mengidentifikasi sumber ketidaksesuaian dalam data. Ini memungkinkan analisis yang lebih terfokus dibanding hanya melihat uji chi-square secara keseluruhan.

   Contoh kasus : Dilakukan sebuah survey dan diperoleh data dengan kriteria sebagai berikut

   Gender	Democrat	Republican	Independent	Total
   Female	495	272	590	1357
   Male	330	265	498	1093
   Total	825	537	1088	2450

   Maka diminta apakah ada hubungan antara gender dengan identifikasi partai politik?

   Langkah pertama yakni dengan menghitungan chi-square secara keseluruhan

   # Data Observasi
   data_matrix <- matrix(c(495, 272, 590, 330, 265, 498), nrow = 2, byrow = TRUE)
   colnames(data_matrix) <- c("Democrat", "Republican", "Independent")
   rownames(data_matrix) <- c("Female", "Male")
   
   # Uji Chi-Square
   chi_test <- chisq.test(data_matrix)
   
   # Hasil
   list(
     Chi_Square = chi_test$statistic,
     P_Value = chi_test$p.value,
     Decision = ifelse(chi_test$p.value < 0.05, "Tolak H0 (Ada hubungan)", "Gagal Tolak H0 (Tidak ada hubungan)")
   )

   ## $Chi_Square
   ## X-squared 
   ##  12.56926 
   ## 
   ## $P_Value
   ## [1] 0.00186475
   ## 
   ## $Decision
   ## [1] "Tolak H0 (Ada hubungan)"

   Langkah selanjutnya adalah dengan melakukan partisi chi–square democrat vs republician dan partisi democrat + republic vs independet. Terakhir tinggal di bandingkan dengan nilai chi square tabel

   # Data Observasi
   data_matrix <- matrix(c(495, 272, 330, 265), nrow = 2, byrow = TRUE)
   colnames(data_matrix) <- c("Democrat", "Republican")
   rownames(data_matrix) <- c("Female", "Male")
   # Uji Chi-Square Partisi 1
   chi_test1 <- chisq.test(data_matrix)
   # Data Partisi 2
   data_matrix2 <- matrix(c(767, 590, 595, 498), nrow = 2, byrow = TRUE)
   colnames(data_matrix2) <- c("Dem+Rep", "Independent")
   rownames(data_matrix2) <- c("Female", "Male")
   # Uji Chi-Square Partisi 2
   chi_test2 <- chisq.test(data_matrix2)
   # Hasil
   list(Chi_Square_Partisi1 = chi_test1, Chi_Square_Partisi2 = chi_test2)

   ## $Chi_Square_Partisi1
   ## 
   ##  Pearson's Chi-squared test with Yates' continuity correction
   ## 
   ## data:  data_matrix
   ## X-squared = 11.178, df = 1, p-value = 0.0008279
   ## 
   ## 
   ## $Chi_Square_Partisi2
   ## 
   ##  Pearson's Chi-squared test with Yates' continuity correction
   ## 
   ## data:  data_matrix2
   ## X-squared = 0.98267, df = 1, p-value = 0.3215

   Dengan ini diperoleh bahwa gender mempengaruhi pemilihan democrat dan republic. Hal ini ditunjukkan jelas dengan nilai pvalue yang sangat kecil. Sedangkan pilihan independent tidak berpengaruh signifikan berdasarkan gender.

2. Uji Likelihood Ratio

   Ini merupakan alternatif dari pengujian chisquare.

   \[ G^2 = 2 \sum_{i=1}^{k} O_i \ln\left(\frac{O_i}{E_i}\right) \]

   Dimana:

   - \(O_i\) = Frekuensi yang diamati pada kategori \(i\)

   - \(E_i\) = Frekuensi yang diharapkan pada kategori \(i\) berdasarkan distribusi yang dihipotesiskan,

   - \(k\) = Jumlah kategori atau kelompok yang ada dalam data,

   - \(\ln\) = Fungsi logaritma natural.

   Perhitungannya manual :

   Langkah 1: Hitung Frekuensi Ekspektasi \[\hat{\mu}_{11} = \frac{50 \times 50}{100} = 25\]
   \[\hat{\mu}_{12} = \frac{50 \times 50}{100} = 25\]
   \[\hat{\mu}_{21} = \frac{50 \times 50}{100} = 25\]
   \[\hat{\mu}_{22} = \frac{50 \times 50}{100} = 25\] Langkah 2 : Hitung Statistik Uji G \[ G^2 = 2 \left[ 40 \ln\left(\frac{40}{25}\right) + 10 \ln\left(\frac{10}{25}\right) + 20 \ln\left(\frac{20}{25}\right) + 30 \ln\left(\frac{30}{25}\right) \right] \] \[ = 17.26 \]

   Langkah 3 : Bandingkan dengan Distribusi Chi Square

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

\[ \chi^2_{0.05, 1} = 3.841 \]

Karena \[G^2 = 17.26 > 3.841\] maka tolak \[H_0\] dan simpulkan bahwa ada hubungan signifikan antara status merokok dan kanker paru-paru. Demikian juga dengan perhitungan menggunakan R

# Data Observasi
data <- matrix(c(40, 10, 20, 30), nrow = 2, byrow = TRUE)
colnames(data) <- c("Langsung", "Tunda")
rownames(data) <- c("Olahraga", "Jarang")

# Hitung Frekuensi Ekspektasi
data_expected <- chisq.test(data)$expected

# Hitung Statistik G²
G2 <- 2 * sum(data * log(data / 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 to Reject H0")
)

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

3. Uji Exact Fisher

   Digunakan untuk menguji hubungan antara dua variabel kategori dalam tabel kontingensi kecil, terutama ukuran 2x2, ketika frekuensi sel rendah. Uji ini memberikan nilai p yang akurat tanpa mengandalkan distribusi chi-square. \[ P = \frac{ \binom{a + b}{a} \binom{c + d}{c} }{ \binom{n}{a + c} } \] Cocok untuk sampel kecil dan ketika asumsi uji chi-square tidak dipenuhi.

   data <- matrix(c(40, 10, 20, 30), nrow = 2, byrow = TRUE)
   fisher.test(data)

   ## 
   ##  Fisher's Exact Test for Count Data
   ## 
   ## data:  data
   ## p-value = 8.309e-05
   ## alternative hypothesis: true odds ratio is not equal to 1
   ## 95 percent confidence interval:
   ##   2.265787 16.390918
   ## sample estimates:
   ## odds ratio 
   ##   5.881258

2.4.2.4 Analisis Residual

Dalam tabel kontingensi, analisis residual digunakan untuk mengidentifikasi sel yang menyumbang paling banyak terhadap hubungan antara variabel. Jika nilai residual 0 maka variabel baris dan kolom tidak punya hubungan kuat atau polanya tidak jelas. Dalam pengerjaannnya ada 2 macam analisis untuk residual, yakni

1. Pearson Residual

   Selisih antara frekuensi yang diobservasi (O) dan diharapkan (E) dalam sebuah sel, yang kemudian dibagi dengan akar dari nilai yang diharapkan. Rumusnya adalah:

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

   Nilai residual ini menunjukkan seberapa besar deviasi antara data aktual dengan ekspektasi berdasarkan hipotesis nol. Semakin besar nilainya (positif atau negatif), semakin besar kontribusi sel tersebut terhadap nilai total chi-square, dan semakin “tidak sesuai” sel tersebut dengan model nol.

2. Standardized Residual

   Standardized Residual adalah Pearson Residual yang telah disesuaikan dengan mempertimbangkan varians dari distribusi. Dalam konteks tabel kontingensi, residual ini dihitung untuk memperbaiki bias dan ketidakseimbangan dalam data. Rumus umumnya adalah:

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

   Residual ini memberikan ukuran yang lebih tepat tentang keanehan suatu sel, dengan memperhitungkan posisi sel dalam struktur tabel. Nilai residual yang lebih besar dari ±2 sering dianggap menunjukkan bahwa sel tersebut menyimpang secara signifikan dari yang diharapkan.

Contoh perhitungannya

# Data Observasi
observed <- matrix(c(40, 10, 20, 30), nrow = 2, byrow = TRUE)

# Hitung nilai ekspektasi
expected <- chisq.test(observed)$expected
 
# Pearson Residual
pearson_residual <- (observed- expected) / sqrt(expected)
 
# Standardized Residual
 row_sum <- rowSums(observed)
 col_sum <- colSums(observed)
 total_sum <- sum(observed)

 # Menampilkan hasil
standardized_residual <- (observed- expected) / sqrt(expected * (1- row_sum/total_sum)*(1-col_sum/total_sum))
 
list(
    Pearson_Residual = pearson_residual,
    Standardized_Residual = standardized_residual
)

## $Pearson_Residual
##           [,1]      [,2]
## [1,]  1.825742 -2.236068
## [2,] -1.825742  2.236068
## 
## $Standardized_Residual
##           [,1]      [,2]
## [1,]  4.082483 -5.000000
## [2,] -3.333333  4.082483

Kemudian dapat juga mendeteksi outlier. Outlier sendiri adalah data yang jauh dari penyebaran data lainnya. Contoh perhitungannya sebagai berikut

# Data Observasi
observed <- matrix(c(40, 10, 20, 30), nrow = 2, byrow = TRUE)
colnames(observed) <- c("Sukses", "Gagal")
rownames(observed) <- c("Grup A", "Grup B")

# Hitung nilai ekspektasi
expected <- chisq.test(observed)$expected
 
# Pearson Residual
pearson_residual <- (observed- expected) / sqrt(expected)

# Standardized Residual
row_sum <- rowSums(observed)
col_sum <- colSums(observed)
total_sum <- sum(observed) 
standardized_residual <- (observed- expected) / sqrt(expected * (1- row_sum/total_sum)*(1-col_sum/total_sum))

# Menampilkan hasil
list(
 Observed = observed,
 Expected = expected,
 Pearson_Residual = pearson_residual,
 Standardized_Residual = standardized_residual
)

## $Observed
##        Sukses Gagal
## Grup A     40    10
## Grup B     20    30
## 
## $Expected
##        Sukses Gagal
## Grup A     30    20
## Grup B     30    20
## 
## $Pearson_Residual
##           Sukses     Gagal
## Grup A  1.825742 -2.236068
## Grup B -1.825742  2.236068
## 
## $Standardized_Residual
##           Sukses     Gagal
## Grup A  4.082483 -5.000000
## Grup B -3.333333  4.082483

Jika nilai | r | > 3, maka data tersebut dianggap sebagai outlier. Pada data ini dapat dilihat bahwa terdapat outlier jika diuji menggunakan standardized residual, namun dengan pearson tidak terdapat outlier yang signifikan.

2.5 Review

Pada Bab 2 ini dibahas beberapa analisis terhadap tabel kontingensi 2 arah.

1. Tabel Kontingensi dua arah merupakan tabel silang yang berisi data frekuensi dari treatment. Metode ini adalah analisis yang sederhana dari data kategori. Dari tabel ini dapat dicari beberapa peluang diantaranya peluang bersama (joint probabilities), peluang marginal dan peluang bersyarat.
2. Pada Tabel Kontingensi ini dapat ditentukan ukuran asosiasi untuk menentukan hubungan antar variabel X dan Y. Ukuran asosiasi ini terdiri dari Beda Peluang, Resiko Relatif dan Odd Rasio.
3. Inferensi pada tabel kkotingensi 2 arah bertujuan utnuk mengevaluasi hubungan variabel dalam tabel kontingensi. Terdiri banyak metode untuk estimasi point dan intervalnya
4. Uji Proporsi digunakan untuk membandingkan proporsi kejadian antar kedua kejadian
5. Uji Asosiasi dilakukan untuk mengukur dan mengevaluasi hubungan antar variabel
6. Uji Independensi dilakukan utnuk menentukan apakah ada hubungan statistikal antara kedua variabel ketagorikal. Pengujiannya dapat dilakukan dengan chi square, uji likelihood rasio, exact fisher dan lainnya.
7. Analisis Residual bertujuan untuk mengidentifikan sel yang menyumbangkan hubungan paling banyak serta melihat outlier.

3 Tabel Kontingensi Tiga Arah

3.1 Definisi Umum

Merupakan tabel yang digunakan muntuk meyajikan dan menganalisis 3 variabel kategori secara langsung. Dalam banyak situasi, biasanya hubungan dari variabel X dan Y dapat diperngaruhi oleh variabel lainnya, misalkan Z. Oleh karena itu variabel Z ini disebut sebagai variabel kontrol atau kovariat.

Contohnya

Hubungan antara kebiasaan begadang dengan seringnya lupa pada mahasiswa dan pekerja.

Kebiasaan orang-orang sekarang adalah bekerja ketika mepet dengan deadlinenya. Hal ini jadi salah satu pemicu utama banyak orang-orang menjadikan begadang menjadi kebiasaan. Padahal sebenarnya kebiasaan begadang merupakan kebiasaan buruk yang mungkin saja memberikan pengaruh yang buruk. Kebiasaan ini telah menjalar ke semua kelompok manusia, terutama mahasiswa dan orang yang bekerja.

Mahasiswa sering begadang untuk menyelesaikan tuntutan akademik seperti mengerjakan tugas atau sibuk dengan kegiatan keorganisasian atau bahkan sibuk dengan hiburan hingga lupa istirahat. Sementara itu, pekerja begadang karena beban kerja tinggi tekanan dan lainnya. Hal ini bisa saja menjadi salah satu alasan turunya daya ingat dan konsentrasi.

Menurut Walker (2017), kebiasaan begadang merupakan kebiasaan yang membawa efek buruk, karena tidur memberikan fungsi yang baik untuk otak, termasuk memori dan konsentrasi. Selain itu, dikuatkan oleh Curcio (2006) menyatakan bahwa kurang tidur dapat memengauhi kemampuan belajar dan daya ingat.

Untuk itu, pada kasus ini akan dianalisis dari 120 orang responden untuk mengetahui apakah kebiasaan begadang secara signifikan memengauhi seseorang menjadi lebih sering lupa. Variabel yang digunakan adalah sebagai berikut

X = Kebiasaan begadang (Ya atau Tidak)

Y = Sering Lupa (Ya atau Tidak)

Z = Status (Mahasiswa atau Pekerja)

Untuk hasil survei dapat diperoleh sebagai berikut

Status	Begadang	Sering Lupa YA	Sering Lupa Tidak
Mahasiswa	Ya	30	10
Mahasiswa	Tidak	8	12
Pekerja	Ya	18	7
Pekerja	Tidak	10	25

Pada data ini, variabel status bekerja sebagai variabel confounder. Confounder adalah variabel yang bisa memengaruhi baik variabel independen (X) maupun variabel dependen (Y). Hal ini terkadang bisa memberikan hubungan antar variabel dependen dan independent padahal sebenarnya hubungan terjadi karena confounder yang bukan sebab utama dari kasus tersebut.

Pada kasus tadi, variabel pekerja dan mahasiswa ini dapat menjadi confouder karena dapat memengaruhi kebiasaan begadang maupun sering lupa. Jika hal ini tidak dikontrol, hubungan antar keduanya bisa tampak lebih kuat maupun lebih lemah. oleh karena itu perlu adanya kontrol terhadap variabel Z ini agar dapat memastikan bahwa kebiasaan begadang memang menjadi faktor utama dari sering lupa.

3.2 Tabel Parsial dan Marginal

Tabel parsial merupakan tabel yang menyajikan hubungan natara X dan Y pada setiap kategori Z. INi memungkinkan analisis hubungan bersyarat dimana efek Z dikendalikan.

Tabel Marginal adalah tabel yang diperoleh dengan mengabaikan Z yaitu dengan menjumlahkan semua kateogir Z. Ini memberikan gambaran umum antara X dan Y tanpa mempertimbangkan Z. Hal ini yang dapat menyebabkan distorsi interpretasi seperti dalan Simpson’ Paradox.

Simpson paradox adalah enomena dalam analisis data dimana tren atau hubungan yang terlihat pada kelompok kecil berbalik arah atau hilang ketika data digabungkan secara keseluruhan.

Oleh karena itu, analisis yang mempertimbangkan tabel parsial lebih diutamakan untuk mendapatkan kesimpulan yang lebih akurat dalam penelitian Ilmiah.

Pada kasus sebelumya, kita dapat membentuk tabel frekuensi parsial dan marginalnya sebagai berikut

1. Tabel Parsial

   Tabel ini merupakan tabel yang menyajikan data dari dua variabel dari tabel kontingensi tiga arah dengan satu variabel lainnya sebagai kontrol.

   # Instal Package
   library(knitr)
   # Data awal
   data <- array(
     c(30,10,8,12,18,7,10,25),
     dim = c(2,2,2),
     dimnames = list(
       Y = c("Sering Lupa", "Tidak Sering Lupa"),
       X = c("Begadang", "Tidak Begadang"),
       Z = c("Mahasiswa", "Bekerja")
     )
   )
   
   data

   ## , , Z = Mahasiswa
   ## 
   ##                    X
   ## Y                   Begadang Tidak Begadang
   ##   Sering Lupa             30              8
   ##   Tidak Sering Lupa       10             12
   ## 
   ## , , Z = Bekerja
   ## 
   ##                    X
   ## Y                   Begadang Tidak Begadang
   ##   Sering Lupa             18             10
   ##   Tidak Sering Lupa        7             25

2. Tabel Marginal

   Tabel yang menampilkan jumlah total observasi untuk setiap variabel dengan mengabaikan variabel lainnya dalam tabel kontingensi tiga arah.

   margin_YX <- apply(data, c(1,2),sum)
   margin_YZ <- apply(data, c(1,3),sum)
   margin_XZ <- apply(data, c(2,3),sum)
   
   cat("\n ##Marginal Y dan X \n")

   ## 
   ##  ##Marginal Y dan X

   kable(data, format = "html", booktabs = TRUE)

   	Begadang.Mahasiswa	Tidak Begadang.Mahasiswa	Begadang.Bekerja	Tidak Begadang.Bekerja
   Sering Lupa	30	8	18	10
   Tidak Sering Lupa	10	12	7	25

   kable(margin_YX, format = "html", booktabs = TRUE)

   	Begadang	Tidak Begadang
   Sering Lupa	48	18
   Tidak Sering Lupa	17	37

3. Tabel Peluang Bersama, Marginal dan Bersyarat

   Tabel peluang bersama didefinisikan sebagai tabel yang berisi peluang bersama dari setiap sel. Tabel peluang Marginal merupakan sum dari peluang bersama sedangkan tabel peluang bersyarat merupakan tabel yang dihitung berdasarkan jumlah total dalam setiap kategori. Contohnya :

   data3 <- array(c(30, 10, 8, 12, 18, 7, 10, 25),
     dim = c(2, 2, 2),
     dimnames = list(
     Begadang = c("Ya", "Tidak"),
     Pelupa = c("Ya", "Tidak"),
     Status = c("Mahasiswa", "Pekerja")
   ))
   # Hitung probabilitas bersama
   total <- sum(data3)
   joint_prob <- data3 / total
   ftable(joint_prob)

   ##                 Status  Mahasiswa    Pekerja
   ## Begadang Pelupa                             
   ## Ya       Ya            0.25000000 0.15000000
   ##          Tidak         0.06666667 0.08333333
   ## Tidak    Ya            0.08333333 0.05833333
   ##          Tidak         0.10000000 0.20833333

   # Hitung probabilitas marginal
   marginal_X <- apply(joint_prob, 1, sum)
   marginal_Z <- apply(joint_prob, 3, sum)
   # Tampilkan hasil
   marginal_X

   ##    Ya Tidak 
   ##  0.55  0.45

3.3 Ukuran Asosiasi

Ukuran asosiasi adalah ukuran yang digunakan untuk mengukur hubungan antara dua variabel, terutama dalam data kategorik. Ada beberapa yang akan diujikan yakni beda peluang, resiko relatif dan odd ratio. Sama halnya dengan tabel kontingensi dua arah, perhitungan ukuran asosiasi ini menggunakan formula yang sama dengan tabel peluang yang telah diberikan sebelumnya.

3.3.1 Menghitung Ukuran Asosiasi Mahasiswa

Ukuran asosiasi digunakan untuk melihat interpretasi dari tabel mahasiswa.

a <-data["Sering Lupa", "Begadang", "Mahasiswa"]
b <-data["Tidak Sering Lupa", "Begadang", "Mahasiswa"]
c <-data["Sering Lupa", "Tidak Begadang", "Mahasiswa"]
d <-data["Tidak Sering Lupa", "Tidak Begadang", "Mahasiswa"]

p1 <- a / (a + b)
p2 <- c / (c + d)

RD_mhs <- p1 - p2
RR_mhs <- p1 / p2
OR_mhs <- (a*d)/(b*c)

ukuran_asosiasi_mhs <- data.frame(
  Statistik = c("RD", "RR", "OR"),
  Nilai = c(RD_mhs, RR_mhs, OR_mhs)
)

kable(ukuran_asosiasi_mhs, format = "html", booktabs = TRUE, col.names = c("Ukuran Asosiasi", "Nilai"))

Ukuran Asosiasi	Nilai
RD	0.350
RR	1.875
OR	4.500

Dengan ini diperoleh bahwa

1. Mahasiswa yang biasa begadang memiliki peluang lupa lebih besar sebanyak 35% daripada mahasiswa yang tidak biasa begadang.
2. Mahasiswa yang biasa begadang memiliki 1,8x kemungkinan lebih sering lupa daripada mahasiswa yang tidak biasa begadang.
3. Mahasiswa yang biasa begadang memiliki kemungkinan 4,5 kali lebih sering lupa daripada mahasiswa yang tidak biasa begadang.

3.3.2 Ukuran Asosiasi Untuk Orang Bekerja

a <-data["Sering Lupa", "Begadang", "Bekerja"]
b <-data["Tidak Sering Lupa", "Begadang", "Bekerja"]
c <-data["Sering Lupa", "Tidak Begadang", "Bekerja"]
d <-data["Tidak Sering Lupa", "Tidak Begadang", "Bekerja"]

p1 <- a / (a + b)
p2 <- c / (c + d)

RD_bkr <- p1 - p2
RR_bkr <- p1 / p2
OR_bkr <- (a*d)/(b*c)

ukuran_asosiasi_bkr <- data.frame(
  Statistik = c("RD", "RR", "OR"),
  Nilai = c(RD_bkr, RR_bkr, OR_bkr)
)

kable(ukuran_asosiasi_bkr, format = "html", booktabs = TRUE, col.names = c("Ukuran Asosiasi", "Nilai"))

Ukuran Asosiasi	Nilai
RD	0.4342857
RR	2.5200000
OR	6.4285714

Dengan ini diperoleh bahwa

1. Pekerja yang biasa begadang memiliki peluang lupa lebih besar sebanyak 43% daripada pekerja yang tidak biasa begadang.
2. Pekerja yang biasa begadang memiliki 2,52x kemungkinan lebih sering lupa daripada pekerja yang tidak biasa begadang.
3. Pekerja yang biasa begadang memiliki kemungkinan 6,42 kali lebih sering lupa daripada Pekerja yang tidak biasa begadang.

3.3.3 Ukuran Asosiasi untuk tabel marginal

a <-margin_YX["Sering Lupa", "Begadang"]
b <-margin_YX["Tidak Sering Lupa", "Begadang"]
c <-margin_YX["Sering Lupa", "Tidak Begadang"]
d <-margin_YX["Tidak Sering Lupa", "Tidak Begadang"]

p1 <- a / (a + b)
p2 <- c / (c + d)

RD_marginal <- p1 - p2
RR_marginal <- p1 / p2
OR_marginal <- (a*d)/(b*c)

ukuran_asosiasi_marginal <- data.frame(
  Statistik = c("RD", "RR", "OR"),
  Nilai = c(RD_marginal, RR_marginal, OR_marginal)
)

kable(ukuran_asosiasi_marginal, format = "html", booktabs = TRUE, col.names = c("Ukuran Asosiasi", "Nilai"))

Ukuran Asosiasi	Nilai
RD	0.4111888
RR	2.2564103
OR	5.8039216

Dengan ini diperoleh bahwa

1. Orang yang biasa begadang memiliki peluang lupa lebih besar sebanyak 41% daripada orang yang tidak biasa begadang.
2. Orang yang biasa begadang memiliki 2,25x kemungkinan lebih sering lupa daripada orang yang tidak biasa begadang.
3. Orang yang biasa begadang memiliki kemungkinan 5,8 kali lebih sering lupa daripada orang yang tidak biasa begadang.

3.3.4 Conditional Independent

Mengacu pada situasi di mana dua variabel bebas secara statistik jika diketahui nilai dari variabel ketiga. Dalam kata lain, dua variabel tampak memiliki hubungan, tetapi hubungan itu sebenarnya dijelaskan sepenuhnya oleh variabel ketiga.

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

\[ \frac{n_{ijk}}{n_{k++}} = \left( \frac{n_{i+k}}{n_{k++}} \right) \cdot \left( \frac{n_{+jk}}{n_{k++}} \right) \]

Selain itu, dilakukan pengujian terhadap conditional Independent untuk membuktikan apakah sebuah variabel memiliki hubungan signifikan setelah dikendalikan oleh variabel kontrol. Contoh perhitungannya :

library(knitr)
library(kableExtra)

data <- array(
 c(30, 10, 8, 12, 18, 7, 10, 25),
 dim = c(2, 2, 2),
 dimnames = list(
 "Begadang" = c("Ya", "Tidak"),
 "Pelupa" = c("Ya", "Tidak"),
 "Status" = c("Mahasiswa", "Pekerja")
 )
)
for (z in dimnames(data)$Status) {
 df <- as.data.frame.matrix(data[,,z])
 colnames(df) <- c("Ya", "Tidak")
 rownames(df) <- c("Ya", "Tidak")
 print(kable(df, format = "pandoc", booktabs = TRUE,
              caption = paste("Tabel untuk status =", z)) %>%
      kable_styling(full_width = FALSE, position = "center"))
}

## 
## 
## Table: Tabel untuk status = Mahasiswa
## 
##          Ya   Tidak
## ------  ---  ------
## Ya       30       8
## Tidak    10      12

## 
## 
## Table: Tabel untuk status = Pekerja
## 
##          Ya   Tidak
## ------  ---  ------
## Ya       18      10
## Tidak     7      25

# Uji Chi-Square untuk masing-masing strata usia
chisq_muda <- chisq.test(data[,,"Mahasiswa"])
chisq_tua <- chisq.test(data[,,"Pekerja"])
chisq_muda

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data[, , "Mahasiswa"]
## X-squared = 5.6071, df = 1, p-value = 0.01789

chisq_tua

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data[, , "Pekerja"]
## X-squared = 9.375, df = 1, p-value = 0.0022

Dengan ini dapat dilihat bahwa terdapat hubungan signifikan antara kebiasaan begadang dengan frekuensi sering lupa setelah dikontrol oleh variabel status.

3.4 Inferensi Tabel Kontingensi 3 Arah

3.4.1 Pengujian Statistik untuk melihat independensi Bersyarat

Indenpendensi bersyarat merupakan konsep penting dalam analisis tabel kontingensi 3 arah. Ini merujuk pada kondisi dimana dua variabel pada setiap level variabel ketiga. Secara matematis dapat ditulisankan sebagai :

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

Meskipun mereka X dan Y independent, tapi tidak menutup kemungkinan bahwa mereka tidak independet pada tabel marginalnya yang ini disebut sebagai simpson paradox.

Uji Mantel-Haenszel digunakan untuk menguji hubungan dua variabel sambil mengontrol variabel ketiga. Jadi disini dilihat apakah hubungan kebiasaan begadang dan sering lupa tetap signifikan setelah dikontrol berdasarkan status. Rumus dari CMH adalah sebagai berikut

\[ \chi^2_{\text{CMH}} = \frac{\left[ \sum_{k=1}^K \left( a_k - \frac{(a_k + b_k)(a_k + c_k)}{n_k} \right) \right]^2} {\sum_{k=1}^K \frac{(a_k + b_k)(c_k + d_k)(a_k + c_k)(b_k + d_k)}{n_k^2(n_k - 1)}} \]

Metode ini berguna dalam pengendalian faktor rancu ketika tabel kontingensi berlapis, mampu menguji hipotesis dua variabel akibat variabel ketiga serta bias dalam pembaurannya. Contoh perhitungannya dapat dilakukan dengan soal sebelumnya. Menggunakan fungsi mantelhaen.test pada R dan diperoleh hasil sebagai berikut

library(vcdExtra)
mantelhaen.test(data, correct = FALSE)

## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  data
## Mantel-Haenszel X-squared = 17.69, df = 1, p-value = 2.6e-05
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##   2.412071 12.089198
## sample estimates:
## common odds ratio 
##               5.4

Pada pengujian diatas, diperoleh beberapa kesimpulan :

1. p value = 2.6 * e-5, yang berarti ada bukti yang kuat bahwa kebiasaan begadang memliki hubungan yang signifikan dengan frekuensi lupa bahkan setelah dikontrol dengan berdasarkan status mahasiswa atau pekerja. Jadi dapat disimpulkan bahwa variabel X dan Y tetap signifikan berpengaruh
2. common odd ratios = 5,4 yang berarti individu yang memiliki kebiasaan begadang 5,4 kali lebih mungkin mengalami sering lupa daripada individu yang tidak punya kebiasaan begadang.

3.4.2 Odd Rasio Bersama

Odd rasio bersama ini diperoleh ketika terdapat tabel kontingensi dengan ukuran (2 x 2 x k) yang hal ini menyebabkan akan ada K odd rasio bersyarat. Jika odd rasio bersyarat ini relatif sama dan memiliki arah yang sama, maka bisa ditentukan nilai tunggal rasio yang disebut dengan odd rasio bersama. Formula odd rasio bersama ini adalah

penaksir

\[ \theta_{\text{MH}} = \frac{ \sum_{k=1}^{K} \frac{n_{11k} \cdot n_{22k}}{n_{\cdot \cdot k}} }{ \sum_{k=1}^{K} \frac{n_{12k} \cdot n_{21k}}{n_{\cdot \cdot k}} } \] Dengan:

• \(n_{ijk}\) = frekuensi sel baris \(i\), kolom \(j\), strata ke-\(k\)
• \(n_{\cdot \cdot k} = \sum_{i,j} n_{ijk}\) = total frekuensi pada strata ke-\(k\)

standar error

Standar error dari \(\log(\theta_{\text{MH}})\) didekati dengan:

\[ \sigma^2 \left[ \log(\theta_{\text{MH}}) \right] = \frac{ \sum\limits_{k=1}^{K} (n_{12k} + n_{21k}) \cdot \left( \frac{n_{12k} \cdot n_{21k}}{n_{\cdot \cdot k}^2} \right) } { 2 \left( \sum\limits_{k=1}^{K} \frac{n_{12k} \cdot n_{21k}}{n_{\cdot \cdot k}} \right)^2 } \]

Dengan:

• \(n_{ijk}\) = frekuensi sel baris \(i\), kolom \(j\), strata ke-\(k\)
• \(n_{\cdot \cdot k} = \sum_{i,j} n_{ijk}\)

Dengan \(w_k\) adalah bobot strata \(k\).

Interval kepercayaan \(100(1 - \alpha)\%\) untuk odds ratio Mantel-Haenszel diberikan oleh:

\[ \log(\theta_{\text{MH}}) \pm z_{\alpha/2} \cdot \sigma\left[ \log(\theta_{\text{MH}}) \right] \]

Sehingga, bentuk interval kepercayaan untuk \(\theta_{\text{MH}}\) adalah:

\[ \left( \exp\left[ \log(\theta_{\text{MH}}) - z_{\alpha/2} \cdot \sigma \right], \quad \exp\left[ \log(\theta_{\text{MH}}) + z_{\alpha/2} \cdot \sigma \right] \right) \]

Sebagai contoh, tersedia data berikut

City	Smooking	Kanker (YA)	Kanker (NO)
Beijing	Smoker	126	100
Beijing	Non Smoker	35	61
Shanghai	Smoker	908	688
Shanghai	Non Smoker	497	807
Shenyang	Smoker	913	747
Shenyang	Non Smoker	336	598
Nanjing	Smoker	235	172
Nanjing	Non Smoker	58	121
Harbin	Smoker	402	308
Harbin	Non Smoker	121	215
Zhengzou	Smoker	182	156
Zhengzou	Non Smoker	72	98
Taiyuun	Smoker	60	99
Taiyuun	Non Smoker	11	43
Nanchang	Smoker	104	89
Nanchang	Non Smoker	21	36

Diketahui data hasil peneltiain 8 provinsi di China mengenai kebiasaan merokok (X), kanker paru (Y) dan provinsi (Z). Hitung odd rasio bersamanya.

library(epitools) # Paket untuk menghitung odds ratio
# Data dari 8 kota dalam bentuk array
lung_cancer_data <- array(
c(126, 35, 61, 100,
 908, 497, 688, 807,
 913, 336, 747, 598,
 235, 58, 172, 121,
 402, 121, 308, 215,
 182, 72, 156, 98,
 60, 11, 99, 43,
 104, 21, 89, 36),
 dim = c(2, 2, 8),
 dimnames = list(
 Smoking = c("Smokers", "Nonsmokers"),
 Lung_Cancer = c("Yes", "No"),
 City = c("Beijing", "Shanghai", "Shenyang", "Nanjing", "Harbin",
 "Zhengzhou", "Taiyunn", "Nanchang")
 )
)
# Menghitung odds ratio bersama
mh_test <- mantelhaen.test(lung_cancer_data, correct = FALSE)
print(mh_test)

## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  lung_cancer_data
## Mantel-Haenszel X-squared = 308.89, df = 1, p-value < 2.2e-16
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  2.056753 2.469706
## sample estimates:
## common odds ratio 
##          2.253791

Dari output tersebut dapat dilihat bahwa odd rasio bersamanya adalah 2,25 yang mengartikan bahwa orang yang merokok memiliki peluang 2,25x lebih besar untuk terkena kanker daripada orang yang tidak merokok.

3.4.3 Uji Homogenitas Odd Rasio

Homogenitas ini terjadi jika odd rasio pada setiao tabel parsial bernilai sama

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

Jika odd rasio berbeda di setiap kontrol (Z) maka terdapat interaksi antara X dan Y dalam hubungannya terhadap Z dan berlaku sebaliknya. Hipotesis dari pengujian ini adalah

H0 : odd rasio sama pada seluruh strata

H1 : Setidaknya ada satu odd rasio yang berbeda

Untuk menguji homogenitas, digunakan uji Breslow-Day. Ini ditujukan untuk melihat apakah kebiasaan begadang dan frekuensi lupa berbeda secara signifikan antara dua kelompok mahasiswa dan bekerja. Rumus dari Breslow Day adalah sebagai berikut

\[ Q_{\text{BD}} = \sum_{k=1}^K \frac{\left( \hat{\text{OR}}_k - \hat{\text{OR}}_{\text{MH}} \right)^2}{\text{Var}(\hat{\text{OR}}_k)} \]

Stastitik ini mengikuti distribusi chi square dengan derajat bebas = k - 1 sehingga pengambilan keputusannya dalapt dilakukan dengan perbandingan terhadap nilai chi square tabel.

library(DescTools)
library(vcdExtra)
breslow_test <- BreslowDayTest(data)
print(breslow_test)

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

Pada hasil di atas diperoleh p-value bernilai 0,6652 yang mana keputusannya adalah menerima H0. Jadi antara kebiasaan begadang dan frekuensi lupa pada kelompok mahasiswa dan pekerja relatif sama alias homogen.

Dari beberapa contoh soal mengenai kasusu kebiasaan begadang terhadap frekuensi lupa dengan kontrol status sebagai berikut

1. Individu yang memiliki kebiasaan begadang lebih sering mengalami lupa daripada individu yang tidak biasa begadang, baik itu mahasiswa maupun pekerja.
2. kebiasaan begadang memliki hubungan yang signifikan dengan frekuensi lupa bahkan setelah dikontrol dengan berdasarkan status mahasiswa atau pekerja.
3. antara kebiasaan begadang dan frekuensi lupa pada kelompok mahasiswa dan pekerja relatif sama alias homogen.

3.5 Review

Pada Bab ini dicobakan beberapa pengerjaan dengan penjelasan sebagai berikut

1. Tabel Kontingensi 3 Arah memberikan 3 variabel yakni variabel prediktor (X), variabel Respon (Y) dan Variabel yang dikontrol (Z). Analisis ini dikontrol oleh Z untuk melihat apakah ada pengaruh dari variabel lain selain variabel X dan Y.
2. Tabel kontingensi 3 arah juga memiliki yang namanya peluang bersama, peluang marginal dan peluang bersyarat serta juga memilki ukuran asosiasi yakni beda peluang, resiko relatif dan odd rasio.
3. Dalam Tabel Kontingensi 3 Arah terdapat sebuah ukuran yang disebut dengan conditional independent. Ukuran ini bertujuan untuk melihat bagaimana hubungan antara X dan Y dikontrol dengan variabel Z.
4. Pengujian Independensi Bersyarat adalah lanjutan dari bahasan sebelumnya yakni menguji apakah hubungan X dan Ysignifikan dengan kontrol dari variabel Z,. Pengujian ini menggunakan statistik uji Cochran Mantel Hanzel (CMH)
5. Pengujian odd rasio bersama bertujuna untuk melihat apakah ada asosiasi antara variabel X dan Y bersyarat variabel Z
6. Uji Homogenitas odd rasio dilakukan untuk menguji apakah odd rasio berbeda signifikan dari kontrol yang diberikan variabel Z. Pengujian ini dilakukan dengan statistik uji Breslow-Day.

4 Generalized Linear Model

4.1 Definisi Umum

Generalized Linear Model (GLM) adalah suatu kerangka statistik yang memperluas model regresi linear dengan memungkinkan variabel respon mengikuti distribusi yang tidak harus normal, seperti binomial, Poisson, atau gamma. GLM menghubungkan nilai harapan dari variabel respon dengan kombinasi linear dari variabel prediktor melalui suatu fungsi link. Model ini digunakan untuk menganalisis berbagai jenis data, terutama ketika asumsi normalitas tidak terpenuhi.

4.2 Komponen GLM

GLM memiliki 3 komponen utama yakni

1. Komponen Acak

   Nilai pengamatan variabel Y yang saling bebas dan berdistribusi keluarga eksponensial. Jadi variabel responnya mengikuti distribusi-distribusi yang termasuk ke dalam keluarga eksponen, contohnya distribusi normal dan poisson mengikuti kaidah exponential family. Untuk menentukan apakah sebuah distribusi mengikuti keluarga eksponensial, dapat dilakukan dengan pengujian atau pembuktian rumus berikut

   \[ f(y; \theta) = \exp\left(a(y) \cdot b(\theta) + c(\theta) + d(y)\right) \]

   Jika variabel acak dengan distribusi tertentu dapat dinyatakan dalam bentuk di atas, maka dinyatakan bahwa distribusi tersebut merupakan bagian dari exponential family. Contoh \[ P(Y = y) = p^y (1 - p)^{1 - y} \]

   Dengan substitusi \(p = \frac{e^\theta}{1 + e^\theta}\), maka bentuk exponential family-nya menjadi: \[ P(Y = y) = \exp \left( y \cdot \log\left( \frac{\pi}{1 - \pi} \right) + n \cdot \log(1 - \pi) + \log \binom{n}{y} \right) \]

   Dengan ini dapat dilihat bahwa distribusi binomial mengikuti kaidah keluarga eksponen

2. Komponen sistematik

   Satu set parameter B dan variabel explanatory X yang membentuk kombinasi linear : \[ \eta = \alpha + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p \]

   Bisa dikatakan bahwa ini adalah bentuk dari model linearnya.

3. Fungsi Penghubung (Link Function)

   Merupakan fungsi penghubung antara komponen sistematik dengan nilai ekspektasi dari komponen acak.

4.3 GLM untuk Respon Biner

Untuk nilai respon biner dapat dianalisis melalui model regresi logistik. Regresi Logistik hampir sama dengan model linear yang mana menginput kombinasi koefisien untuk menghasilkan prediksi, namun bedanya model ini hanya untuk data biner yang nanti prediksinya akan mengikuti kurva S (S-Shape). Sebagai contoh analisis terhadap kemungkinana seorang siswa diterima di perguruan tinggi tertentu berdasarkan nilai A, B dan C. Outputnya adalah perkiraan probabilitas seorang siswa diterima dengan menganalisis variabel variabel penentunya.

Secara teknis, Regresi Logistik lebih mudah diterapkan dalam metode pembelajaran. Regresi ini menggunakan fungsi sigmoid dalam memetakan prediksi dan probabilitasnya. Fungsi Sigmoid sendiri adalah kurva berbentuk S yang mengubah setiap nilai riil menjadi rentang antara 0 dan 1. Formula fungsi sigmoid sebagai berikut

\[ \mu = \frac{1}{1 + e^{-\eta}} \]

Jika hasil fungsi lebih dari ambang batas, maka model akan diperkirakan masuk ke kelompok tertentu. Misal ambang batas yang ditetapkan adalah 0.5. Ketika hasil nya kurang dari 0.5, maka dianggap masuk kelas bawah (biner = 0) dan sebaliknya.

Regresi logistik memiliki sejumlah keunggulan yang membuatnya populer dalam penerapan pembelajaran mesin. Salah satu alasan utamanya adalah kemudahan dalam implementasi dan efisiensi komputasi. Proses pembelajaran mesin umumnya melibatkan tahap pelatihan dan pengujian, di mana model dilatih untuk mengenali pola dari data masukan dan mengaitkannya dengan keluaran yang sesuai. Regresi logistik dapat dilatih dengan cepat dan tidak membutuhkan sumber daya komputasi yang besar, sehingga sangat cocok untuk digunakan terutama pada tahap awal pengembangan model. Selain itu, regresi logistik juga mudah untuk ditafsirkan, menjadikannya pilihan yang tepat bagi peneliti atau praktisi yang menginginkan hasil analisis yang transparan dan dapat dijelaskan secara logis.

Keunggulan lainnya terletak pada kemampuannya untuk bekerja secara efektif pada data yang dapat dipisahkan secara linear, yaitu ketika dua kelas dalam data dapat dipisahkan dengan garis lurus. Dalam konteks klasifikasi biner, regresi logistik mampu memetakan probabilitas suatu kejadian menggunakan fungsi logistik atau sigmoid, yang menghasilkan output dalam rentang 0 hingga 1. Fungsi ini memungkinkan prediksi probabilistik yang berguna dalam banyak aplikasi nyata. Selain itu, model ini juga memberikan informasi yang bernilai tentang pengaruh masing-masing variabel prediktor terhadap hasil, baik dari segi besar pengaruh maupun arah hubungannya, apakah positif atau negatif. Hal ini membuat regresi logistik tidak hanya berguna untuk prediksi, tetapi juga untuk eksplorasi dan pemahaman struktur data.

Simulasi Struktur Data dan Visualisasi

Sebuah perusahaan sedang meneliti hubungan antara skor psikotes pelamar dengan kemungkinan diterima bekerja. Peneliti mensimulasikan data untuk 50 pelamar, di mana skor psikotes (variabel x) diasumsikan berkisar dari sangat rendah (-7) hingga sangat tinggi (+7).

library(ggplot2)
#Simulasi Data Logistik
set.seed(28)
n <- 50
x <- seq(-7,7,length.out = n)
log.odds <- (-0.5) + 0.5 * x
prob <- 1 / (1 + exp(-log.odds))
y <- rbinom(n,1,prob)

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

##           x y       prob
## 1 -7.000000 0 0.01798621
## 2 -6.714286 0 0.02069111
## 3 -6.428571 0 0.02379294
## 4 -6.142857 0 0.02734679
## 5 -5.857143 0 0.03141437
## 6 -5.571429 0 0.03606454

#Visualisasi 
plot (x,y,pch = 16, col ="gray60",
      xlab = "X", ylab = "Probabilities",
      main = "Simulasi Regresi Logistik dengan kurva sigmoid")
lines(x,prob, col = "blue", lwd = 2)
abline( h=0.5 , col = "red", lty = 2)

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

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

Model ini sangat efektif digunakan pada data yang dapat dipisahkan secara linear, yaitu ketika dua kelas dapat dipisahkan oleh sebuah garis lurus. Regresi logistik menggunakan fungsi logistik atau sigmoid untuk memetakan nilai input menjadi probabilitas antara 0 dan 1. Secara matematis, fungsi link logit dituliskan sebagai:
\[ g(\mu) = \log\left(\frac{\mu}{1 - \mu}\right) \]
yang kemudian dihubungkan dengan prediktor linear sebagai berikut:
\[ \log\left(\frac{\mu}{1 - \mu}\right) = X \cdot \beta \]
Persamaan ini memungkinkan model untuk mengukur pengaruh variabel prediktor terhadap probabilitas keluaran serta memberikan informasi arah dan kekuatan hubungannya.

Dengan ini diperoleh fungsi inverse linknya sebagai berikut

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

4.3.1 Estimasi Parameter

Metode estimasi pada GLM umumnya menggunakan Maximum Likelihood Estimation dengan fungsi loglikelihood fungsi untuk regresi logistik :

\[ \mathcal{L}(\beta) = \sum_{i=1}^{n} \left[ y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \right] \] dengan : \[ \pi_i = \frac{\exp(\mathbf{x}_i^T \beta)}{1 + \exp(\mathbf{x}_i^T \beta)} \] Estimasi ini dilakukan dengan iterasi Newton-Raphson.

Contoh perhitungan dalam R

Pada contoh ini, kami akan memanfaatkan data bawaan mtcars dari paket dasar R, dengan terlebih dahulu mengubah variabel target menjadi variabel biner. Tujuan pemodelan adalah untuk memprediksi jenis transmisi mobil (apakah otomatis (0) atau manual (1)) dengan menggunakan variabel prediktor seperti konsumsi bahan bakar (mpg) dan tenaga mesin (hp).

data(mtcars)
df <- mtcars
df$am <- factor(df$am, labels = c("Automatic", "Manual"))
summary(df[, c("am", "mpg", "hp")])

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

#Estimasi Model Regresi Logistik
model_real <- glm(am ~ mpg + hp, data = df, family = binomial)
summary(model_real)

## 
## Call:
## glm(formula = am ~ mpg + hp, family = binomial, data = df)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -33.60517   15.07672  -2.229   0.0258 *
## mpg           1.25961    0.56747   2.220   0.0264 *
## hp            0.05504    0.02692   2.045   0.0409 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 43.230  on 31  degrees of freedom
## Residual deviance: 19.233  on 29  degrees of freedom
## AIC: 25.233
## 
## Number of Fisher Scoring iterations: 7

exp(coef(model_real))

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

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

library(ggplot2)
ggplot(df, aes(x =mpg, y = hp, color = pred_class, shape = am)) +
  geom_point(size = 3) +
  labs(title = "Prediksi Regresi Logistik: Transmisi Mobil",
       x="Mil per Galon (mpg)",
       y="Tenaga Kuda (hp)",
       color = "Prediksi", shape = "Fakta") +
  theme_minimal()

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

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

Confusion matrix tersebut menunjukkan performa model dalam memprediksi jenis transmisi mobil berdasarkan variabel mpg dan hp. Dari total observasi, model berhasil mengklasifikasikan 16 mobil bertransmisi otomatis dan 10 mobil bertransmisi manual dengan benar. Terdapat kesalahan klasifikasi sebanyak 3 unit untuk masing-masing kategori, yakni 3 mobil manual diprediksi sebagai otomatis, dan sebaliknya. Secara umum, ini menunjukkan bahwa model cukup andal dalam membedakan dua jenis transmisi berdasarkan fitur yang digunakan.

4.4 GLM untuk respon cacah (counting)

Regresi Poisson adalah metode statistik yang umum digunakan untuk menganalisis data cacah atau data hitung, yaitu data yang menyatakan jumlah kejadian dalam suatu periode waktu atau wilayah tertentu. Model ini sangat cocok untuk variabel respons yang berupa bilangan bulat non-negatif seperti jumlah kecelakaan lalu lintas, jumlah panggilan telepon, atau jumlah penyakit pada populasi. Dalam regresi Poisson, diasumsikan bahwa variabel respons mengikuti distribusi Poisson, di mana rata-rata kejadian (parameter lambda) bergantung pada variabel prediktor melalui fungsi link logaritmik. Dengan kata lain, log dari nilai harapan jumlah kejadian dimodelkan sebagai kombinasi linier dari variabel-variabel independen, sehingga hubungan antara prediktor dan respons menjadi eksponensial.

Salah satu keunggulan utama regresi Poisson adalah kemampuannya untuk menangani data dengan karakteristik variansi yang sama dengan rata-rata, yang merupakan ciri khas distribusi Poisson. Model ini memungkinkan peneliti untuk memahami dan memprediksi faktor-faktor yang memengaruhi frekuensi kejadian, serta menilai signifikansi pengaruh tiap variabel prediktor secara simultan. Namun, dalam praktiknya, sering ditemukan kasus overdispersi, yaitu variansi data lebih besar dari rata-rata, yang dapat menyebabkan model Poisson standar kurang tepat. Dalam situasi seperti ini, metode alternatif seperti regresi negatif binomial bisa digunakan sebagai solusi untuk mengatasi masalah overdispersi tersebut. Regresi Poisson tetap menjadi alat analisis yang penting dalam ilmu sosial, kesehatan masyarakat, dan bidang lain yang memerlukan pemodelan data hitung.

Regresi ini cocok digunakan pada data respon cacah. Selain itu regresi poisson mengikuti keluarga eksponensial dengan

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

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

Fungsi Link kanoniknya adalah

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

Hubungan antara prediktor dan respons dalam bentuk linear kemudian dinyatakan sebagai:

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

Sehingga bentuk eksponensial dari model regresi Poisson menjadi:

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

Model ini memungkinkan kita untuk menafsirkan efek dari setiap variabel prediktor terhadap jumlah kejadian dalam skala logaritmik, serta memperkirakan nilai ekspektasi dari variabel cacah berdasarkan kombinasi linier prediktor.

Dalam Estimasi parameternya, dilakukan juga dengan Maximum Likelihood Estimation (MLE) dengan fungsinya sebagai berikut

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

Nilai 𝛽 dapat diperoleh melalui metode numerik seperti iterasi Newton-Raphson

Contoh perhitungannya

Sebuah perusahaan ritel ingin memahami bagaimana efektivitas iklan online memengaruhi jumlah pelanggan yang datang ke toko fisik mereka. Variabel x merepresentasikan skor efektivitas iklan online, yang diperoleh melalui analisis internal terhadap berbagai indikator kampanye pemasaran digital.

Untuk mensimulasikan hubungan ini, peneliti menghasilkan 200 data acak skor iklan (x) dari distribusi normal. Lalu, mereka mengasumsikan bahwa jumlah pelanggan yang datang ke toko dalam sehari (y) mengikuti distribusi Poisson

set.seed(28)
n <- 200
x <- rnorm(n)
lambda <- exp(0.3 + 0.6 * x)
y <- rpois(n, lambda)
data <- data.frame(y, x)
head(data)

##   y           x
## 1 0 -1.90215722
## 2 1 -0.06429479
## 3 0 -1.33116707
## 4 1 -1.81999167
## 5 4  0.16266969
## 6 2  0.53139634

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

## 
## Call:
## glm(formula = y ~ x, family = poisson, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.20773    0.06738   3.083  0.00205 ** 
## x            0.61181    0.05946  10.290  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 315.07  on 199  degrees of freedom
## Residual deviance: 206.92  on 198  degrees of freedom
## AIC: 552.44
## 
## Number of Fisher Scoring iterations: 5

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

Dari data simulasi tersebut, data dilihat bahwa intercept dan variabel prediktor(efektivitas iklan) memberikan pengaruh perubahan yang signifikan terhadap variabel respon(jumlah konsumen).

4.4.1 Diagnostik dan Overdispersion

Overdispersion terjadi ketika varians lebih besar daripada rata-ratanya karena pada distribusi poisson, nilai rata rata = varians. Jika nilanya lebih dari satu (>1) maka terjadi overdispersion dan negative binomial regression adalah alternatif solusi modelnya. Untuk mendeteksinya digunakan rumusan berikut

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

## [1] 0.916608

Pada pengujian tersebut, diperoleh nilai yang kurang dari 1, maka tidak terjadi overdispersion dan data simulasi ini cocok menggunakan regresi poisson.

Contoh Perhitungan 2

Regresi Poisson digunakan untuk memodelkan data cacah (count data), yaitu data respons berupa bilangan bulat non-negatif. Dalam contoh ini, kita akan menggunakan dataset warpbreaks dari R base package yang berisi jumlah patahan benang pada mesin tenun.

data("warpbreaks")
head(warpbreaks)

##   breaks wool tension
## 1     26    A       L
## 2     30    A       L
## 3     54    A       L
## 4     25    A       L
## 5     70    A       L
## 6     52    A       L

summary(warpbreaks)

##      breaks      wool   tension
##  Min.   :10.00   A:27   L:18   
##  1st Qu.:18.25   B:27   M:18   
##  Median :26.00          H:18   
##  Mean   :28.15                 
##  3rd Qu.:34.00                 
##  Max.   :70.00

poisson_model <- glm(breaks ~ wool + tension, data = warpbreaks, family = poisson)
summary(poisson_model)

## 
## Call:
## glm(formula = breaks ~ wool + tension, family = poisson, data = warpbreaks)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.69196    0.04541  81.302  < 2e-16 ***
## woolB       -0.20599    0.05157  -3.994 6.49e-05 ***
## tensionM    -0.32132    0.06027  -5.332 9.73e-08 ***
## tensionH    -0.51849    0.06396  -8.107 5.21e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 297.37  on 53  degrees of freedom
## Residual deviance: 210.39  on 50  degrees of freedom
## AIC: 493.06
## 
## Number of Fisher Scoring iterations: 4

exp(coef(poisson_model))

## (Intercept)       woolB    tensionM    tensionH 
##  40.1235380   0.8138425   0.7251908   0.5954198

warpbreaks$predicted <- predict(poisson_model, type = "response")
library(ggplot2)
ggplot(warpbreaks, aes(x =tension, y=breaks, color = wool)) +
  geom_jitter(width = 0.2, alpha = 0.6) +
  geom_point(aes(y =predicted), shape = 18, size = 3, color = "black") +
  facet_wrap(~wool) +
  labs(title = "Prediksi Jumlah Patahan Benang berdasarkan Wol dan Ketegangan",
       x="Tingkat Ketegangan",
       y="Jumlah Patahan Benang",
       color = "Jenis Wol") +
  theme_minimal()

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

Output regresi Poisson di atas menunjukkan bahwa variabel kategorikal wool dan tension berpengaruh signifikan terhadap jumlah patahan benang (breaks). Dari koefisien nya, terlihat bahwa kedua variabel ini memberikan koefisien negatif yang berarti semakin banyak penggunaan kedua variabel ini, maka akan menurunkan rata-rata jumlah patahan secara signifikan.

4.5 Inferensi GLM

Inferensi statistik dalam GLM membutuhkan pemahaman terhadap ekspektasi dan varians dari estimator model, terutama untuk mengembangkan uji uji seperti likelihood dan interval kepercayaan.

1. Ekspektasi Estimator, menunjukkan bahwa apakah suatu estimator adalah tak bias.
2. Varians Estimator, menunjukkan presisi. Dalam GLM varians tidak konstan sehingga diperlukan adanya sebauh variabel yang disebut dengan parameter dispersi.

4.5.1 Ekspektasi dan Varians dalam GLM

Distribusi dalam GLM biasanya mengikuti keluarga eksponensial dengan bentuk

\[ f(y;\theta) = \exp\left(a(y)\cdot b(\theta) + c(\theta) + d(y)\right) \]

Untuk menentukan ekspektasi dan varians dalam GLM, harus memenuhi bentuk berikut

1. Ekspektasi

\[ \mu = b'(\theta) = \mu = \frac{ -c'(\theta) }{ b'(\theta) } \]

2. Varians

\[ \text{Var}(Y) = b''(\theta) = \phi V(\mu) = \frac{b''(\theta) \cdot c'(\theta) - c''(\theta) \cdot b'(\theta)}{[b'(\theta)]^3} \]

4.5.2 Metode penaksiran parameter

Metode penaksiran dapat dilakukan dengan beberapa metode.

1. Maximum Likelihood Estimation

   Memaksimumkan fungsi likelihood. namun karena GLM tidak eksplisit maka digunakan metode numeric.

2. Newton-Raphson yang menggunakan turunan pertama fungsi likelihood yakni score vektor (gradien) dan turunan kedua, yakni Hessian matrix.

   \[ \hat{\boldsymbol{\beta}} = \boldsymbol{\beta}^{(t)} - \left[ \mathbf{H}(\boldsymbol{\beta}^{(t)}) \right]^{-1} \mathbf{U}(\boldsymbol{\beta}^{(t)}) \]

3. Fisher Scoring, mengganti Hessian matrix pada Newton-Raphson dengan matriks informasi Fisher. Metode ini dapat digunakan untuk menentukan titik maksimum dari suatu fungsi sejalan dengan permasalahan dalam menentukan taksiran parameter MLE.

4. IRLS, memodifikasi Fisher sehingga estimasi mirip dengan Least Square.

4.5.3 Diagnostik Model GLM

Diagnostik ini dilakukan untuk mengevaluasi apakah model sudah tepat. Ini dapat dilakukan dengan beberapa analisis seperti uji formal atau dengan grafik nilai prediksi vs nilai aktual.

Pada uji formal, dapat dilakukan dengan beberapa metode seperti

1. Statistik Devians

   Statistik devians ini mengukur apakah ada model yang lebih baik. Jika nilai devians besar, maka model keungkinan besar tidak cocok. Pada pengujiannya, model dibandingkan dengan saturated model dengan statistik berikut :

   \[ D = 2 \sum_{i=1}^{n} \left[ y_i \ln\left(\frac{y_i}{\hat{\mu}_i}\right) + (1 - y_i)\ln\left(\frac{1 - y_i}{1 - \hat{\mu}_i}\right) \right] \]

   Keterangan: - \(y_i\) = data aktual pada observasi ke-\(i\) - \(\hat{\mu}_i\) = nilai prediksi dari model pada observasi ke-\(i\) - \(n\) = jumlah observasi - Jika \(y_i = 0\) atau \(y_i = 1\), maka bagian log tersebut disesuaikan agar tidak terjadi log(0)

2. Uji Chi-Square Pearson.

   Menguji apakah model lebih baik daripada tidak ada model sama sekali. Jika nilai signifikan maka model lebih baik daripada tanpa model.

   \[ X^2 = \sum_{i=1}^{k} \frac{(O_i - E_i)^2}{E_i} \]

   Keterangan: - \(O_i\) = frekuensi yang diamati pada kategori ke-\(i\) - \(E_i\) = frekuensi yang diharapkan pada kategori ke-\(i\) - \(k\) = jumlah kategori - Nilai \(X^2\) dibandingkan dengan nilai kritis pada distribusi chi-square dengan \(k - 1\) derajat kebebasan

4.5.4 Detail Metode Estimasi dan Inferensi Regresi Logistik

Seperti di bahasain sebelumnya bahwa regresi logistik berfokus pada pembuatan model estimasi untuk data respon biner. Estimasi menggunakan MLE karena model tidak linear dalam parameternya.

Untuk estimasi menggunakan Metode Newton-Raphson yang mana metode ini mencari nilai koefisien untuk memaksimalkan fungsi log likelohood pada model regresi. Metode Log-Likelihood untuk n observasi :

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

Model regresi logistik untuk probabilitas :

\[ p_i = \frac{1}{1 + e^{-x_i^\top \beta}} \]

Dalam proses pengerjaannya, ada proses iterasi yang menggunakan formula berikut.

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

\[ W^{(t)} = \text{diag}(\pi_i(1 - \pi_i)) \]

Estimasi parameter pada model regresi logistik dilakukan dengan MLE. Metode Newton-Raphson adalah metode untuk memaksimalkan log-likelihood yang langkahlangkahnya merupakan penurunan. Penurunan pertama merupakan scoring dan penurunan kedia berupa matrix Hessian.

Sebagai contoh berikut adalah simulasi Estimasi MLE dengan Newton-Raphson

# Data simulasi
set.seed(28)
x <- rnorm(100)
beta_true <- c(-1, 2)
X <- cbind(1, x)
eta <- X %*% beta_true
p <- 1 / (1 + exp(-eta))
y <- rbinom(100, 1, p)

Newton-Raphson Iterasi manual

beta <- c(0, 0)
tol <- 1e-6
max_iter <- 100
for (i in 1:max_iter) {
  eta <- X %*% beta
  pi_hat <- 1 / (1 + exp(-eta))
  W <- diag(as.numeric(pi_hat * (1- pi_hat)))
  z <- eta + solve(W) %*% (y- pi_hat)
  beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
  if (sum(abs(beta_new- beta)) < tol) break
  beta <- beta_new
}
beta

##         [,1]
##   -0.8943083
## x  2.1833646

Dari syntax ini diperoleh nilai beta untuk intercept dan beta untuk variabel prediktor. Selanjutnya akan diperiksa apakah estimasi nilai beta ini memberikan pengaruh yang signifikan dalam analisis

4.5.4.1 Inferensi Parameter

Uji Inferensi ini berguna untuk mengevaluasi apakah taksirannya signifikan yang menandakan apakah variabel predikstor tersebut memang berpengaruh terhadap prediksi.

1. Uji Wald

   Uji Wald adalah metode statistik yang digunakan untuk menguji signifikansi parameter dalam model regresi, baik itu regresi linier maupun model regresi lainnya, seperti regresi logistik atau regresi Poisson. Uji ini menguji apakah suatu parameter dalam model (misalnya koefisien regresi) berbeda secara signifikan dari suatu nilai tertentu (biasanya nol). Statistik uji Wald adalah

   \[ W = \left( \frac{\hat{\beta}}{\text{SE}(\hat{\beta})} \right)^2 \]

   Keterangan: - \(\hat{\beta}\) = estimasi koefisien regresi - \(\text{SE}(\hat{\beta})\) = standar error dari \(\hat{\beta}\) - Nilai \(W\) mengikuti distribusi chi-square dengan derajat kebebasan 1 (untuk satu koefisien).

   Nilai \(W\) yang besar menunjukkan bahwa koefisien \(\hat{\beta}\) berbeda signifikan dari nol. Bandingkan nilai \(W\) dengan nilai kritis dari distribusi chi-square pada taraf signifikansi tertentu (misal: 0.05). Jika \(W\) lebih besar dari nilai kritis, maka koefisien tersebut signifikan secara statistik.

   Contoh Perhitungannya

   Sebuah studi dilakukan untuk meneliti apakah tingkat stres (diukur dalam bentuk skor x, dari distribusi normal) berpengaruh terhadap kemungkinan seseorang mengalami insomnia (y), yaitu 1 jika insomnia terjadi dan 0 jika tidak.

   Simulasi dilakukan dengan 100 responden

   set.seed(28)
   n <- 100
   x <- rnorm(n)
   log_odds <--0.5 + 1.2 * x
   p <- 1 / (1 + exp(-log_odds))
   y <- rbinom(n, 1, p)
   data <- data.frame(x, y)
   model <- glm(y ~ x, data = data, family = binomial)
   summary(model)

   ## 
   ## Call:
   ## glm(formula = y ~ x, family = binomial, data = data)
   ## 
   ## Coefficients:
   ##             Estimate Std. Error z value Pr(>|z|)    
   ## (Intercept)  -0.4280     0.2327  -1.839   0.0659 .  
   ## x             1.2162     0.2943   4.133 3.58e-05 ***
   ## ---
   ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   ## 
   ## (Dispersion parameter for binomial family taken to be 1)
   ## 
   ##     Null deviance: 135.37  on 99  degrees of freedom
   ## Residual deviance: 110.72  on 98  degrees of freedom
   ## AIC: 114.72
   ## 
   ## Number of Fisher Scoring iterations: 4

   beta_hat <- coef(model)["x"]
   se_beta <- summary(model)$coefficients["x", "Std. Error"]
   
   Z <- beta_hat / se_beta
   Z

   ##        x 
   ## 4.132766

   Wald_stat <- Z^2
   Wald_stat

   ##        x 
   ## 17.07975

   p_value <- 1- pchisq(Wald_stat, df = 1)
   p_value

   ##            x 
   ## 3.584238e-05

   Kesimpulannya Dengan uji Wald diperoleh nilai p-value < alpha. Disimpulkan bahwa keofisiennya signifikan yang berarti prediktor(stress) memberikan pengaruh(insomnia).

2. Uji likelihood Ratio

   Uji Likelihood Ratio (LR) adalah metode statistik yang digunakan untuk membandingkan dua model yang berbeda, yaitu model yang lebih kompleks (model alternatif) dan model yang lebih sederhana (model nol atau model pembanding). Uji ini menguji sejauh mana model yang lebih kompleks dapat menjelaskan data dengan lebih baik dibandingkan model yang lebih sederhana.

   Contoh

   # Model null
   model_null <- glm(y ~ 1, data = data, family = binomial)
   # Likelihood ratio test
   anova(model_null, model, test = "Chisq")

   ## Analysis of Deviance Table
   ## 
   ## Model 1: y ~ 1
   ## Model 2: y ~ x
   ##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
   ## 1        99     135.37                          
   ## 2        98     110.72  1    24.65 6.873e-07 ***
   ## ---
   ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

   Dengan hasil ini dinyatakan pula bahwa model null juga memberikan pengaruh yang signifikan. Ini biasa juga digunakan untuk memlih model berdasarkan nilai deviance. Jika nilai deviancenya besar dan positif, maka model kompleks memang lebih baik, namun kalau nilainya kecil itu menandakan baik model kompleks ataupun sederhana tetap memberikan pengaruh yang tidak jauh. Oleh karena itu model sederhana juga baik.

   Evaluasi model dengan

   AIC(model)

   ## [1] 114.7213

   BIC(model)

   ## [1] 119.9317

   Kedua nilai ini digunakan untuk dibandingkan dengan nilai AIC dan BIC model lain. Semakin kecil nilainya maka semakin baik modelnya.

4.5.4.2 Contoh regresi logistik multiple

Misal kita melakukan pemodeltan terhadap data mtcars yang ada di R studio. Kita ingin memodelkan kemungkinan mobil memiliki efisiensi bahan bakar berdasarkan beberapa prediktor.

Pada tahap awal, dataset mtcars diambil dari pustaka bawaan R dan dilakukan transformasi terhadap variabel mpg (mile per gallon). Tujuannya adalah mengubah mpg menjadi variabel kategorik biner bernama mpg_high, yaitu 1 jika nilai mpg lebih tinggi dari median dan 0 jika tidak. Ini memungkinkan kita memodelkan kemungkinan mobil termasuk kategori efisien atau tidak. Kemudian, hanya beberapa variabel prediktor yang dipilih, yaitu wt (berat mobil), hp (horsepower/tenaga), dan am (tipe transmisi), untuk digunakan dalam analisis regresi logistik.

library(dplyr)

data(mtcars)
mtcars <- mtcars %>%
 mutate(mpg_high = ifelse(mpg > median(mpg), 1, 0)) %>%
 select(mpg_high, wt, hp, am) # pilih beberapa prediktor
head(mtcars)

##                   mpg_high    wt  hp am
## Mazda RX4                1 2.620 110  1
## Mazda RX4 Wag            1 2.875 110  1
## Datsun 710               1 2.320  93  1
## Hornet 4 Drive           1 3.215 110  0
## Hornet Sportabout        0 3.440 175  0
## Valiant                  0 3.460 105  0

Model regresi logistik dibentuk untuk memprediksi probabilitas bahwa mpg_high bernilai 1 (efisien) berdasarkan wt, hp, dan am. Fungsi glm() digunakan dengan keluarga distribusi binomial, karena targetnya bersifat biner. Hasil summary(model_logit) memberikan informasi tentang nilai koefisien regresi, standar error, nilai z, dan p-value untuk setiap prediktor. Ini membantu kita memahami pengaruh masing-masing variabel terhadap kemungkinan efisiensi bahan bakar.

model_logit <- glm(mpg_high ~ wt + hp + am, data = mtcars, family = binomial)
summary(model_logit)

## 
## Call:
## glm(formula = mpg_high ~ wt + hp + am, family = binomial, data = mtcars)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  6.244e+02  4.856e+05   0.001    0.999
## wt          -1.847e+02  1.875e+05  -0.001    0.999
## hp          -8.852e-02  2.045e+03   0.000    1.000
## am          -6.038e+01  1.044e+05  -0.001    1.000
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4.4236e+01  on 31  degrees of freedom
## Residual deviance: 3.4224e-09  on 28  degrees of freedom
## AIC: 8
## 
## Number of Fisher Scoring iterations: 25

Untuk menguji apakah model dengan prediktor (model_logit) secara signifikan lebih baik daripada model tanpa prediktor (model_null), digunakan uji Likelihood Ratio. Fungsi anova() dengan opsi test = "Chisq" membandingkan kedua model. Jika hasil uji menunjukkan nilai p yang signifikan, maka dapat disimpulkan bahwa menambahkan prediktor secara keseluruhan memang memberikan informasi tambahan yang signifikan dalam menjelaskan efisiensi bahan bakar.

model_null <- glm(mpg_high ~ 1, data = mtcars, family = binomial)
anova(model_null, model_logit, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: mpg_high ~ 1
## Model 2: mpg_high ~ wt + hp + am
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        31     44.236                          
## 2        28      0.000  3   44.236 1.344e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Uji Wald digunakan untuk menguji signifikansi masing-masing koefisien regresi secara individu. Fungsi tidy() dari paket broom menyajikan hasil koefisien, error standar, nilai z, dan p-value untuk masing-masing prediktor dengan cara yang lebih ringkas dan rapi. Hasil uji ini membantu mengetahui prediktor mana yang berkontribusi signifikan terhadap model.

library(broom)
tidy(model_logit)

## # A tibble: 4 × 5
##   term         estimate std.error  statistic p.value
##   <chr>           <dbl>     <dbl>      <dbl>   <dbl>
## 1 (Intercept)  624.       485580.  0.00129     0.999
## 2 wt          -185.       187531. -0.000985    0.999
## 3 hp            -0.0885     2045. -0.0000433   1.00 
## 4 am           -60.4      104421. -0.000578    1.00

Odds ratio dihitung dengan mengubah koefisien regresi (dalam log odds) ke skala odds yang lebih mudah dipahami dengan menggunakan eksponensial. Selain itu, juga dihitung confidence interval 95% untuk masing-masing odds ratio dengan pendekatan Wald. Interpretasi dari odds ratio memungkinkan kita mengetahui seberapa besar peluang efisiensi bahan bakar berubah untuk setiap satu unit perubahan prediktor, sambil mempertahankan prediktor lain konstan.

# Odds Ratio dengan CI Wald
coefs <- summary(model_logit)$coefficients
OR <- exp(coefs[, "Estimate"])
CI_lower <- exp(coefs[, "Estimate"]- 1.96 * coefs[, "Std. Error"])
CI_upper <- exp(coefs[, "Estimate"] + 1.96 * coefs[, "Std. Error"])
odds_ratio_table <- cbind(OR, CI_lower, CI_upper)
round(odds_ratio_table, 3)

##                        OR CI_lower CI_upper
## (Intercept) 1.552814e+271        0      Inf
## wt           0.000000e+00        0      Inf
## hp           9.150000e-01        0      Inf
## am           0.000000e+00        0      Inf

Uji Hosmer-Lemeshow digunakan untuk menilai kecocokan model terhadap data. Dalam uji ini, prediksi model dibagi ke dalam beberapa grup (misalnya 10 grup), lalu dibandingkan antara frekuensi aktual dan prediksi di masing-masing grup. Hasil uji memberikan p-value; jika tidak signifikan (p > 0.05), maka model dianggap cukup baik dalam menjelaskan data yang diamati dan tidak ada perbedaan signifikan antara prediksi dan observasi.

library(ResourceSelection)
hoslem.test(mtcars$mpg_high, fitted(model_logit), g=10)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  mtcars$mpg_high, fitted(model_logit)
## X-squared = 1.7112e-09, df = 2, p-value = 1

Kesimpulan Model logistik menunjukkan bahwa variabel predikstor tidak memberikan pengaruh signifikan. Hal ini dibuktikan dengan p value yang besar.

4.5.4.3 Pengujian Likelihood dalam regresi logistik secara manual

Pengujian likelihood ratio dalam regresi logistik digunakan untuk membandingkan dua model: model null dan model full. Model null adalah model yang hanya mencakup intercept (tanpa variabel prediktor), sedangkan model full mencakup satu atau lebih variabel prediktor yang diduga berpengaruh terhadap variabel respons. Tujuan pengujian ini adalah untuk mengetahui apakah penambahan variabel prediktor dalam model full memberikan peningkatan signifikan terhadap kecocokan model dibanding model null. Statistik uji dihitung dengan rumus

\[ G^2 = -2(\ell_0 - \ell_1) \]

• \(\ell_0\): log-likelihood dari model null (tanpa prediktor)
  
• \(\ell_1\): log-likelihood dari model full (dengan prediktor)

Jika nilai ini lebih besar dari nilai kritis (atau p-value < 0.05), maka hipotesis nol ditolak.

Hipotesis nol (H₀) dalam pengujian ini menyatakan bahwa model null sudah cukup baik menjelaskan data, artinya variabel prediktor dalam model full tidak memberikan tambahan informasi yang signifikan. Sebaliknya, hipotesis alternatif (H₁) menyatakan bahwa model full memberikan peningkatan kecocokan yang signifikan. Jika H₀ ditolak, maka disimpulkan bahwa setidaknya satu prediktor dalam model full berpengaruh secara signifikan terhadap variabel respons, sehingga model full lebih layak digunakan daripada model null. Interpretasi dari hasil pengujian ini membantu menentukan apakah penambahan variabel dalam model logistik memang memberikan manfaat nyata dalam menjelaskan variabilitas data.

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

Untuk nilai model tersebut dapat dihitung dengan fungsi log likelihood pada persamaan diatas.

Contoh perhitungannya :

Misal kita punya data

Model null dari model ini dapat dihitung dengan cara

\[ p = \frac{1 + 0 + 1 + 0}{4} = 0.5 \]

\[ \ell_0 = \sum_{i=1}^{4} \left[ y_i \log(0.5) + (1 - y_i) \log(0.5) \right] = 4 \cdot \log(0.5) = -2.7726 \]

Kemudian misalkan hasil estimasi: \(\hat{\beta}_0 = -1\), \(\hat{\beta}_1 = 1\)

Maka hitung \(\hat{p}_i\):

\[ \hat{p}_i = \frac{1}{1 + e^{-(\hat{\beta}_0 + \hat{\beta}_1 x_i)}} \]

Dengan rumus tersebut, maka dapat diperoleh nilai loglik componentnya, tinggal dijumlahkan nilai tersebut sesuai formula pada fungsi log likelihoodnya:

\[ \ell_1 = -0.3133 + (-0.6931) + (-0.1269) + (-1.3133) = -2.4466 \]

Setelah diperoleh nilai l0 dan l1, maka dapat dihitung nilai likelihoodnya

\[ G^2 = -2(\ell_0 - \ell_1) = -2(-2.7726 + 2.4466) = 0.652 \]Nantinya nilai dari likelihood ini dibandingkan dengan chi square tabel. Jika keputusannya adalah tolak H0 maka model full lebih berpengaruh secara signifikan daripada model null.

4.5.5 Detail metode Estimasi dan Inferensi Regresi Poisson

Model ini digunakan untuk data yang bersifat cacah dan menggunakan pendekatan MLE dengan metode IRLS. Model poisson berbentuk

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

Misal sampel acak Y berdistribusi poisson dan rata rata bergantung pada vairabel bebas X, maka model regresinya akan menjadi

\[ \eta_i = \log(\lambda_i) = x_i^\top\beta \]

Untuk estimasi parameternya menggunakan pendekatan Maximum Likelihood Estimation (MLE) dengan metode Iteratively Reweighted Least Square (IRLS). Dengan model regresi sebelumya dapat dicari fungsi Log Likelihood yang dimaksimalkan

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

Serta formulasi iteratif IRLS :

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

Dengan :

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

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

dan

\[ \eta_i = \log(\lambda_i) = x_i^\top \beta \]

Simulasi data :

Sebuah perusahaan digital marketing ingin mengetahui apakah tingkat aktivitas online pengguna (x) memengaruhi jumlah klik iklan (y) yang mereka lakukan dalam satu minggu. Aktivitas online diukur berdasarkan skor yang diperoleh dari kombinasi frekuensi akses dan durasi online.

set.seed(28)
n <- 100
x <- rnorm(n)
X <- cbind(1, x) # Tambah intercept
beta_true <- c(0.5, 0.8)
eta <- X %*% beta_true
lambda <- exp(eta)
y <- rpois(n, lambda)

# Inisialisasi IRLS manual
beta <- c(0, 0)
tol <- 1e-6
max_iter <- 100
for (i in 1:max_iter) {
  eta <- X %*% beta
  lambda <- exp(eta)
  W <- diag(as.numeric(lambda))
  z <- eta + (y- lambda) / lambda
  beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
  if (sum(abs(beta_new- beta)) < tol) {
    cat("Konvergen pada iterasi ke-", i, "\n")
    break
  }
  beta <- beta_new
}

## Konvergen pada iterasi ke- 9

beta # hasil estimasi

##        [,1]
##   0.5085490
## x 0.8049463

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

## 
## Call:
## glm(formula = y ~ x, family = poisson)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.50855    0.08495   5.986 2.14e-09 ***
## x            0.80495    0.06562  12.266  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 275.55  on 99  degrees of freedom
## Residual deviance: 121.89  on 98  degrees of freedom
## AIC: 335.84
## 
## Number of Fisher Scoring iterations: 5

Metode IRLS ini memberikan cara iteratif untuk menghitung estimasi MLE nya. Dapat dilihat pada hasil syntax diatas sebagai berikut

1. Estimasi untuk lambda = B0 = 0,5, yang berarti log rata rata kejadian adalah 1.57. Selanjutnya 0.80 mengartikan setiap kenaikan 1 x menambah 0.80 y.
2. Model konvergen pada iterasi ke-9

Pengujian Hipotesis

Menggunakan uji Wald, likelihood ratio serta evaluasinya dengan AIC dan BIC

# Koefisien dan standar error
coef_val <- coef(model)[2]
se_val <- summary(model)$coefficients[2, 2]
wald_z <- coef_val / se_val
wald_chisq <- wald_z^2
p_value <- 1- pchisq(wald_chisq, df = 1)
cat("Z:", wald_z, "\nChi-Square:", wald_chisq, "\np-value:", p_value)

## Z: 4.132766 
## Chi-Square: 17.07975 
## p-value: 3.584238e-05

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

## Analysis of Deviance Table
## 
## Model 1: y ~ 1
## Model 2: y ~ x
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        99     73.111                     
## 2        98    110.721  1   -37.61

AIC(model)

## [1] 114.7213

BIC(model)

## [1] 119.9317

Dari perhitungan tersebut diperoleh beberapa hal berikut

1. Dengan uji wald dapat dilihat bahwa variabel x(tingkat aktivitas online) memengaruhi signifikan y(jumlah klik iklan)
2. Berdasarkan nilai deviencenya, model full lebih baik dari model null-nya. Hal ini dibuktikan dengan nilai devianse yang bernilai negatif(-37,61 karena deviance full lebih rendah).

4.6 Review

1. GLM merupakan metode analisis yang memungkinkan regresi linear membuat variabel respon tidak harus mengikuti distribusi normal. Bisa saja binomial, poisson, gamma dan lainnya. GLM memiliki 3 komponen utama yakni komponen acak (harus berdistribusi keluarga eksponensia), Komponen sistematik sebagai bentuk linear dan fungsi link untuk menghubungkan komponen sistematik dengan nilai ekspekstasi komponen acak.
2. Generalized Linear Model ketika data responnya berbentuk biner, maka dapat dianalisis dengan menggunakan regresi logistik. Dalam estimasi parameternya menggunakan Maksimum Like Lihood dengan pemaksimalnnya menggunakan metode Newton Raphson. Kemudian setelah ditemukannya estimasi parameter (beta), maka diuji untuk melihat apakah estimasi ini memberikan pengaruh yang signifikan. Pengujian dapat dilakukan dengan metode wald. Kemudian dilakukan uji likelihood ratio untuk membandingkan model null (tanpa predikstor) dengan model full.
3. Pada Kasus Regresi Logistik Multiple, pengujian kecocokan model dapat dmenggunakan uji hosmer lemeshow (prediktor kurang dair 10). Jika nilai p valeu lebih dari taraf signifikansi, maka model dianggap cukup baik dalam menjelaskan data.
4. Untuk variabel data frekuensi (cacah) digunakan regresi poisson. Dalam estimasi paramternya menggunakan maksimum likelihood dengna pengoptimalan menggunakan metode IRLS. Kemudian dilakukan pengujian signifikansi parameter dengan uji wald, perbandingan model dengan likelihood ratio test dan AIC

5 Regresi Logistik

5.1 Definisi Umum

Regresi logistik adalah suatu metode dalam analisis statistik yang digunakan untuk memodelkan hubungan antara satu variabel dependen (respon) yang bersifat kategorik biner, yaitu hanya memiliki dua kemungkinan hasil, seperti “ya” atau “tidak”, “sukses” atau “gagal”, “1” atau “0”, dengan satu atau lebih variabel independen (prediktor) yang dapat berskala nominal, ordinal, interval, maupun rasio. Metode ini merupakan bagian dari model regresi general linear model yang dimodifikasi untuk menangani keluaran yang tidak bersifat kontinu dan tidak mengikuti distribusi normal.

Berbeda dengan regresi linear yang menggunakan fungsi linear untuk memprediksi nilai respon, regresi logistik menggunakan fungsi logit (log-odds) untuk memodelkan probabilitas kejadian dari salah satu kategori. Dengan kata lain, regresi logistik tidak secara langsung memprediksi nilai kategorik, tetapi memprediksi probabilitas bahwa suatu peristiwa (kategori 1) akan terjadi berdasarkan nilai dari variabel-variabel prediktor yang diberikan. Probabilitas tersebut kemudian ditransformasikan melalui fungsi logistik (sigmoid) agar nilainya berada dalam rentang 0 hingga 1.

5.2 Regresi Logistik dengan berbagai skala pengukuran Variabel Prediktor

Dalam penerapannya, regresi logistik ini dapat digunakan ketika variabel predikstornya memiliki skala pengukuran yang berbeda. Skala pengukuran tersebut meliputi :

• nominal, variabel yang tidak memiliki urutan antar kategorinya, hanya sebagai pembeda dari satu variabel dengan variabel lainnya.

• ordinal, variabel yang memiliki urutan tapi jarak antar kategorinya tidak dapat dihitung

• rasio, variabel kontinu yang memiliki nol absolut dan rasio bermakna.

5.2.1 Simulasi Data

Untuk lebih jelas dilakukan simulasi sebagai berikut. Misal kita punya data 500 observasi dan akan dilihat pengaruh setiap variabel terhadap peluang sukses. Pada simulasi ini akan dibandingkan antara perlakuan ordinal ketika ditreatment dengan nominal dan ditreatment dengan numerik.

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

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

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

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

## # A tibble: 6 × 4
##   success gender education  income
##     <int> <chr>  <chr>       <dbl>
## 1       1 Male   Master       72.4
## 2       1 Male   Master       31.9
## 3       1 Female Bachelor     42.5
## 4       0 Male   Bachelor     49.3
## 5       1 Male   HighSchool   63.9
## 6       1 Female HighSchool   56.7

5.2.1.1 Eksplorasi Data

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

## # A tibble: 2 × 3
##   success Jumlah Rata2_Income
##     <int>  <int>        <dbl>
## 1       0    121         44.0
## 2       1    379         51.0

Distribusi jumlah sukses dan tidak sukses, serta rata rata pendapatan pada masing masing grup.

5.2.1.2 Treatment Variabel Ordinal

Treatment sebagai Nominal (Dummy)

Pada treatment ini variabel ordinal didefinisikan sebagai variabel nomimal sehingga variabel Highschool akan membandingkan dengan baselinenya.

sim_data_nominal <- sim_data %>%
  mutate(
    education = factor(education, levels = c("HighSchool", "Bachelor", "Master", "PhD"))
  )

model_nominal <- glm(success ~ gender + education + income, data = sim_data_nominal, family = binomial)
summary(model_nominal)

## 
## Call:
## glm(formula = success ~ gender + education + income, family = binomial, 
##     data = sim_data_nominal)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.500233   0.416402  -1.201  0.22963    
## genderMale        -0.549040   0.221714  -2.476  0.01327 *  
## educationBachelor -0.408949   0.241178  -1.696  0.08996 .  
## educationMaster    1.164894   0.391783   2.973  0.00295 ** 
## educationPhD       1.100415   0.513165   2.144  0.03200 *  
## income             0.038608   0.008266   4.671    3e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 553.37  on 499  degrees of freedom
## Residual deviance: 499.95  on 494  degrees of freedom
## AIC: 511.95
## 
## Number of Fisher Scoring iterations: 5

Interpretasi Koefisien

Intercept merupakan log-odds dasar untuk individu berjenis kelamin perempuan (Female, berpendidikan HighSchool dan income = 0. Individu ini menjadi baseline untuk interpretasi dari variabel lainnya sehingga setiap pengujian variabel lainnya akn dibandingkan dengan nilai yang ada di variabel baseline.

Pada Goodness of Fit dapat dilihat bahwa terdapat penurunan null deviance terhadap residual deviance. Hal ini menunjukkan bahwa model membawa informasi yang cukup baik. Kemudian ada nilai AIC yang mana semakin kecil nilai AIC, maka model semakin baik dan kompleks

Variabel income,gender dan education signifikan pada alpha 10% dalam meningkatkan peluang sukses .

Treatment sebagai Rasio ( Numeric )

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

## 
## Call:
## glm(formula = success ~ gender + education_numeric + income, 
##     family = binomial, data = sim_data_numeric)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.959587   0.448869  -2.138  0.03253 *  
## genderMale        -0.572377   0.218642  -2.618  0.00885 ** 
## education_numeric  0.364025   0.120454   3.022  0.00251 ** 
## income             0.036216   0.008013   4.520 6.19e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 553.37  on 499  degrees of freedom
## Residual deviance: 515.19  on 496  degrees of freedom
## AIC: 523.19
## 
## Number of Fisher Scoring iterations: 4

Dari simulasi tersebut diperoleh bahwa semua variabel memberikan pengaruh yang signifikan terhadap peluang sukses

Perbandingan Model

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

## $AIC_Nominal
## [1] 511.9455
## 
## $AIC_Numeric
## [1] 523.19

5.2.1.3 Goodness of fit

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

list(
 McFadden_R2_Nominal = r2_nominal,
 McFadden_R2_Numeric = r2_numeric
)

## $McFadden_R2_Nominal
## 'log Lik.' 0.09655104 (df=6)
## 
## $McFadden_R2_Numeric
## 'log Lik.' 0.06900278 (df=4)

Pada dua pengujian diatas dapat dilihat bahwa model dengan variabel ordinal dilakukan pendekatan dengan nominal lebih baik daripada pendekatan numerik. Hal ini dibuktikan dengan nilai AIC yang lebih rendah pada model nominal dan nilai McFadden Ryang memberikan nilai lebih besar dair pendekatan numerik.

Visualisasi Prediksi

sim_data_nominal <- sim_data_nominal %>% mutate(predicted = predict(model_nominal, type = "response"))
sim_data_numeric <- sim_data_numeric %>% mutate(predicted = predict(model_numeric, type = "response"))
 
# Plot untuk model nominal
sim_data_nominal %>%
  ggplot(aes(x = income, y =predicted, color = education)) +
  geom_point(alpha = 0.6) +
  labs(title = "Prediksi Probabilitas (Ordinal sebagai Nominal)", x= "Income", y="Prediksi Probabilitas") +
  theme_minimal()

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

Ringkasan Koefisien Model

library(knitr)
library(kableExtra)

fmt <- ifelse(knitr::is_latex_output(), "latex", "html")

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

Ringkasan Koefisien Model dengan Ordinal sebagai Nominal

term	estimate	std.error	statistic	p.value
(Intercept)	-0.5002334	0.4164019	-1.201323	0.2296258
genderMale	-0.5490396	0.2217143	-2.476338	0.0132738
educationBachelor	-0.4089492	0.2411782	-1.695631	0.0899559
educationMaster	1.1648941	0.3917827	2.973317	0.0029460
educationPhD	1.1004148	0.5131649	2.144369	0.0320033
income	0.0386077	0.0082655	4.670929	0.0000030

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

Ringkasan Koefisien Model dengan Ordinal sebagai Numeric

term	estimate	std.error	statistic	p.value
(Intercept)	-0.9595872	0.4488688	-2.137790	0.0325338
genderMale	-0.5723768	0.2186416	-2.617877	0.0088479
education_numeric	0.3640247	0.1204543	3.022098	0.0025103
income	0.0362157	0.0080130	4.519639	0.0000062

Tabel tersebut memberikan estimasi koefisien, standar error, nilai z dan p value. Koefisien postif berarti meningkatkan log odds peluang sukses sedangkan keofisien negatif menurukan.

Dengan ini diperoleh kesimpulannya bahwa

1. Variabel gender perempuan memiliki peluang sukses lebih tinggi dari laki laki karena estimate gender Malenya bernilai negatif.
2. Variabel Pendidikan, ketika dianggap nominal maka dibandingkan dengan highschool dan diperoleh bahwa pendidikan sarjana tidak memberikan efek signifikan dengan tingkat signifikansinya adalah 5%.
3. Variabel income ketika pendapatannya tinggi maka peluang suksesnya makin besar
4. Secara perhitungan, pendekatan orndinal sebagai nominal lebih baik daripada pendekatan sebagai numerik

5.3 Pemilihan Model dan Evaluasi

5.3.1 Membangun Model

Dalam regresi logistik, memilih model yang tepat merupakan langkah penting untuk memperoleh prediksi probabilitas yang akurat terhadap suatu kejadian biner. Terdapat dua pendekatan utama yang umum digunakan dalam membangun model, yaitu pendekatan Confirmatory dan Exploratory.

Confirmatory

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

Ciri-ciri utama pendekatan ini antara lain:

• Model dikembangkan berdasarkan landasan teori atau hasil penelitian terdahulu.

• Fokus utamanya adalah menguji apakah suatu pengaruh memang signifikan secara statistik, bukan semata-mata mencari model dengan kinerja terbaik.

• Peneliti biasanya memulai dengan model lengkap, kemudian mengevaluasi apakah menambahkan atau menghapus suatu efek (seperti interaksi antar variabel) memberikan perbedaan signifikan dalam kualitas model.

• Pengujian signifikansi dilakukan dengan membandingkan model yang mengandung efek tertentu dengan model yang tidak, misalnya menggunakan Likelihood Ratio Test.

Contoh penerapan:
Jika suatu teori menyatakan bahwa variabel x1 dan x2 memengaruhi kemungkinan seseorang membeli produk, maka model logistik langsung dibangun dengan kedua variabel tersebut. Selanjutnya, diuji apakah keberadaan x2 dalam model benar-benar memberikan kontribusi yang signifikan.

Explanatory

Pendekatan Eksploratori digunakan ketika peneliti belum memiliki teori yang jelas atau ingin menggali kemungkinan hubungan antar variabel yang belum diketahui.

Ciri-ciri pendekatan ini meliputi:

• Model dikembangkan secara bertahap dengan tujuan menemukan kombinasi variabel prediktor yang optimal.

• Pemilihan variabel didasarkan pada pertimbangan statistik seperti AIC, deviance, atau log-likelihood.

Proses seleksi variabel dapat dilakukan dengan beberapa metode, yakni Forward Selection dengan dimulai dari model kosong, lalu variabel ditambahkan satu per satu jika terbukti signifikan. Backward Elimination dengan dimulai dari model lengkap, kemudian variabel yang tidak signifikan dihapus. Stepwise Selection adalah gabungan keduanya, di mana variabel dapat masuk atau keluar model secara dinamis berdasarkan kriteria statistik.

Tujuan utama dari pendekatan ini adalah memperoleh model yang parsimonious, yaitu model yang cukup sederhana namun tetap memiliki kemampuan prediksi yang baik.

Pemilihan antara pendekatan Confirmatory dan Exploratory bergantung pada tujuan analisis: Jika tujuannya untuk menguji hipotesis tertentu, gunakan pendekatan Confirmatory.Jika tujuannya untuk menemukan model terbaik berdasarkan pola dalam data, gunakan pendekatan Exploratory.

Dalam praktiknya, kedua pendekatan ini sering saling melengkapi: teori digunakan sebagai fondasi awal, kemudian pendekatan eksploratori diterapkan untuk menyempurnakan struktur model.

5.3.2 Simulasi Perhitungan

5.3.2.1 Input library dan data

Seseorang ingin meneliti apakah karakteristik individu seperti tingkat stres (x1), status merokok (x2), dan aktivitas fisik (x3) memengaruhi kemungkinan terkena insomnia (y). Data disimulasikan menggunakan model regresi logistik.

library(knitr)
library(dplyr)
library(ggplot2)
library(MASS)
library(caret)
library(pROC)
library(DescTools)
set.seed(28)

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

##   y          x1 x2         x3
## 1 0 -1.90215722  0 -0.8897054
## 2 0 -0.06429479  0 -2.1269942
## 3 0 -1.33116707  0 -1.3989726
## 4 0 -1.81999167  1 -0.5385668
## 5 0  0.16266969  1  1.3377003
## 6 0  0.53139634  1  0.4151210

5.3.2.2 Pembentukan Model

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

## 
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.5749     0.2461  -2.336  0.01949 *  
## x1            1.3974     0.2366   5.907 3.49e-09 ***
## x2           -0.7872     0.3727  -2.112  0.03467 *  
## x3            0.4651     0.1787   2.603  0.00925 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 244.35  on 199  degrees of freedom
## Residual deviance: 185.01  on 196  degrees of freedom
## AIC: 193.01
## 
## Number of Fisher Scoring iterations: 5

Dari hasil tersebut dapat dilihat bahwa semua variabel X memberikan pengaruh yang signifikan .

5.3.2.3 Metode Stepwise

null_model <- glm(y ~ 1, data = df, family = binomial)

step_forward <- step(null_model, direction = "forward", scope = formula(model_full), trace = FALSE)

step_backward <- step(model_full, direction = "backward", trace = FALSE)

step_both <- step(null_model, direction = "both", scope = formula(model_full), trace = FALSE)

AIC(model_full, step_forward, step_backward, step_both)

##               df      AIC
## model_full     4 193.0117
## step_forward   4 193.0117
## step_backward  4 193.0117
## step_both      4 193.0117

Diperoleh nilai AIC yang sama untuk setiap metode seleksi. Untuk itu, analisis selanjutnya boleh memakaimetode manapun. Namun disini kita akan mencoba untuk menggunakan metode stepwise selection.

5.3.2.4 Evaluasi Model : ROC dan AUC

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

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

auc(roc_obj)

## Area under the curve: 0.8225

5.3.2.5 Pseudo R Squared

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

##   CoxSnell Nagelkerke   McFadden 
##  0.2567110  0.3639846  0.2428283

5.3.2.6 Tabel Klasifikasi dan Evaluasi

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

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 123  33
##          1  17  27
##                                          
##                Accuracy : 0.75           
##                  95% CI : (0.684, 0.8084)
##     No Information Rate : 0.7            
##     P-Value [Acc > NIR] : 0.06955        
##                                          
##                   Kappa : 0.3557         
##                                          
##  Mcnemar's Test P-Value : 0.03389        
##                                          
##             Sensitivity : 0.4500         
##             Specificity : 0.8786         
##          Pos Pred Value : 0.6136         
##          Neg Pred Value : 0.7885         
##              Prevalence : 0.3000         
##          Detection Rate : 0.1350         
##    Detection Prevalence : 0.2200         
##       Balanced Accuracy : 0.6643         
##                                          
##        'Positive' Class : 1              
## 

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

## Sensitivity Specificity 
##   0.4500000   0.8785714

Pada seluruh pengujian diatas, diperoleh bahwa metode seleksi terbaik adalah stepwise selection. Kemudian modenya dapat memprediksi dengan akurasi 75% (cukup baik) Nilai R square psuedo cukup baik juga. Setiap variabel x memberikan pengaruh yang signifikan Kemudian pada klasifikasinya cukup baik dengan 123 posisif ditebak benar, 27 negatif ditebak benar, 33 positif malah ditebak negatif dan 17 positif, malah dianggap negatif yang mana positif adalah mengalami insomnia sedangkan negatif adalah tidak mengalami.

5.3.3 Metode Perbandingan Model dalam Regresi Logistik

Pada bagian ini akan dibahas mengenai cara membandingkan model regresi logistik menggunakan ukuran yakni Deviance AIC (Akaike Information Criterion), Likelihood-Ratio serta menjelaskan prinsip Parsimony dalam pemilihan model

5.3.3.1 Simulasi dan input data

Seorang peneliti ingin mengetahui apakah faktor usia (x1), jenis kelamin (x2), dan tingkat aktivitas harian (x3) berpengaruh terhadap kemungkinan seseorang mengalami tekanan darah tinggi (y). Data disimulasikan menggunakan model regresi logistik dengan 300 observasi.

library(MASS)
library(broom)
library(DescTools)
set.seed(28)

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

5.3.3.2 Pembuatan Model

model1 <- glm(y ~ x1, data = data, family = binomial)
model2 <- glm(y ~ x1 + x2, data = data, family = binomial)
model3 <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)

5.3.3.3 Perbandingan AIC dan Ceviance

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

##     Model      AIC Deviance
## 1 Model 1 299.0623 295.0623
## 2 Model 2 295.6804 289.6804
## 3 Model 3 268.8393 260.8393

Model 3 memiliki nilai AIC dan deviance paling rendah dibandingkan Model 1 dan Model 2, menandakan bahwa ia merupakan model dengan fit terbaik dan kompleksitas paling efisien. Penurunan AIC yang cukup besar dari Model 2 ke Model 3 menunjukkan peningkatan signifikan dalam kinerja model. Oleh karena itu, Model 3 lebih disarankan untuk digunakan dalam analisis lebih lanjut.

5.3.3.4 Likelihood Ratio Test

anova(model1, model2, test = "LRT")

## Analysis of Deviance Table
## 
## Model 1: y ~ x1
## Model 2: y ~ x1 + x2
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       298     295.06                       
## 2       297     289.68  1   5.3819  0.02035 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(model2, model3, test = "LRT")

## Analysis of Deviance Table
## 
## Model 1: y ~ x1 + x2
## Model 2: y ~ x1 + x2 + x3
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       297     289.68                          
## 2       296     260.84  1   28.841 7.857e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dengan likelihood ratio test ini dapat dilihat pengaruh dari variabelnya. Ketika membandingkan model 1 dan 2 yang mana model 2 menambahkan variabel x2, model mampu menjelaskan data dengan lebih baik yang berarti bahwa model 2 memiliki pengaruh signifikan dalam pemodelan.

Selanjutnya ketika uji model 2 dan 3 dengan perbandingkan menambah variabel X3, juga memberikan pengaruh yang signifikan sehingga model 3 menjadi model yang paling cocok untu kmenjealskan data.

5.3.3.5 Prinsip Parsimony

Model yang lebih kompleks cenderung menghasilkan nilai AIC dan deviance yang lebih rendah. Namun, model yang lebih sederhana sering kali lebih mudah untuk diinterpretasikan. Jika penurunan nilai AIC tidak terlalu signifikan, maka lebih disarankan memilih model yang lebih sederhana. Prinsip parsimony penting dalam pemodelan karena dapat mencegah terjadinya overfitting terhadap data.

Rumus AIC dinyatakan sebagai:

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

di mana log𝐿 adalah nilai log-likelihood dari model, dan 𝑘 adalah jumlah parameter dalam model. AIC digunakan untuk mengevaluasi keseimbangan antara kesesuaian model (melalui log-likelihood) dan kompleksitas model (melalui jumlah parameter). AIC yang lebih kecil menunjukkan model yang lebih baik secara keseluruhan, karena AIC memberikan penalti terhadap model yang terlalu kompleks meskipun memiliki likelihood tinggi.

Deviance didefinisikan dengan rumus:

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

Deviance mengukur sejauh mana model saat ini menyimpang dari model saturasi (model yang cocok sempurna dengan data). Semakin kecil nilai deviance, semakin baik model dalam memprediksi data aktual.

Sementara itu, statistik Likelihood Ratio (G²) dihitung dengan rumus: \[ G^2 = -2 \left( \log L_0 - \log L_1 \right) \] di mana log𝐿₀ adalah log-likelihood model yang lebih sederhana, dan log𝐿₁ adalah log-likelihood model yang lebih kompleks. G² digunakan untuk menguji apakah penambahan variabel ke dalam model memberikan peningkatan signifikan terhadap kecocokan model. Jika nilai G² cukup besar dan p-value yang menyertainya kecil, maka model yang lebih kompleks dianggap lebih unggul secara statistik.

5.3.3.6 Evaluasi Tabel Klasifikasi dan Akurasi Model

pred_prob <- predict(model3, type = "response")
pred_class <- factor(ifelse(pred_prob >= 0.5, 1, 0))
conf_matrix <- confusionMatrix(pred_class, data$y, positive = "1")
conf_matrix

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 202  41
##          1  16  41
##                                          
##                Accuracy : 0.81           
##                  95% CI : (0.761, 0.8528)
##     No Information Rate : 0.7267         
##     P-Value [Acc > NIR] : 0.0005172      
##                                          
##                   Kappa : 0.4714         
##                                          
##  Mcnemar's Test P-Value : 0.0014785      
##                                          
##             Sensitivity : 0.5000         
##             Specificity : 0.9266         
##          Pos Pred Value : 0.7193         
##          Neg Pred Value : 0.8313         
##              Prevalence : 0.2733         
##          Detection Rate : 0.1367         
##    Detection Prevalence : 0.1900         
##       Balanced Accuracy : 0.7133         
##                                          
##        'Positive' Class : 1              
## 

Dari output tersebut dapat dilihat akurasinya sebesar 81% yang menyatakan ini sangat baik.

5.3.3.7 Sensitivitas dan Spesifisitas

Sensitivitas adalah kemampuan model untuk mendekteksi kelas positif dengan tepat

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

Spesifisitas adalah kemampuan odel untuk mendeteksi kelas negatif dengan tepat.

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

Dengan : TP adalah jumlah kasus positif yang diprediksi benar oleh model FN adalah jumlah kasus positif yang salah diprediksi oleh model TN adalah jumlah kasus negatif yang diprediksi benar oleh model FP adalah jumlah kasus negatif yang salah diprediksi oleh model

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

## Sensitivity Specificity 
##   0.5000000   0.9266055

Hasilnya menyatakan bahwa model masih cukup baik untuk mendeteksi kelas positif namun sangat baik mendekateksi kelas negatif.

Dengan ini dapat dilihat beberapa metode evaluasi model yakni

• Deviance yang semakin kecil berarti semakin baik kecocokan modelnya,

• AIC yang rendah memmberikan keseimbangan antara kecocokan dan kompleksitas

• Likelihood reasio test untuk mengevaluasi apakah model kompleks signifikan,

• tabel klasifikasi membantuk menilai kinerja prediksi aktual dan prediksi model, serta

• prinsip parsimony mengutamakan model sederhana jika performanya mirip

5.3.4 Detail ROC

Kurva ROC adalah alat visual yang digunakan untuk mengevaluasi performa klasifikasi data biner. Kurva ini cocok untuk data yang proporsi kelasnya seimbang. Kurva ini dianggap ideal jika memiliki bentuk

1. Naik tajam secara vertiakl (menuju sensitivitas = 1)
2. Bergerak horizontal menuju 1 sprecifity = 1
3. Area dibawah kurva (AUC/Concordance index) mendekati 1 (Ketika bernilai 0,5 maka dianggap tidak lebih baik dari tebak acak, ketika lebih dari 0,7 maka dianggap cukup baik).

Kurva ROC ini memmberikan perbandingan peforma beberapa model klasifikasi serta memberikan pilihan optimal berdasarkan kebutuuhan aplikasi. Misal lebih baik menghindari false negatif atau positif?

Contoh

Seorang peneliti ingin mengkaji apakah variabel usia (x1), jenis kelamin (x2), dan pola tidur (x3) dapat memprediksi kemungkinan seseorang mengalami gangguan kecemasan (y). Data disimulasikan menggunakan model regresi logistik, dan performa model dievaluasi menggunakan kurva ROC. Hasil area under the curve (AUC) menunjukkan seberapa baik model membedakan antara individu yang mengalami kecemasan dan yang tidak.

library(pROC)
set.seed(28)
x1 <- rnorm(200)
x2 <- rbinom(200, 1, 0.5)
x3 <- rnorm(200)
lin_pred <--1 + 1.5 * x1- 0.7 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(200, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
pred <- predict(model, type = "response") 
roc_obj <- roc(data$y, pred)
plot(roc_obj)

auc(roc_obj)

## Area under the curve: 0.8718

Simulasi Pemilihan Threshold Optimal

Untuk memilih threshold terbaik, kita bisa mengevaluasi sensitivitas dan spesifisitas pada berbagai cut-off.

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

results$Sensitivity <- sapply(thresholds, function(t) {
 pred_class <- ifelse(pred >= t, 1, 0)
 cm <- table(Pred = pred_class, Obs = data$y)
 TP <- cm["1", "1"]
 FN <- cm["0", "1"]
 TP / (TP + FN)
})

results$Specificity <- sapply(thresholds, function(t) {
 pred_class <- ifelse(pred >= t, 1, 0)
 cm <- table(Pred = pred_class, Obs = data$y)
 TN <- cm["0", "0"]
 FP <- cm["1", "0"]
 TN / (TN + FP)
})

print(results)

##    Threshold Sensitivity Specificity
## 1       0.10   0.9615385   0.5405405
## 2       0.15   0.9423077   0.6486486
## 3       0.20   0.8653846   0.6959459
## 4       0.25   0.8076923   0.7500000
## 5       0.30   0.7307692   0.8108108
## 6       0.35   0.6538462   0.8445946
## 7       0.40   0.6346154   0.8851351
## 8       0.45   0.5576923   0.9189189
## 9       0.50   0.5576923   0.9256757
## 10      0.55   0.5000000   0.9324324
## 11      0.60   0.4615385   0.9459459
## 12      0.65   0.4038462   0.9527027
## 13      0.70   0.3076923   0.9797297
## 14      0.75   0.2884615   0.9864865
## 15      0.80   0.2115385   0.9864865
## 16      0.85   0.1730769   0.9932432
## 17      0.90   0.1153846   1.0000000

Cut-off optimal bisa dipilih berdasarkan maksimum dari Sensitivity + Specificity Atau mempertimbangkan trade-off sesuai tujuan aplikasi (misalnya: jika False Negative harus dihindari, maka prioritaskan sensitivitas tinggi).

5.3.5 Precision-Recall Curve (PR Curve)

Kurva Precision-Recall (PR) adalah alat evaluasi performa model klasifikasi, khususnya sangat berguna saat bekerja dengan data yang tidak seimbang (class imbalance). Ada dua komponen yakni

Precision (Presisi), yakni proporsi hasil prediksi positif memang benar bernilai positif.

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

Recall (Sensitivitas), yakni proporsi berhasil memprediksi positif ketika data tersebut memang positif

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

PR Curve ini menunjukkan bagaimana presisi berubah saat recall meningkat. Ideealinya kita ingin kedua nilai ini tinggi namun karena ada trade off maka model yang dianggap baik itu adalah model yang PR Curvenya melengkung ke pojok kanan atas. Luas kurva mendekati satu berarti model sangat baik.

contoh

Seorang peneliti ingin mengkaji apakah variabel usia (x1), jenis kelamin (x2), dan pola tidur (x3) dapat memprediksi kemungkinan seseorang mengalami gangguan kecemasan (y). Data disimulasikan menggunakan model regresi logistik, dan performa model dievaluasi menggunakan PR Curve.

library(PRROC)
set.seed(28)
x1 <- rnorm(200)
x2 <- rbinom(200, 1, 0.5)
x3 <- rnorm(200)
lin_pred <--1 + 1.5 * x1- 0.7 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(200, 1, p)
data <- data.frame(y = y, x1, x2, x3)

model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
prob <- predict(model, type = "response")

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

PR Curve cocok untuk analisis data yang kelas positifnya jauh lebih jarang daripada kelas negatif seperti deteksi penipuan dan diagnosis penyakit langka. Dibanding dengan ROC Curve, PR Curve ini fokus pada kelas positif saja.

5.3.6 Pseudo R-squared pada Regresi Logistik

Pseudo R-squared pada regresi logistik bertujuan untuk memberikan ukuran seberapa baik model menjelaskan variasi dalam data, sebagai analog dari R-squared pada regresi linear. Karena regresi logistik memodelkan probabilitas kejadian dari variabel biner, maka R-squared konvensional tidak dapat digunakan, sehingga dikembangkan beberapa versi Pseudo R-squared seperti McFadden, Cox & Snell, dan Nagelkerke. Nilai ini digunakan untuk menilai kualitas model, membandingkan beberapa model, serta mengevaluasi peningkatan performa saat menambahkan variabel baru. Meskipun tidak memiliki interpretasi sekuat R-squared linear, Pseudo R-squared tetap berguna dalam menilai kecocokan model logistik secara umum.

Contoh

Seorang peneliti ingin mengkaji apakah variabel usia (x1), jenis kelamin (x2), dan pola tidur (x3) dapat memprediksi kemungkinan seseorang mengalami gangguan kecemasan (y). Data disimulasikan menggunakan model regresi logistik, dan performa adalah menggunakan Psuedo R Squared

Simulasi

set.seed(28)
n <- 300
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)
lin_pred <--1 + 1.2 * x1- 0.6 * x2 + 0.8 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)

model logistik dan null

model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
model_null <- glm(y ~ 1, data = data, family = binomial)

Like lihood dan Rumus

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

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

dengan L0 adalah likelihood model null (tanpa prediktor ) dan LM adlaah likelihood penuh.

Perhitungan Manual R Square

logL0 <- logLik(model_null)
logLM <- logLik(model)
L0 <- exp(logL0)
LM <- exp(logLM)
 
n <- nobs(model)
 
cox_snell <- 1- (L0 / LM)^(2 / n)
mcfadden <- 1- (as.numeric(logLM) / as.numeric(logL0))
 
r2 <- data.frame(
 R2_Cox_Snell = cox_snell,
 R2_McFadden = mcfadden
)
r2

##   R2_Cox_Snell R2_McFadden
## 1    0.2618645    0.258827

Perhitungan menggunakan package

Ketika menggunakan Package, diperoleh beberapa nilai pseudonya. Package yang dapat digunakan ada package rcompanion, DecsTools, dan pscl.

Menggunakan pscl

library(pscl)
pR2(model)

## fitting null model for pseudo-r2

##          llh      llhNull           G2     McFadden         r2ML         r2CU 
## -130.4196517 -175.9638404   91.0883775    0.2588270    0.2618645    0.3791889

Menggunakan rcompanion

library(rcompanion)
nagelkerke(model)

## $Models
##                                               
## Model: "glm, y ~ x1 + x2 + x3, binomial, data"
## Null:  "glm, y ~ 1, binomial, data"           
## 
## $Pseudo.R.squared.for.model.vs.null
##                              Pseudo.R.squared
## McFadden                             0.258827
## Cox and Snell (ML)                   0.261865
## Nagelkerke (Cragg and Uhler)         0.379189
## 
## $Likelihood.ratio.test
##  Df.diff LogLik.diff  Chisq    p.value
##       -3     -45.544 91.088 1.2787e-19
## 
## $Number.of.observations
##           
## Model: 300
## Null:  300
## 
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
## 
## $Warnings
## [1] "None"

Menggunakan DecsTools

library(DescTools)
PseudoR2(model, which = "all")

##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.2588270       0.2360950       0.2618645       0.3791889       0.2329100 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.4314536       0.3017814       0.4586520       0.2954017     268.8393033 
##             BIC          logLik         logLik0              G2 
##     283.6544332    -130.4196517    -175.9638404      91.0883775

Dengan ini diperoleh beberapa nilai R Square, namun biasanya digunakan McFadden R Square dan Cox & Snell R Square. Umumnya nilai R Square mendekati 1 berarti model memiliki kekuatan prediktif yang baik. Namun menurut McFadden (1977), nilai ini dianggap sebagai model dengan kecocokan baik ketika nilainya lebih dari 0,2. Sedangkan nilai Cox & Snell tidak memiliki kriteria khusus sehingga nilai ini lebih cocok untuk membandingkan satu model dengna model lainnya.

5.4 Review

1. Regresi Logistik adalah suatu metode analisis yang digunakan untuk memodelkan hubungan variabel yang bersifat kategori, ketika variabel responnya adalah biner. Model ini tidak memprediksi nilai numerik, tapi memprediksi peluang bahwa suatu peristiwa termasuk ke dalam kelompok mana.
2. Regresi Logistik dapat digunakan ketika variabel prediktor memiliki skala pengukuran yang berbeda. Pada model ini, kamu bisa menggabungkan variabel prediktor nominal dengan numerik, bahkan bisa melakukan pendekatan kepada variabel ordinal. Perndekatan tersebut meliputi pendekatan ordinal sebagai nominal dan ordinal sebagai numerik. Nantinya, model ini dibandingkan untuk melihat pendekatan mana yang cocok menggunakan AIC (semakin kecil semakin baik), atau McFadden (semakin besar semakin baik).
3. Dalam membangun model dapat menggunakan 2 pendekatan, yakni confirmatory dan explanatory. Pendekatan ini digunakan berdasarkan tujuan penelitian. Jika tujuannya adalah menguji hipotesis tertentu, maka pendekatan confirmaoty lebih cock namun ika ingin melihat model mana yang terbaik berdasarkan pola data, maka lebih baik explanatory.
4. Dalam menyeleksi variabel prediktor bisa digunakan dengan metode stepwise baik itu forward, bakcward maupun gabungannya.
5. Untuk membandingkan model dapat digunakan ukuran Deviance AIC yang mana semakin kecil nilai AIC maka semakin baik modelnya. Kemudian dapat menggunakanLikelihood test untuk melihat vairabel yang signifikan ketika dimasukkan ke model serta prinsip parsimony. Dengan pengujian ini dapat dilihat juga kemampuan model untuk mendeteksi suatu data (sensitivitas dan spesifisitas)
6. Untuk mengevaluias performa klasifikasi data biner, maka dapat dilakukan dengan ROC. Namun jika datanya tidak seimabng, maka lebih baik dievaluiasi dengan PR Curve meskipun hanya fokus pada nilai positif saja. Selain itu untuk melihat seberapa baik model menjelaskan variasi data, dapat digunakan nilai R Square yang dapat dihitung dengan banyak metode. NKita disini menggunakan McFadden (lebih dari 0,2 sudah baik) dan Cox Snell yang cocok untuk perbandingan model juga.

6 Distribusi Multinomial

6.1 Definisi Umum

Distribusi multinomial adalah generalisasi dari distribusi binomial untuk lebih dari dua kategori. Distribusi ini menggambarkan probabilitas hasil dari sejumlah percobaan independen yang masing-masing menghasilkan satu dari k kemungkinan kategori, dengan probabilitas tetap untuk setiap kategori dan total percobaan sebanyak n. Jika X₁, X₂, …, Xₖ adalah jumlah kejadian dalam masing-masing kategori dan p₁, p₂, …, pₖ adalah probabilitas masing-masing kategori (dengan ∑pᵢ = 1), maka vektor (X₁, X₂, …, Xₖ) mengikuti distribusi multinomial dengan parameter n dan p.

\[ P(X_1 = x_1, \ldots, X_k = x_k) = \frac{n!}{x_1! \cdot x_2! \cdots x_k!} \cdot p_1^{x_1} \cdot p_2^{x_2} \cdots p_k^{x_k}, \quad\\ \text{dengan } \sum_{i=1}^{k} x_i = n, \sum_{i=1}^{k} p_i = 1 \]

6.1.1 Contoh Kasus dan Perhitungan Manual

Sebuah survei dilakukan terhadap 10 mahasiswa tentang minuman favorit mereka: kopi, teh, atau susu. Diketahui bahwa secara umum 50% mahasiswa menyukai kopi, 30% teh, dan 20% susu. Hitung probabilitas bahwa dari 10 mahasiswa, 5 memilih kopi, 3 teh, dan 2 susu.

Langkah pertama : Deklarasikan hal yang diketahui

\[ \text{Diketahui } n = 10, \; p_1 = 0.5, \; p_2 = 0.3, \; p_3 = 0.2 \\ \text{Ingin dicari: } P(X_1 = 5, X_2 = 3, X_3 = 2) \]

Langkah Kedua : Gunakan rumus dan Hitung

\[ P(5,3,2) = \frac{10!}{5! \cdot 3! \cdot 2!} \cdot (0.5)^5 \cdot (0.3)^3 \cdot (0.2)^2 \]

\[ = \frac{3628800}{120 \cdot 6 \cdot 2} \cdot 0.03125 \cdot 0.027 \cdot 0.04 \]

\[ = 2520 \cdot 0.00003375 = 0.08505 \]

Langkah Ketiga : Kesimpulan

Dengan perhitungan diperolh bahwa peluang 10 mahasiswa tersebut memperoleh 5 kopi, 3 teh dan 2 susu adalah 8,5%.

6.1.2 Contoh kasus dan Perhitungan R

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

## [1] 0.08505

Dengan ini diperleh bahwa pelaung mahasiswa memiilih 5 kopi 3 teh dan 2 susu adalah 8,5%.

6.2 Multinomial Logistic Regression

Multinomial logistic regression adalah metode statistik yang digunakan untuk memodelkan hubungan antara satu variabel kategorik nominal dengan lebih dari dua kategori dan satu atau lebih variabel independen. Berbeda dengan regresi logistik biner yang hanya membedakan dua kelas, multinomial logit digunakan ketika variabel dependen tidak memiliki urutan alami antar kategorinya, misalnya dalam kasus memilih jenis transportasi (mobil, motor, jalan kaki). Model ini sangat berguna ketika kita ingin memahami bagaimana variabel-variabel prediktor mempengaruhi probabilitas suatu observasi termasuk ke dalam salah satu dari beberapa kategori. \[ \textbf{Model Logit untuk Kategori ke-}j: \\ \log\left(\frac{P(Y = j)}{P(Y = k)}\right) = \beta_{j0} + \beta_{j1}X_1 + \beta_{j2}X_2 + \dots + \beta_{jp}X_p, \quad\\ \text{untuk } j = 1, 2, \dots, k-1 \]

\[ \textbf{Fungsi Probabilitas Kategori ke-}j: \\ P(Y = j) = \frac{e^{\eta_j}}{\sum_{l=1}^{k} e^{\eta_l}}, \quad\\ \\ \text{dengan } \eta_j = \beta_{j0} + \beta_{j1}X_1 + \dots + \beta_{jp}X_p \]

6.2.1 Baseline category logit model

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

\[ \log\left(\frac{\pi(j)}{\pi(c)}\right), \quad \text{untuk } j = 1, \dots, c-1 \]

Dengan pi(j) adalah probabilitas respon berada dalam kategori j sedangkan pi(c) adalah probabilitas respon berada di kategori acuan. Maka, terdapat sebanyak (𝑐 − 1) fungsi logit.

Catatan: Kategori baseline bisa ditentukan secara eksplisit, tetapi default di R adalah kategori terakhir.

Model Regresi Jika terdapat satu prediktor 𝑥, maka bentuk umum model logit-nya adalah:

\[ \log\left(\frac{\pi(j)}{\pi(c)}\right) = \alpha_j + \beta_j x \]

Contoh Kasus: 3 Kategori Respon Misalkan respon 𝑌 memiliki tiga kategori: 𝑌 ∈ {1,2,3}, dan kita gunakan kategori ke-3 sebagai baseline. Maka

\[ \log\left(\frac{\pi_1}{\pi_3}\right) = \alpha_1 + \beta_1 x \] \[ \log\left(\frac{\pi_2}{\pi_3}\right) = \alpha_2 + \beta_2 x \]

Terdapat dua model logit, satu untuk perbandingan kategori 1 dengan 3, dan satu lagi untuk kategori 2 dengan 3. Relasi Antar Kategori Jika ingin menghitung logit antara kategori 1 dan 2

\[ \log\left(\frac{\pi_1/\pi_3}{\pi_2/\pi_3}\right) = \log\left(\frac{\pi_1}{\pi_3}\right) - \log\left(\frac{\pi_2}{\pi_3}\right) \] \[ (\alpha_1 + \beta_1 x) - (\alpha_2 + \beta_2 x) = (\alpha_1 - \alpha_2) + (\beta_1 - \beta_2)x \] Model baseline-category logit digunakan unruk respon kategori yang lebih dari dua. sehingga menghasilkan fungsi logit terhadap satu baseline.

6.2.2 Estimasi Parameter

Estimasi dilakukan menggunakan metode maximum likelihood dengan algoritma iteratif seperti Newton Raphson. Log-likelihood:

\[ \ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} \sum_{j=1}^{K} y_{ij} \log(\pi_{ij}) \]

6.2.3 Contoh Kasus dan Simulasi Perhitungan

6.2.3.1 Contoh kasus 1

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

• Device: Jenis perangkat yang dipilih (Laptop, Desktop, Tablet)

• Age: Usia karyawan

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

• Experience: Tahun pengalaman kerja

Tujuannya adalah: Mengetahui bagaimana usia, divisi, dan pengalaman kerja memengaruhi pilihan perangkat kerja.

Simulasi

set.seed(28)
n <- 150

Department <- sample(c("IT", "Marketing", "HR"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 35, sd = 7))
Experience <- round(pmax(rnorm(n, mean = 7, sd = 3), 0))
 
# Simulasikan Device berdasarkan probabilitas berbeda per Department
Device <- sapply(Department, function(dep) {
 if (dep == "IT") {
 sample(c("Laptop", "Desktop", "Tablet"), size = 1, prob = c(0.6, 0.2, 0.2))
 } else if (dep == "Marketing") {
 sample(c("Laptop", "Desktop", "Tablet"), size = 1, prob = c(0.3, 0.3, 0.4))
 } else {
 sample(c("Laptop", "Desktop", "Tablet"), size = 1, prob = c(0.4, 0.4, 0.2))
 }
})

df <- data.frame(Device = factor(Device), Age, Department = factor(Department), Experience)
df$Device <- relevel(df$Device, ref = "Laptop") # baseline
head(df)

##    Device Age Department Experience
## 1  Laptop  31         IT          8
## 2  Laptop  36         IT          8
## 3  Laptop  37         IT          8
## 4  Laptop  34         IT         10
## 5 Desktop  21  Marketing         11
## 6  Tablet  31  Marketing          9

Estimasi Model

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

## # weights:  18 (10 variable)
## initial  value 164.791843 
## iter  10 value 141.736231
## final  value 141.098538 
## converged

summary(model_mnlogit)

## Call:
## multinom(formula = Device ~ Age + Department + Experience, data = df)
## 
## Coefficients:
##         (Intercept)         Age DepartmentIT DepartmentMarketing Experience
## Desktop    1.301547 -0.03860602   -2.2684307         -0.08200637 0.02670575
## Tablet    -2.292099  0.01684437   -0.5505138          1.34794203 0.10732146
## 
## Std. Errors:
##         (Intercept)        Age DepartmentIT DepartmentMarketing Experience
## Desktop    1.111799 0.02879988    0.6118749           0.4667083 0.06984485
## Tablet     1.305083 0.03064649    0.5930516           0.5421861 0.07647548
## 
## Residual Deviance: 282.1971 
## AIC: 302.1971

Nilai P-Value dan Interpretasi

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

##         (Intercept)    Age DepartmentIT DepartmentMarketing Experience
## Desktop      0.2417 0.1801       0.0002              0.8605     0.7022
## Tablet       0.0790 0.5826       0.3533              0.0129     0.1605

Interpretasi

Koefisien untuk kategori “Desktop” dan “Tablet” dibandingkan dengan baseline “Laptop”. Jika Nilai p-value kecil (< 0.05) menunjukkan variabel tersebut signifikan memengaruhi preferensi perangkat.

Prediksi dan Validasi

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

##          Actual
## Predicted Laptop Desktop Tablet
##   Laptop      40       9     13
##   Desktop     18      27     12
##   Tablet      10      10     11

Kesimpulan

Model regresi logistik multinomial tersebut menjelaskan hasil signifikansi dari setiap variabel serta tabel klasifikasi prediksinya. Dari tabel dapat dilihat bahwa 40 laptop berhasil dianggap sebagai laptop, namun 18 lainnya dianggap seagai desktop dan 10 lainnya dianggap tablet.

6.2.3.2 Contoh kasus kedua

Pada kasus ini dicoba menggunakan set data yang ada di R, yakni data iris.

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

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

summary(model)

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

Interpretasi koefisien

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

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

Disini dapat dilihat bahwa p valuenya lebih dari taraf signifikansi, berarti semua variabel tidak memberikan pengaruh yang signifikan sehingga ini menjelaskan bahwa variabel tersebut tidak menjelaskan pengaruhnya terhadap probabilitas kelas.

Prediksi dan Visualisasi

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

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

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

Meskipun p value tadi memberikan keputusan yang tidak signifikan, namun pada pengklasifikasiannya diperoleh hasil yang baik. ini bisa saja terjadi mungkin karena adanya asumsi lain seperti multikolinearitas dan lainnya.

Oleh karena itu, jika tujuannya adalah mengklasifikasikan sebuah data, maka model ini sudah cocok. Tapi jika tujuannya adalah untuk membuat model, maka ini belum menjadi model yang baik.

6.3 Review

Distribusi Multinomial adalah distribusi kategori yang variabel responnya lebih dari 2 pilihan.

Multinomial logistic regression efektif digunakan untuk klasifikasi kategori nominal. Model ini memberikan estimasi log-odds terhadap baseline dan prediksi kategori baru berdasarkan kombinasi prediktor.

Proses analisisnya dimulai dengan input data dan estimasi parameter. Kemudian menggunakan statistik z untuk melihat signifikansi parameter yang diestimasi tadi. Apakah variabel memberikan pengaruh terhadap model? Terakhir diperiksa hasil prediksi menggunakan confus matrix (tabel pengklasifikasian) untuk melihat hasil klasifikasi data.

7 Regresi Logistik Ordinal

7.1 Definisi Umum

Regresi logistik ordinal adalah metode statistik yang digunakan untuk memodelkan hubungan antara variabel respon kategorik yang memiliki urutan (ordinal) dengan satu atau lebih variabel prediktor. Tidak seperti regresi logistik multinomial yang tidak mempertimbangkan urutan kategori, regresi logistik ordinal memanfaatkan informasi tentang peringkat atau tingkat kategori respon, misalnya tingkat kepuasan dari “sangat tidak puas” hingga “sangat puas”.

Model yang paling umum digunakan adalah cumulative logit model dengan asumsi proportional odds, di mana logit dari probabilitas kumulatif untuk setiap ambang batas kategori dimodelkan sebagai fungsi linear dari prediktor. Asumsi ini menyatakan bahwa pengaruh setiap prediktor adalah konstan di seluruh ambang kategori, yang menjadikan model lebih sederhana dan efisien dibandingkan model tanpa asumsi tersebut. Regresi logistik ordinal cocok digunakan ketika tujuan analisis adalah untuk memahami atau memprediksi peluang relatif suatu observasi berada di atau di bawah suatu tingkat kategori.

7.2 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_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p, \quad \\ \text{untuk } j = 1, 2, \dots, J - 1 \]

Interpretasi Keofesien

Koefisien 𝛽 menjelaskan efek 𝑥 terhadap kemungkinan berada pada kategori yang lebih rendah atau sama. Jika 𝛽>0: semakin besar 𝑥, semakin tinggi peluang berada di kategori rendah. Jika 𝛽<0: semakin besar 𝑥, semakin besar peluang berada di kategori tinggi. Odds ratio:

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

7.3 Contoh Data

Misal kita memiliki data fiktif tingkat kepuasan pelanggan (1: Tidak Puas, 2: Cukup, 3: Puas) terhadap kecepatan layanan (1-10) :

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

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

Estimasi Model

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

## Call:
## polr(formula = satisfaction ~ speed, data = df, Hess = TRUE)
## 
## Coefficients:
##       Value Std. Error t value
## speed 0.845    0.09771   8.648
## 
## Intercepts:
##                  Value  Std. Error t value
## Tidak Puas|Cukup 1.2120 0.3990     3.0379 
## Cukup|Puas       4.3249 0.5332     8.1113 
## 
## Residual Deviance: 241.3823 
## AIC: 247.3823

Nilai P value

(ctable <- coef(summary(model_ord)))

##                      Value Std. Error  t value
## speed            0.8449704  0.0977117 8.647587
## Tidak Puas|Cukup 1.2120113  0.3989650 3.037889
## Cukup|Puas       4.3249415  0.5332006 8.111284

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

##                      Value Std. Error  t value p value
## speed            0.8449704  0.0977117 8.647587  0.0000
## Tidak Puas|Cukup 1.2120113  0.3989650 3.037889  0.0024
## Cukup|Puas       4.3249415  0.5332006 8.111284  0.0000

Dari hasil p value ini dapat dilihat bahwa setiap variabel memberikan pengaruh yang signifikan terhadap model

Prediksi Peluang

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

##    Tidak Puas      Cukup      Puas
## 1 0.046849132 0.47815236 0.4749985
## 2 0.020677576 0.30126016 0.6780623
## 3 0.008988488 0.16041632 0.8305952
## 4 0.003881085 0.07667453 0.9194444
## 5 0.001670897 0.03460008 0.9637290

Hasil ini memberikan kejelasan peluang ketika tingkat pelayanannya bernilai 5 sampai 9. Sebagai contoh, ketika kecepatan pelayanannya 5 maka peluang konsumen memberikan penilaian puas adalah 0,4749. Atau contoh selanjutnya ketika kecepatan pelayanannya 9, maka peluang konsumen memberikan penilaian tidak puas adalah 0,0016.

Goodness-of-Fit dan Proportional Odds

Model cumulative logit mengasumsikan efek prediktor sama untuk setiap cutoff. Jika efek predikstor ini tidak sama, maka dapat dipertimbangkan untuk menggunakan model lain seperti Generalized Ordinal Model, Adjacent-Category Logit dan Continuation Ratio (Sequential) logit karena model ini cdapat digunakan ketika asumsi proportional odds tidak terpenuhi.

7.4 Asumsi Paralelisme dalam Regresi Logistik Ordinal

Dalam Model Regresi Logistik Ordinal, metode yang sering diguakan adalah cumulative logit dengan asumsi proporsional odds. Asumsi ini juga disebut sebagai Asumsi Paralelisme. Asumsi ini menyatakan bahwa koeisien setiap katefori kumulatif variabel respon sama.

Bentuk umum model:

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

Jika menggunakan modal cumulative logit ini, maka asumsi paralelisme harus terpenuhi. Jika tidak bisa menyebabkan salah interpretasi . Pengujiannya dapat menggunakan uji likelihood ratio test atau uji brant.

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

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

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

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

Dari pengujian ini dapat dilihat nilai probability adalah 0,44. Pada pengujian brant, probability ini sama dengan nilai p value. Sehingga dapat dilihat bahwa kita gagal menolak H0 yang berarti data memenuhi asumsi paralelisme.

7.5 Review

Regresi Logistik Ordinal adalah metode statistik yang digunakan pada variabel respon ordinal untuk memprediksi. Umumnya, model ini menggunakan model cumulative logit model yang harus memenuhi asumsi proportional odds. Jika model tidak memeuhi asumsi proportional odds, maka dapat dilakukan dengan model lain seperti Adjacent Category Logit, Contionuation Ratio (Sequential) Logit, Generalized Ordinal Model dan lainnya. Dalam pengujian, diperoleh prediksi peluang kejadian sebuah kegiatan. Kemudian juga dilakukan pengujian asumsi proportional odds (paralelisme) menggunakan library brant.

8 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 mana yang independen, karena seluruh variabel diperlakukan secara simetris. Model ini lebih cocok 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 serta interaksi di antara variabel-variabel tersebut. Misalnya, dalam tabel kontingensi tiga arah (misalnya: jenis kelamin, status merokok, dan penyakit paru), model log-linier dapat membedakan apakah interaksi dua variabel cukup menjelaskan data, ataukah diperlukan interaksi tiga arah untuk menangkap struktur asosiasinya. Penyesuaian model dapat dilakukan menggunakan metode likelihood ratio test untuk membandingkan model sederhana dengan model lebih kompleks.

Di sisi lain, regresi logistik adalah pendekatan paling umum ketika terdapat satu variabel kategorik yang secara eksplisit dianggap sebagai variabel dependen (misalnya, kejadian penyakit: ya/tidak), dan satu atau lebih variabel kategorik atau numerik sebagai prediktor. Model ini memodelkan logit dari probabilitas kejadian (yaitu log odds), dan sangat berguna dalam studi observasional dan eksperimental untuk menjelaskan atau memprediksi peluang suatu outcome. Regresi logistik juga memiliki ekstensi untuk outcome kategorik lebih dari dua kelas, seperti regresi logistik multinomial dan regresi logistik ordinal.

Tabel Kontingensi

Tabel Kontingensi menyajikan frekuensi dari kombinasi kategori variabel. Contoh

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

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

Model Log Linear

Model loglinear memodelkan logaritma dari ekspektasi frekuensi sel dalam tabel kontingensi.

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

Model ini cocok untuk menyelidiki asosiasi dan independensi antar variabel. Tidak ada membedakan secara jelas mana yang variabel respon mana penjelas serta umum digunakan pada tabel multi dimensi.

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

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

Model Regresi Logistik

Ini digunakan untuk data respon biner dengan

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

Model ini bertujuan untuk memprediksi peluang suatu outcome dan sering digunakan pada studi observasional atau eksperimental.

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

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

Perbandingan Ketika Pendekatan

Aspek	Tabel Kontingensi	Model Loglinear	Regresi Logistik
Tujuan	Deskripsi frekuensi	Deteksi asosiasi	Prediksi probabilitas
Variabel dependen	Tidak ada	Tidak ada (simetris)	Ada (eksplisit)
Distribusi	Tidak diasumsikan	Poisson (frekuensi sel)	Binomial (probabilitas)
Bentuk Model	Tidak ada	GLM: log(mu) ~ efek	GLM: logit(p) ~ prediktor
Cocok untuk	Eksplorasi awal	Tabel > 2 variabel	Studi prediktif

8.1 Konsep Awal

8.1.1 Tabel Kontingensi dan Model Log Linear

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

# Contoh tabel 2x2
matrix(c(30, 20, 50, 70), nrow=2,
       dimnames = list(Obat = c("Timolol", "Placebo"),
                       Serangan = c("Ya", "Tidak")))

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

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

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

8.1.2 Model Saturated

Model saturated atau model penuh menyertakan seluruh efek utama dan interaksi sehingga cocok sempurna terhadap data dan tidak mengasumsikan independensi antar variabel

Contoh formulasi untuk tabel 2x2:

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

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

Model saturated dapat dipasang dengan loglm dari package {MASS}:

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

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

8.1.3 Model Independent

Model independen mengasumsikan bahwa tidak ada interaksi antara variabel:

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

Model ini menguji hipotesis bahwa variabel X dan Y saling independen.

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

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

8.1.4 Odd Rasio dan Interpretasi

Odd Rasio untuk tabel 2x2 :

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

Interpretasi nilai OR:

• OR = 1: Tidak ada asosiasi

• OR > 1: Asosiasi positif

• OR < 1: Asosiasi negatif

8.1.5 Estimasi Parameter

Dalam model saturated, estimasi dilakukan dengan pembatasan seperti sum-to-zero dan estimasi parameter juga dilakukan dengan iterative proportional fitting (IPF)

# Estimasi odds ratio dan log-odds
logOR <- log((data[1,1] * data[2,2]) / (data[1,2] * data[2,1]))
logOR

## [1] -0.4183685

8.1.6 Model lebih sederhana dan 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:
##  ~Obat + Serangan 
## Model 2:
##  ~Obat * Serangan 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   2.087576  1                                    
## Model 2   0.000000  0   2.087576         1         0.1485
## Saturated 0.000000  0   0.000000         0         1.0000

8.1.7 Studi Kasus Kepercayaan terhadap Surga

Studi dari Agresti (2019) membahas hubungan antara kebahagiaan dan kepercayaan terhadap kehidupan akhirat.

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

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

loglm(~ Kebahagiaan + Surga, data = data_survey)

## Call:
## loglm(formula = ~Kebahagiaan + Surga, data = data_survey)
## 
## Statistics:
##                        X^2 df  P(> X^2)
## Likelihood Ratio 0.8911136  2 0.6404675
## Pearson          0.8836760  2 0.6428538

P value untuk kedua uji jauh lebih besar dari 0,05 yang berarti tidak ada bukti untuk menolak hipotesis null. Artinya tidak ada hubungan secara statstik antara kebahagiaan dengan kepercayaan terhadap surga.

8.2 Model Log Linear 2 arah

8.2.1 Model log linear pada tabel kontingensi

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

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

8.2.2 Perbedaan utama antara model log linear dengan 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).

8.2.3 Estimasi Parameter Log Linear

Sistem Persamaan Log Linear dan constraint sum to zero

\[ \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} \\ \\ \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 \end{aligned} \]

Rumus estimasi parameter dengan sum to zero constraint

\[ \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} \left[ \log \mu_{12} - \log \mu_{11} - \log \mu_{22} + \log \mu_{21} \right] \]

8.2.4 Analisis Data Tabel Kontingensi (Contoh Kasus)

Diberikan data:

Merokok	Sakit	Sehat
Ya	30	20
Tidak	10	40

8.2.4.1 Bentuk model Log Linear

model log linearnya adalah

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

dengan constraint sum to zero :

\[ \sum_i \lambda^A_i = 0 \]

\[ \sum_j \lambda^B_j = 0 \]

\[ \sum_{i,j} \lambda^{AB}_{ij} = 0 \]

8.2.4.2 Estimasi Parameter Model (Manual, Sum to Zero)

Misal

• A1 = Merokok (Ya), A2 = Tidak

• B1 = Sakit, B2 = Sehat

Observasi:

• n11=30, n12=20n

• n21=10, n22=40n

\[ \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} \\ \\ \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 \end{aligned} \]

Langkah-langkah:

1. Hitung rata-rata log frekuensi sel:

\[ \begin{aligned} \lambda &= \frac{1}{4} \sum_{i=1}^2 \sum_{j=1}^2 \log(n_{ij}) \\ &= \frac{1}{4} \big[ \log(30) + \log(20) + \log(10) + \log(40) \big] \\ &= 3.0971 \end{aligned} \]

2. Efek utama A (Merokok):

\[ \begin{aligned} \lambda^A_1 &= \frac{1}{2} \big[ (\log 30 + \log 20) - (\log 10 + \log 40) \big] \\ &= \frac{1}{2} \big[ (3.4012 + 2.9957) - (2.3026 + 3.6889) \big] \\ &= \frac{1}{2} (6.3969 - 5.9915) \\ &= \frac{1}{2} (0.4054) \\ &= 0.2027 \\ \lambda^A_2 &= -0.2027 \end{aligned} \]

3. Efek utama B (Status):

\[ \begin{aligned} \lambda^B_1 &= \frac{1}{2} \big[ (\log 30 + \log 10) - (\log 20 + \log 40) \big] \\ &= \frac{1}{2} \big[ (3.4012 + 2.3026) - (2.9957 + 3.6889) \big] \\ &= \frac{1}{2} (5.7038 - 6.6846) \\ &= \frac{1}{2} (-0.9808) \\ &= -0.4904 \\ \lambda^B_2 &= +0.4904 \end{aligned} \]

4. Efek interaksi:

\[ \begin{aligned} \lambda^{AB}_{11} &= \frac{1}{4} \big[ \log(30) - \log(20) - \log(10) + \log(40) \big] \\ &= \frac{1}{4} \big[ 3.4012 - 2.9957 - 2.3026 + 3.6889 \big] \\ &= \frac{1}{4} (1.7918) \\ &= 0.4479 \\ \lambda^{AB}_{12} &= -\lambda^{AB}_{11} = -0.4479 \\ \lambda^{AB}_{21} &= -0.4479 \\ \lambda^{AB}_{22} &= +0.4479 \end{aligned} \]

8.2.4.3 Hitung Odds Ratio dan Interval Kepercayaan

\[ \begin{aligned} \text{OR} &= \frac{n_{11} n_{22}}{n_{12} n_{21}} = \frac{30 \times 40}{20 \times 10} = \frac{1200}{200} = 6 \end{aligned} \]

Log odds ratio:

\[ \log(\text{OR}) = \log(6) = 1.7918 \] Standard error (SE):

\[ \begin{aligned} SE &= \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}} \\ &= \sqrt{\frac{1}{30} + \frac{1}{20} + \frac{1}{10} + \frac{1}{40}} \\ &= \sqrt{0.0333 + 0.05 + 0.1 + 0.025} \\ &= \sqrt{0.2083} \\ &= 0.4564 \end{aligned} \]

95% Confidence Interval for log(OR):

\[ \begin{aligned} \log(\text{OR}) \pm 1.96 \times SE &= 1.7918 \pm 1.96 \times 0.4564 \\ &= (1.7918 - 0.895), (1.7918 + 0.895) \\ &= (0.8968, 2.6868) \end{aligned} \]

Back-transform to get CI for OR:

\[ \begin{aligned} \text{Lower} &= \exp(0.8968) = 2.452 \\ \text{Upper} &= \exp(2.6868) = 14.68 \end{aligned} \]

Jadi, OR = 6 (95% CI: 2.45 – 14.68)

8.2.5 Fitting Model dengan log Linear dengan R

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

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

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

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

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

## 
## Call:
## glm(formula = Freq ~ Merokok + Status, family = poisson, data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.996e+00  1.871e-01  16.013   <2e-16 ***
## MerokokTidak 3.892e-10  2.000e-01   0.000    1.000    
## StatusSehat  4.055e-01  2.041e-01   1.986    0.047 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 21.288  on 3  degrees of freedom
## Residual deviance: 17.261  on 1  degrees of freedom
## AIC: 43.036
## 
## Number of Fisher Scoring iterations: 4

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

## 
## Call:
## glm(formula = Freq ~ Merokok * Status, family = poisson, data = data)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.4012     0.1826  18.629  < 2e-16 ***
## MerokokTidak              -1.0986     0.3651  -3.009  0.00262 ** 
## StatusSehat               -0.4055     0.2887  -1.405  0.16015    
## MerokokTidak:StatusSehat   1.7918     0.4564   3.926 8.65e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.1288e+01  on 3  degrees of freedom
## Residual deviance: 3.9968e-15  on 0  degrees of freedom
## AIC: 27.775
## 
## Number of Fisher Scoring iterations: 3

8.2.6 Interpretasi parameter

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

• Efek “Merokok” dan “Status” menunjukkan perbedaan log frekuensi antar kategori.

• Interaksi signifikan menunjukkan adanya asosiasi antara Merokok dan Status.

Nilai log(6) =1.79 itu sama dengan efek interaksi ouput R

8.2.7 Analisis Data Tabel Kontingensi 2x3 (Contoh Kasus)

Suatu survei dilakukan untuk mengetahui hubungan antara Jenis Kelamin (Laki-laki/Perempuan) dan Kategori BMI (Kurus/Normal/Gemuk):

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

8.2.8 Bentuk Model Log-Linear untuk Tabel 2x3

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

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

8.2.9 Fitting model log linear di R

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

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

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

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

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

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

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

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

8.2.10 Interpretasi

Model tanpa interaksi: - Jika deviance tidak signifikan, maka Jenis Kelamin dan BMI independen.

• Intercept: log frekuensi pada kategori referensi (Laki-laki, Kurus)

• Koefisien JenisKelaminPerempuan: perbedaan log-frekuensi antara Perempuan vs Laki-laki (pada Kurus)

• Koefisien BMI: perbedaan log-frekuensi kategori BMI (Normal/Gemuk) terhadap Kurus (pada Laki-laki)

Model dengan interaksi:

• Jika koefisien interaksi signifikan, berarti ada hubungan/asosiasi antara Jenis Kelamin dan BMI. Artinya distribusi BMI berbeda antara Laki-laki dan Perempuan.

Contoh interpretasi hasil (misal):

• Jika koefisien JenisKelaminPerempuan negatif: proporsi Perempuan pada kategori referensi lebih kecil dibanding Laki-laki.

• Jika koefisien BMI_Normal positif: kemungkinan seseorang Normal lebih tinggi daripada Kurus (pada Laki-laki).

• Jika model interaksi signifikan, pola distribusi BMI berbeda antara Laki-laki dan Perempuan.

8.3 Model Log Linear Tiga Arah

Pada pembahasan sebelumnya, kita telah memahami bahwa salah satu tujuan utama dari penyusunan model log linear adalah untuk mengestimasi parameter-parameter yang menjelaskan hubungan di antara variabel-variabel kategorik.

Pada materi kali ini, kita akan membahas model log linear yang lebih kompleks, yaitu model log linear untuk tabel kontingensi tiga arah. Model ini melibatkan tiga variabel kategorik, sehingga kemungkinan interaksi yang dapat terjadi di dalam model pun menjadi lebih banyak. Dalam konteks ini, interaksi paling tinggi yang dapat dimodelkan adalah interaksi tiga arah, yaitu interaksi yang melibatkan ketiga variabel secara bersamaan.

8.3.1 Model Log-Linear untuk Tabel Tiga Arah

Model log-linear yang melibatkan tiga variabel kategorik (misal: X, Y, dan Z) dapat dibangun dalam berbagai bentuk model, tergantung pada tingkat interaksi yang ingin dimasukkan. Berikut adalah beberapa alternatif model log-linear yang umum digunakan:

1. 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 (X, Y, dan Z).

2. 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 mengakomodasi interaksi dua arah antar variabel tanpa
memasukkan interaksi tiga arah.

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

4. Model Joint Independence

   Independensi antara X & Y:

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

*Independensi antara X & Z*:

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

*Independensi antara Y & Z*:

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

Model Tanpa Interaksi

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

Model ini hanya memasukkan efek utama tanpa interaksi antar variabel.

8.3.2 Pengujian Interaksi dalam Model Log-Linear Tiga Arah

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:

1. Pengujian Interaksi Tiga Arah (XYZ):

   • Bandingkan model saturated dengan model homogenous.

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

8.3.3 Contoh Soal

8.3.3.1 Soal

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

8.3.3.2 Tabel Data Survei

Fundamentalisme	Jenis Kelamin	Mendukung	Menolak	Total
Fundamentalist	Laki-laki	128	32	160
Fundamentalist	Perempuan	123	73	196
Fundamentalist	Total	251	105	356
Moderate	Laki-laki	182	56	238
Moderate	Perempuan	168	105	273
Moderate	Total	350	161	511
Liberal	Laki-laki	119	49	168
Liberal	Perempuan	111	70	181
Liberal	Total	230	119	349

Keterangan:

• Fundamentalisme: Fundamentalist, Moderate, Liberal

• Jenis Kelamin: Laki-laki, Perempuan

• Sikap: Mendukung (Favor), Menolak (Oppose) hukuman mati

8.3.3.3 Analisis Log-Linear untuk Tabel Tiga Arah

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

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

data <- data.frame(
  Fundamentalisme = z.fund,
  Jenis_Kelamin   = x.sex,
  Sikap           = y.fav,
  Frekuensi       = counts

)
data

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

# Membuat tabel 3d
table3d <- xtabs(Frekuensi ~ Fundamentalisme + Jenis_Kelamin + Sikap, data = data)
ftable(table3d)

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

Dari sini kita dapat melihat tabel 3 arah dari studi kasus tersebut. Analisis Log-Linear dilakukan model log-linear dan membandingkan kecocokan model (parsimonious model)

8.3.3.4 Uji Model Interaksi Tiga Arah (Saturated vs Homogenous)

8.3.3.4.1 Penentuan Kategori Referensi

##=============================##
# Penentuan kategori reference
##=============================##
x.sex  <- relevel(x.sex, ref = "2F")
y.fav  <- relevel(y.fav, ref = "2opp")
z.fund <- relevel(z.fund, ref = "3lib")

Bagian ini adalah penetuan baseline dari variabel yang ada yakni pada variabel gender, perempuan digunakan sebagai baseline. Pada variabel sikap hukuman mati, menolak dianggap baseline dan liberal menjadi baseline untuk variabel tingkat fundamentalisme.

Model Saturated Model log-linear saturated memasukkan semua interaksi hingga tiga arah:

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

Dengan ini dibentuk :

# Model saturated
model_saturated <- glm(counts ~ x.sex + y.fav + z.fund +
             x.sex*y.fav + x.sex*z.fund + y.fav*z.fund +
             x.sex*y.fav*z.fund,
             family = poisson(link = "log"))
summary(model_saturated)

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

exp(model_saturated$coefficients)

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

8.3.3.4.2 Ringkasan Model

Model yang digunakan adalah model log-linear saturated dengan semua efek utama, interaksi dua arah, dan interaksi tiga arah. Model ini memodelkan hubungan antara jenis kelamin (x.sex), sikap terhadap hukuman mati (y.fav), dan tingkat fundamentalisme (z.fund) terhadap frekuensi responden.

8.3.3.4.3 Hasil Estimasi Koefisien

Pada output diatas dapat dlihat hasil estimasi perhitungannya. diperoleh bahwa estimasi untuk intercept sebesar 4,24 yang berarti rata rata log jumlah kasus untuk kategori referensi adalah 4,24. Kemudian laki laki memiliki expected count sebesar 0,7 kali perempuan dalam kategoru referensi lainnya. dan dapat dilihat hasil variabel lain.

8.3.3.4.4 Interpretasi signifikansi

Dapat dilihat pada output bahwa intercept, pernyataan mendukung, moderate memberikan pengaruh yang signifikan pada taraf 5%. Sedangkan yang lainnya tidak memberikan pengaruh signifikan pada taraf tersebut.

8.3.3.4.5 Goodness-of-Fit

Residual deviance ≈ 0 menandakan model saturated benar-benar fit terhadap data (seluruh variasi data dijelaskan oleh model). AIC = 100.14 dapat digunakan untuk perbandingan dengan model yang lebih sederhana.

8.3.3.4.6 Kesimpulan

• Model saturated ini sangat fit dengan data, namun tidak semua parameter/interaksi signifikan.

• Efek utama yang paling signifikan adalah:

  • Sikap mendukung hukuman mati (expected count 1.6x lebih tinggi dari yang menolak)

  • Kelompok Moderate (expected count 1.5x lebih tinggi dari Liberal)

• Tidak ditemukan bukti kuat interaksi dua atau tiga arah yang signifikan.

• Model yang lebih sederhana (tanpa interaksi tiga arah) perlu dipertimbangkan untuk model final yang lebih parsimonious.

Catatan interpretasi:

• Nilai exp(coef) menyatakan rasio ekspektasi (expected count ratio) dibandingkan baseline.

• Efek positif → menaikkan expected count; Efek negatif → menurunkan expected count.

• Koefisien signifikan pada p-value < 0.05.

8.3.3.5 Model Homogenous

Model log-linear homogenous memasukkan semua efek utama dan semua interaksi dua arah, tanpa interaksi tiga arah. Secara matematis, model ini dapat dituliskan sebagai berikut:

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

# Homogenous Model
model_homogenous <- glm(counts ~ x.sex + y.fav + z.fund +
              x.sex*y.fav + x.sex*z.fund + y.fav*z.fund,
              family = poisson(link = "log"))
summary(model_homogenous)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     x.sex * z.fund + y.fav * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            4.31096    0.10522  40.972  < 2e-16 ***
## x.sex1M               -0.51575    0.13814  -3.733 0.000189 ***
## y.fav1fav              0.35707    0.12658   2.821 0.004788 ** 
## z.fund1fund           -0.06762    0.14452  -0.468 0.639854    
## z.fund2mod             0.33196    0.13142   2.526 0.011540 *  
## x.sex1M:y.fav1fav      0.66406    0.12728   5.217 1.81e-07 ***
## x.sex1M:z.fund1fund   -0.16201    0.15300  -1.059 0.289649    
## x.sex1M:z.fund2mod    -0.08146    0.14079  -0.579 0.562887    
## y.fav1fav:z.fund1fund  0.23873    0.16402   1.455 0.145551    
## y.fav1fav:z.fund2mod   0.13081    0.14951   0.875 0.381614    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.361  on 11  degrees of freedom
## Residual deviance:   1.798  on  2  degrees of freedom
## AIC: 97.934
## 
## Number of Fisher Scoring iterations: 3

8.3.3.5.1 Uji Hipotesis: Apakah Ada Interaksi Tiga Arah? (Saturated vs Homogenous)

Pengujian ini menggunakan residual deviance dari kedua model (saturated dan homogenous).Jika selisih deviance model homogenous dengan model saturated lebih besar daripada daerah kritis penolaknnya, maka dianggap tidak ada pengaruh dari interaksi 3 arah.

Hipotesis

• H0: Tidak ada interaksi tiga arah (model homogenous sudah cukup)

• H1: Ada interaksi tiga arah (model saturated diperlukan)

Hitung Selisih Deviance

# Deviance antar model
Deviance.model <- model_homogenous$deviance - model_saturated$deviance
Deviance.model

## [1] 1.797977

# Hitung derajat bebas
# Derajat bebas = db model homogenous - db model saturated
derajat.bebas <- (model_homogenous$df.residual - model_saturated$df.residual)
derajat.bebas

## [1] 2

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

Interpretasi Pada taraf nyata 5%, belum cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi tiga arah antara jenis kelamin, fundamentalisme, dan pendapat mengenai hukuman mati.

Catatan:

• Model pengurang adalah model yang lebih lengkap (lebih banyak parameter, df lebih kecil), yaitu model saturated.

• Derajat bebas dihitung dari selisih derajat bebas model homogenous dan saturated.

• Keputusan berdasarkan perbandingan deviance model dengan chi-square tabel.

Dari perhitungan tersebut, maka dapat diuji hipotesisnya yakni :

Hipotesis

H0 : (Tidak ada interaksi tiga arah; model yang terbentuk adalah model homogenous)

H1 : (Ada interaksi tiga arah; model yang terbentuk adalah model saturated)

Tingkat Signifikansi

α = 5%

Statistik Uji

Δ Deviance = Deviance model homogenous − Deviance model saturated

= 1.798 − 0.00

= 1.798

db = db model homogenous − db model saturated

= 2−0

= 2

Daerah Penolakan

Tolak H0 jika Δ Deviance > χ2(0.05,2)

= χ2(0.05,2)

= 5.991

Keputusan

Karena 1.798 < 5.991, maka terima H0

Interpretasi

Pada taraf nyata 5%, belum cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi tiga arah antara jenis kelamin, fundamentalisme, dan pendapat mengenai hukuman mati.

Ingat, dalam membuat selisih deviance, model yang menjadi pengurang adalah model yang lebih lengkap (parameter yang lebih banyak atau derajat bebasnya lebih kecil).

Makin banyak parameter, makin kecil derajat bebasnya, karena:

• db = d banyaknya amatan (atau perkalian dimensi tabel kontingensi, misal 2×2×3=12) dikurangi banyaknya parameter (koefisien, termasuk intercept).

Cek di output R ada berapa banyak coefficients-nya (termasuk intercept) untuk menghitung derajat bebas yang benar.

8.3.3.6 Uji Model Interaksi Dua Arah (Homogenous vs Conditional On X)

8.3.3.6.1 Model Conditional on X

Model log-linear conditional pada X memasukkan efek utama dan interaksi dua arah antara X dengan Y dan X dengan Z, tanpa interaksi antara Y dengan Z maupun interaksi tiga arah.

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

# Conditional Association on X
model_conditional_X <- glm(counts ~ x.sex + y.fav + z.fund +
              x.sex*y.fav + x.sex*z.fund,
              family = poisson(link = "log"))
summary(model_conditional_X)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     x.sex * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          4.23495    0.08955  47.293  < 2e-16 ***
## x.sex1M             -0.52960    0.13966  -3.792 0.000149 ***
## y.fav1fav            0.48302    0.08075   5.982 2.20e-09 ***
## z.fund1fund          0.07962    0.10309   0.772 0.439916    
## z.fund2mod           0.41097    0.09585   4.288 1.81e-05 ***
## x.sex1M:y.fav1fav    0.65845    0.12708   5.181 2.20e-07 ***
## x.sex1M:z.fund1fund -0.12841    0.15109  -0.850 0.395405    
## x.sex1M:z.fund2mod  -0.06267    0.13908  -0.451 0.652274    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   3.9303  on  4  degrees of freedom
## AIC: 96.067
## 
## Number of Fisher Scoring iterations: 4

Pengujian Ada Tidaknya Interaksi Antara Y dan Z (Homogenous Model vs Conditional Association on X)

Hipotesis

H0 : (Tidak ada interaksi antara pendapat hukuman mati (Y) dan fundamentalisme (Z))

H1 : (Ada interaksi antara pendapat hukuman mati (Y) dan fundamentalisme (Z))

Tingkat Signifikansi

α = 5%

Statistik Uji

Δ Deviance = Deviance model conditional on X − Deviance model homogenous

= 3.903 − 1.798

= 2.132

db = db model conditional on X − db model homogenous

= 4 − 2

= 2

Daerah Penolakan

Tolak H0 jika Δ Deviance > χ2(0.05,2)

= 5.991

Keputusan

Karena 2.132 < 5.991, maka terima H0

Kesimpulan

Dengan taraf nyata 5%, belum cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi antara pendapat tentang hukuman mati dan fundamentalisme. Dengan kata lain, model yang terbentuk adalah model tanpa parameter λ

# Pengujian hipotesis

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

## [1] 2.132302

8.3.3.6.2 Pengujian Selisih Deviance (Conditional on X vs Homogenous)

Langkah-langkah uji hipotesis menggunakan residual deviance:

# Selisih deviance antar model
Deviance.model <- model_conditional_X$deviance - model_homogenous$deviance
Deviance.model

## [1] 2.132302

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

## [1] 2

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

## [1] 5.991465

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

## [1] "Terima"

Interpretasi Karena nilai Deviance.model = 2.13 lebih kecil dari nilai kritis chi-square tabel = 5.99 (dengan df = 2, alpha = 0.05), maka keputusan uji adalah “Terima”.

Pada taraf nyata 5%, belum cukup bukti untuk menolak H0, atau dengan kata lain tidak ada interaksi antara pendapat mengenai hukuman mati (Y) dan fundamentalisme (Z). Model yang terbentuk cukup hanya sampai dua interaksi dengan X (conditional on X), sehingga interaksi Y*Z tidak signifikan secara statistik.

8.3.3.7 Uji Model Interaksi Dua Arah (Homogenous vs Conditional On Y)

8.3.3.7.1 Model Conditional on Y

Model log-linear conditional pada Y memasukkan efek utama dan interaksi dua arah antara X dengan Y dan Y dengan Z, tanpa interaksi antara X dengan Z maupun interaksi tiga arah.

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

# Conditional Association on Y
model_conditional_Y <- glm(counts ~ x.sex + y.fav + z.fund +
              x.sex*y.fav + y.fav*z.fund,
              family = poisson(link = "log"))
summary(model_conditional_Y)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     y.fav * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            4.33931    0.09919  43.748  < 2e-16 ***
## x.sex1M               -0.59345    0.10645  -5.575 2.48e-08 ***
## y.fav1fav              0.37259    0.12438   2.996  0.00274 ** 
## z.fund1fund           -0.12516    0.13389  -0.935  0.34989    
## z.fund2mod             0.30228    0.12089   2.500  0.01240 *  
## x.sex1M:y.fav1fav      0.65845    0.12708   5.181 2.20e-07 ***
## y.fav1fav:z.fund1fund  0.21254    0.16205   1.312  0.18966    
## y.fav1fav:z.fund2mod   0.11757    0.14771   0.796  0.42606    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   2.9203  on  4  degrees of freedom
## AIC: 95.057
## 
## Number of Fisher Scoring iterations: 4

Pengujian Ada Tidaknya Interaksi antara X dan Z (Homogenous model vs Conditional Association on Y)

• Hipotesis

  • H0 : (Tidak ada interaksi antara jenis kelamin (X) dan fundamentalisme (Z))

  • H1 : (Ada interaksi antara jenis kelamin (X) dan fundamentalisme (Z))

• Tingkat Signifikansi

  • α = 5%

• Statistik Uji

  • Δ Deviance = Deviance model conditional on Y − Deviance model homogenous

    = 2.9203 − 1.798

    = 1.1223

  • db = db model conditional on Y − db model homogenous

    = 4 − 2

    = 2

• Daerah Penolakan

  • Tolak H0 jika Δ Deviance > χ2(0.05,2)

    = 5.991

• Keputusan

  • Karena 1.1223 < 5.991, maka terima H0

• Kesimpulan

  • Dengan taraf nyata 5%, belum cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi antara jenis kelamin dan fundamentalisme. Dengan kata lain, model yang terbentuk adalah model tanpa parameter λXZ

8.3.3.7.2 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] 1.122315

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

## [1] 2

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

## [1] 5.991465

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

## [1] "Terima"

Interpretasi Karena nilai Deviance.model = 1.12 lebih kecil dari nilai kritis chi-square tabel = 5.99 (df = 2, alpha = 0.05), maka keputusan uji adalah “Terima”.

Kesimpulan: Pada taraf nyata 5%, belum cukup bukti untuk menolak H0. Artinya, tidak ada interaksi antara jenis kelamin (X) dan fundamentalisme (Z) yang signifikan secara statistik. Model tanpa parameter λ(XZ) sudah cukup baik untuk data ini.

8.3.3.8 Uji Model Interaksi Dua Arah (Homogenous VS Conditional On Z)

8.3.3.8.1 Model Conditional on Z

Model log-linear conditional pada Z memasukkan efek utama dan interaksi dua arah antara X dengan Z dan Y dengan Z, tanpa interaksi antara X dengan Y maupun interaksi tiga arah.

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

# Conditional Association on Z
model_conditional_Z <- glm(counts ~ x.sex + y.fav + z.fund +
              x.sex*z.fund + y.fav*z.fund,
              family = poisson(link = "log"))
summary(model_conditional_Z)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * z.fund + 
##     y.fav * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            4.12255    0.10518  39.195  < 2e-16 ***
## x.sex1M               -0.07453    0.10713  -0.696    0.487    
## y.fav1fav              0.65896    0.11292   5.836 5.36e-09 ***
## z.fund1fund           -0.06540    0.15126  -0.432    0.665    
## z.fund2mod             0.33196    0.13777   2.410    0.016 *  
## x.sex1M:z.fund1fund   -0.12841    0.15109  -0.850    0.395    
## x.sex1M:z.fund2mod    -0.06267    0.13908  -0.451    0.652    
## y.fav1fav:z.fund1fund  0.21254    0.16205   1.312    0.190    
## y.fav1fav:z.fund2mod   0.11757    0.14771   0.796    0.426    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.361  on 11  degrees of freedom
## Residual deviance:  29.729  on  3  degrees of freedom
## AIC: 123.87
## 
## Number of Fisher Scoring iterations: 4

Pengujian Ada Tidaknya Interaksi antara X dan Y (Homogenous model vs Conditional Association on Z)

• Hipotesis

  • H0 : (Tidak ada interaksi antara jenis kelamin (X) dan pendapat tentang hukuman mati (Y))

  • H1 : (Ada interaksi antara jenis kelamin (X) dan pendapat tentang hukuman mati (Y))

• Tingkat Signifikansi

  • α = 5%

• Statistik Uji

  • Δ Deviance = Deviance model conditional on Z − Deviance model homogenous

    = 29.729 − 1.798

    = 27.931

  • db = db model conditional on Z − db model homogenous

    = 3 − 2

    = 1

• Daerah Penolakan

  • Tolak H0 jika Δ Deviance > χ2(0.05,1)

    = 3.841

• Keputusan

  • Karena 27.931 > 3.841, maka tolak H0

• Kesimpulan

  • Dengan taraf nyata 5 persen, ada interaksi antara jenis kelamin dan pendapat tentang hukuman mati. Dengan kata lain, model yang terbentuk adalah model dengan parameter λXY.

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

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

## [1] 27.93095

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

## [1] 1

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

## [1] 3.841459

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

## [1] "Tolak"

Interpretasi Karena nilai Deviance.model = 27.93 jauh lebih besar dari nilai kritis chi-square tabel = 3.84 (df = 1, alpha = 0.05), maka keputusan uji adalah “Tolak”.

Kesimpulan: Pada taraf nyata 5%, terdapat bukti yang cukup untuk menolak H0H0. Artinya, ada interaksi yang signifikan antara jenis kelamin (X) dan pendapat tentang hukuman mati (Y). Dengan kata lain, model terbaik yang terbentuk adalah model yang menyertakan parameter interaksi λ(XY).

8.3.3.9 Pemilihan Model Terbaik

8.3.3.9.1 Ringkasan Model Log Linier

Model	Parameter	Deviance	Jumlah Parameter	df	AIC
Saturated	Full Faktor	0.00	12	0	100.14
Homogenous	Tanpa Interaksi XYZ	1.798	10	2	97.934
Conditional on X	Faktor Utama dan Interaksi X	3.9303	8	4	96.067
Conditional on Y	Faktor Utama dan Interaksi Y	2.9203	8	4	95.057
Conditional on Z	Faktor Utama dan Interaksi Z	29.729	9	3	123.87

8.3.3.9.2 Ringkasan Pengujian Interaksi 3 Arah dan 2 Arah

Interaksi	Pengujian	ΔΔ deviance	ΔΔ df	Chi-square Tabel	Keputusan	Keterangan
XYZ	Saturated vs Homogenous	1.798	2	5.991	Tidak Tolak H0	tidak ada interaksi
YZ	Conditional on X vs Homogenous	2.1323	2	5.991	Tidak Tolak H0	tidak ada interaksi
XZ	Conditional on Y vs Homogenous	1.1223	2	5.991	Tidak Tolak H0	tidak ada interaksi
XY	Conditional on Z vs Homogenous	27.931	1	3.841	Tolak H0	ada interaksi

8.3.3.9.3 Kesimpulan Pemilihan Model Terbaik

Dari hasil di atas diketahui bahwa asosiasi yang nyata hanya terdapat antara jenis kelamin dan pendapat mengenai hukuman mati. Sehingga, model terbaik adalah:

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

Model terbaik adalah model log-linear tanpa interaksi tiga arah dan hanya memuat interaksi dua arah antara jenis kelamin dan sikap terhadap hukuman mati.

8.3.4 Model Terbaik

Model terbaik dipilih berdasarkan pengujian interaksi yang signifikan, yaitu hanya interaksi dua arah antara jenis kelamin (X) dan sikap terhadap hukuman mati (Y):

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

# Model Terbaik
bestmodel <- glm(counts ~ x.sex + y.fav + z.fund +
                 x.sex*y.fav,
                 family = poisson(link = "log"))
summary(bestmodel)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav, 
##     family = poisson(link = "log"))
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        4.26518    0.07794  54.721  < 2e-16 ***
## x.sex1M           -0.59345    0.10645  -5.575 2.48e-08 ***
## y.fav1fav          0.48302    0.08075   5.982 2.20e-09 ***
## z.fund1fund        0.01986    0.07533   0.264    0.792    
## z.fund2mod         0.38130    0.06944   5.491 4.00e-08 ***
## x.sex1M:y.fav1fav  0.65845    0.12708   5.181 2.20e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   4.6532  on  6  degrees of freedom
## AIC: 92.79
## 
## Number of Fisher Scoring iterations: 4

Dari summary model diatas terlihat bahwa best model memiliki AIC yang lebih rendah dibandingkan saturated, homogeneous, dan conditional model

8.3.4.1 Interpretasi Koefisien Model Terbaik

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

##                          koef   exp_koef
## (Intercept)        4.26517861 71.1776316
## x.sex1M           -0.59344782  0.5524194
## y.fav1fav          0.48302334  1.6209677
## z.fund1fund        0.01985881  1.0200573
## z.fund2mod         0.38129767  1.4641834
## x.sex1M:y.fav1fav  0.65845265  1.9318008

8.3.4.2 Interpretasi Koefisien Model Terbaik

• exp(λX(Pria))=e xp(−0,593) =0,552 → nilai odds
  Tanpa memperhatikan fundamentalisme dan pendapat mengenai hukuman mati, peluang seseorang berjenis kelamin laki-laki adalah 0,55 kali dibandingkan perempuan.
  Atau, peluang seseorang berjenis kelamin perempuan adalah 1 / (0,55) = 1,81 kali dibandingkan laki-laki.

• exp(λY(dukung)) = exp(0,483) = 1,621 → nilai odds
  Tanpa memperhatikan jenis kelamin dan fundamentalisme, peluang seseorang mendukung hukuman mati adalah 1,621 kali dibandingkan yang menolak.

• exp(λZ(funda)) = exp(0,01986) = 1,02 → nilai odds
  Tanpa memperhatikan jenis kelamin dan pendapat mengenai hukuman mati, peluang seseorang fundamentalist adalah 1,02 kali dibandingkan liberal.

• exp(λZ(moderate)) = exp(0,381) = 1,464 → nilai odds
  Tanpa memperhatikan jenis kelamin dan pendapat mengenai hukuman mati, peluang seseorang moderate adalah 1,464 kali dibandingkan liberal.

• exp(λXY(laki,dukung)) = exp(0,658) = 1,932 → nilai odds ratio
  Tanpa memperhatikan fundamentalisme, odds mendukung hukuman mati (dibandingkan menolak) jika dia laki-laki adalah 1,932 kali dibandingkan odds yang sama jika dia perempuan.

8.3.5 Nilai dugaan model terbaik

# Fitted values dari model terbaik
data.frame(
  Fund = z.fund,
  sex = x.sex,
  favor = y.fav,
  counts = counts,
  fitted = bestmodel$fitted.values
)

##     Fund sex favor counts    fitted
## 1  1fund  1M  1fav    128 125.59539
## 2  1fund  1M  2opp     32  40.10855
## 3  1fund  2F  1fav    123 117.69079
## 4  1fund  2F  2opp     73  72.60526
## 5   2mod  1M  1fav    182 180.27878
## 6   2mod  1M  2opp     56  57.57155
## 7   2mod  2F  1fav    168 168.93257
## 8   2mod  2F  2opp    105 104.21711
## 9   3lib  1M  1fav    119 123.12582
## 10  3lib  1M  2opp     49  39.31990
## 11  3lib  2F  1fav    111 115.37664
## 12  3lib  2F  2opp     70  71.17763

8.3.5.1 Perhitungan Manual Nilai Dugaan (Fitted Value) Model Terbaik

Secara manual, nilai fitted value diperoleh dengan cara sebagai berikut:

\[ \begin{aligned} \hat{\mu}_{111} &= \exp(\lambda + \lambda_{x_{1m}} + \lambda_{y_{1fav}} + \lambda_{z_{fund}} + \lambda_{xy_{1m,1fav}}) \\ &= \exp(4.265 - 0.593 + 0.483 + 0.01986 + 0.658) \\ =\exp(4.833) \\ &=125.595 \end{aligned} \]

\[ \begin{aligned} \hat{\mu}_{112} &= \exp(\lambda + \lambda_{x_{1m}} + \lambda_{y_{1fav}} + \lambda_{z_{2mod}} + \lambda_{xy_{1m,1fav}}) \\ &= \exp(4.265 - 0.593 + 0.483 + 0.381 + 0.658) \\ &= \exp(5.195) \\ &= 180.279 \end{aligned} \]

\[ \begin{aligned} \hat{\mu}_{113} &= \exp(\lambda + \lambda_{x_{1m}} + \lambda_{y_{1fav}} + \lambda_{z_{lib}} + \lambda_{xy_{1m,1fav}}) \\ &= \exp(4.265 - 0.593 + 0.483 + 0 + 0.658) \\ &= \exp(4.813) \\ &= 123.126 \end{aligned} \]

\[ \begin{aligned} \hat{\mu}_{121} &= \exp(\lambda + \lambda_{x_{1m}} + \lambda_{y_{2opp}} + \lambda_{z_{fund}} + \lambda_{xy_{1m,2opp}}) \\ & = \exp(4.265 - 0.593 + 0 + 0.01986 + 0) \\ & = \exp(3.692) \\ & = 40.109 \end{aligned} \]

\[ \begin{aligned} \hat{\mu}_{122} &= \exp(\lambda + \lambda_{x_{1m}} + \lambda_{y_{2opp}} + \lambda_{z_{2mod}} + \lambda_{xy_{1m,2opp}}) \\ & = \exp(4.265 - 0.593 + 0 + 0.381 + 0) \\ & = \exp(4.053) \\ & = 57.572 \end{aligned} \]

\[ \begin{aligned} \hat{\mu}_{123} &= \exp(\lambda + \lambda_{x_{1m}} + \lambda_{y_{2opp}} + \lambda_{z_{lib}} + \lambda_{xy_{1m,2opp}}) \\ &= \exp(4.265 - 0.593 + 0 + 0 + 0) \\ &= \exp(3.672) \\ &= 39.320 \end{aligned} \]

\[ \begin{aligned} \hat{\mu}_{211} &= \exp(\lambda + \lambda_{x_{2f}} + \lambda_{y_{1fav}} + \lambda_{z_{fund}} + \lambda_{xy_{2f,1fav}}) \\ &= \exp(4.265 + 0 + 0.483 + 0.01986 + 0) \\ &= \exp(4.768) \\ &= 117.691 \end{aligned} \]

\[ \begin{aligned} \hat{\mu}_{212} &= \exp(\lambda + \lambda_{x_{2f}} + \lambda_{y_{1fav}} + \lambda_{z_{2mod}} + \lambda_{xy_{2f,1fav}}) \\ &= \exp(4.265 + 0 + 0.483 + 0.381 + 0) \\ &= \exp(5.1295) \\ &= 168.933 \end{aligned} \]

\[ \begin{aligned} \hat{\mu}_{213} &= \exp(\lambda + \lambda_{x_{2f}} + \lambda_{y_{1fav}} + \lambda_{z_{lib}} + \lambda_{xy_{2f,1fav}}) \\ &= \exp(4.265 + 0 + 0.483 + 0 + 0) \\ &= \exp(4.748) \\ &= 115.377 \end{aligned} \]

\[ \begin{aligned} \hat{\mu}_{221} &= \exp(\lambda + \lambda_{x_{2f}} + \lambda_{y_{2opp}} + \lambda_{z_{fund}} + \lambda_{xy_{2f,2opp}}) \\ &= \exp(4.265 + 0 + 0 + 0.01986 + 0) \\ &= \exp(4.285) \\ &= 72.605 \end{aligned} \]

\[ \begin{aligned} \hat{\mu}_{222} &= \exp(\lambda + \lambda_{x_{2f}} + \lambda_{y_{2opp}} + \lambda_{z_{2mod}} + \lambda_{xy_{2f,2opp}}) \\ &= \exp(4.265 + 0 + 0 + 0.381 + 0) \\ &= \exp(4.646) \\ &= 104.217 \end{aligned} \]

\[ \begin{aligned} \hat{\mu}_{223} &= \exp(\lambda + \lambda_{x_{2f}} + \lambda_{y_{2opp}} + \lambda_{z_{lib}} + \lambda_{xy_{2f,2opp}}) \\ &= \exp(4.265 + 0 + 0 + 0 + 0) \\ &= \exp(4.265) \\ &= 71.178 \end{aligned} \]

Fitted Value ini memberikan prediksi banyaknya sebuah kejadian misal prediksi fundamentalisme seorang laki laki yang mendukung hukuman mati adalah 125,59 alias 126. Mendekati nilai aslinya yakni 128 orang.

8.4 Review

Log Linear Model adalah bentuk khusus yang digunakan pada frekeunsi sel tabel kontingensi dan berdistribusi poisson. Pada Model ini dilakukan beberapa kali penentuan signifikansi variabel hingga diperoleh model terbaiknya. Pada tabel kontingensi 2 x 2, akan diperoleh model hingga interaksi antara A dan B yang disebut model saturated. Kemudian diseleksi dengan model tanpa interaksi untuk menentukan model mana yang lebih baik. Pada Tabel kontingensi 2 x 3, sama yakni akan ada interaksi antara A dan B. Kemudian dibandingkan dengan model tanpa interaksi. Jika model tanpa interaksi memberikan hasil yang signifikan maka model tanpa interaksi lebih cocok. Pada tabel kontingensi 3 x 3, akan diperoleh model yang lebih kompleks yakni ada interaksi 2 variabel dan interaksi 3 variabel. Untuk itu dilakukan pengujian saturated(full model) vs homogenous (model tanpa interaksi 3 arah). Jika terdapat variabel yang tidak signifikan, maka dilanjutkan untuk menguji pada interaksi tingkat 2 dan model independent menggunakan selisih nilai deviance. Diakhir dapat dilakukan interpretasi koefisien model untuk melihat kesimpulan model serta dapat dilihat kemampuan model dalam memprediksi menggunakan fitting.

Penutup

Buku ini telah menguraikan berbagai konsep dasar serta penerapan analisis data kategori secara sistematis dan terstruktur. Dimulai dari pengertian data kategori, penyusunan tabel kontingensi, hingga pendekatan inferensial yang relevan, pembaca diarahkan untuk memahami pentingnya analisis ini dalam menjawab pertanyaan-pertanyaan riset yang melibatkan data klasifikasi atau kategorisasi.

Melalui berbagai contoh, penjelasan teori, dan ilustrasi praktis, diharapkan buku ini tidak hanya menambah wawasan konseptual tetapi juga membekali pembaca dengan keterampilan analitis yang aplikatif. Analisis data kategori memiliki kontribusi yang signifikan dalam banyak bidang, karena mampu menjelaskan dinamika hubungan antar variabel yang tidak dapat diukur secara langsung melalui data kuantitatif.

Masih banyak kekurangan buku ini. Namun harapannya buku ini dapat memberikan kontribusi nyata dalam pengembangan literasi statistik di kalangan akademisi maupun praktisi. Pemahaman terhadap data kategori merupakan salah satu langkah penting dalam membentuk pengambilan keputusan yang berbasis data, dan buku ini diharapkan dapat menjadi salah satu rujukan yang terpercaya dalam proses tersebut.

Referensi

Agresti, A. (2013). Categorical Data Analysis (3rd ed.). Hoboken, NJ: Wiley.

Agresti, A. (2019). An Introduction to Categorical Data Analysis (2nd ed.). Hoboken, NJ: Wiley.

Christensen, R. (2005). Categorical Data Analysis Using the SAS System (2nd ed.). Cary, NC: SAS Institute Inc.

Dobson, A. J., & Barnett, A. G. (2008). An Introduction to Generalized Linear Models (3rd ed.). Boca Raton, FL: Chapman & Hall/CRC.

Haryadi, S., & Marwan. (2020). Pengantar Model Linier Umum (Generalized Linear Model). Bandung: Alfabeta.

Kleinbaum, D. G., & Klein, M. (2010). Logistic Regression: A Self-Learning Text (3rd ed.). New York, NY: Springer.

Kurnia, A., & Setyawan, W. (2019). Analisis Log-Linear untuk Data Kategori. Malang: UB Press.

Menard, S. (2002). Applied Logistic Regression Analysis (2nd ed.). Thousand Oaks, CA: SAGE Publications.

Modul ADK, Pertemuan 14. (2021). Universitas Padjadjaran.

Utami, R. P., & Nurkholis, M. (2018). Regresi Multinomial: Teori dan Aplikasi dengan SPSS dan R. Surabaya: Airlangga University Press.

Widodo, A. (2016). Analisis Tabel Kontingensi dalam Penelitian Sosial (Edisi Revisi). Yogyakarta: Pustaka Pelajar.