Pendahuluan

Data kategori adalah tipe data yang terdiri dari kategori atau kelompok, bukan angka yang memiliki nilai kontinu. Data ini banyak digunakan dalam berbagai bidang untuk memahami pola, hubungan, dan tren yang tidak bisa diukur langsung menggunakan skala numerik. Contoh dari data kategori adalah jenis kelamin (laki-laki/perempuan), status pernikahan (menikah/belum menikah), tingkat pendidikan (rendah/menengah/tinggi), dan preferensi konsumen (setuju/netral/tidak setuju).

Tujuan Analisis Data Kategori

Analisis data kategori memiliki beberapa tujuan utama yang membantu dalam pengambilan keputusan berbasis data. Beberapa tujuan utama dari analisis ini meliputi:

Mengidentifikasi Pola Tren

Analisis data kategori memungkinkan peneliti untuk mengidentifikasi pola yang tersembunyi dalam data. Misalnya, dalam riset pemasaran, analisis ini bisa digunakan untuk memahami preferensi konsumen berdasarkan kelompok demografis tertentu.

Menganalisis Hubungan Antarvariabel

Pada banyak penelitian, hubungan antara dua atau lebih variabel kategori menjadi fokus utama. Misalnya, dalam penelitian kesehatan, analisis ini dapat membantu memahami hubungan antara pola hidup dan penyakit tertentu.

Membantu dalam Pengambilan Keputusan

Dengan memahami pola dan hubungan dalam data kategori, pengambil keputusan dapat merancang strategi yang lebih tepat sasaran. Misalnya, dalam kebijakan publik, analisis ini bisa digunakan untuk merancang kebijakan sosial berdasarkan segmentasi masyarakat.

Mengembangkan Model Prediktif

Banyak model prediktif yang menggunakan data kategori sebagai input, seperti regresi logistik untuk memprediksi kejadian biner (misalnya, apakah seseorang memiliki risiko tinggi terhadap penyakit atau tidak).

Definisi dan Ruang Lingkup Analisis Data Kategori

Analisis data kategori adalah teknik statistik yang digunakan untuk menganalisis variabel yang bersifat kategorik.

Nominal vs Ordinal

  • Nominal: Tidak memiliki urutan (contoh: warna, jenis kelamin).
  • Ordinal: Memiliki urutan (contoh: tingkat pendidikan, tingkat kepuasan).

Data Biner vs Multikategori

  • Biner: Memiliki dua kategori (contoh: sukses/gagal, sehat/sakit).
  • Multikategori: Memiliki lebih dari dua kategori.

Perbedaan dengan Data Kuantitatif

Data kategorik berbeda dengan data kuantitatif karena tidak dapat diukur pada skala numerik yang kontinu.

Manfaat Analisis Data Kategori dalam Berbagai Bidang

Analisis data kategori memiliki banyak manfaat dalam berbagai bidang, baik di sektor akademik maupun industri. Berikut adalah beberapa bidang utama di mana analisis ini diterapkan:

Ilmu Sosial dan Psikologi

Dalam ilmu sosial dan psikologi, analisis data kategori digunakan untuk memahami perilaku manusia, opini publik, dan faktor sosial lainnya. Misalnya:

  • Survei kepuasan pelanggan dengan skala Likert (sangat puas hingga sangat tidak puas).
  • Studi mengenai hubungan antara faktor sosial-ekonomi dan tingkat pendidikan seseorang.
  • Penelitian tentang pengaruh terapi psikologis terhadap kelompok pasien tertentu.

Kesehatan dan Kedokteran

Di bidang kesehatan, analisis data kategori sangat penting dalam epidemiologi dan studi klinis. Contoh aplikasinya adalah:

  • Mengkategorikan pasien berdasarkan status kesehatan (sehat/sakit/kronis).
  • Menilai efektivitas pengobatan berdasarkan jenis terapi yang diterima pasien.
  • Menganalisis faktor risiko penyakit berdasarkan gaya hidup dan karakteristik pasien.

Pemasaran dan Bisnis

Dalam pemasaran dan bisnis, data kategori digunakan untuk memahami preferensi pelanggan, segmentasi pasar, dan efektivitas strategi pemasaran. Beberapa aplikasi meliputi:

  • Menentukan preferensi pelanggan terhadap merek tertentu (menyukai/tidak menyukai/netral).
  • Analisis segmentasi pelanggan berdasarkan usia, jenis kelamin, dan lokasi geografis.
  • Studi tentang loyalitas pelanggan berdasarkan tingkat kepuasan layanan.

Pendidikan

Dalam dunia pendidikan, analisis data kategori berguna untuk mengevaluasi metode pengajaran, tingkat kepuasan mahasiswa, dan efektivitas kurikulum. Contohnya:

  • Survei kepuasan mahasiswa terhadap metode pengajaran dosen.
  • Analisis hubungan antara latar belakang sosial-ekonomi dan prestasi akademik.
  • Studi tentang efektivitas program pembelajaran berbasis teknologi.

Kebijakan Publik dan Pemerintahan

Pemerintah sering menggunakan analisis data kategori untuk memahami kebutuhan masyarakat dan merancang kebijakan yang lebih efektif. Beberapa penerapannya antara lain:

  • Analisis tingkat kepuasan masyarakat terhadap layanan publik.
  • Studi tingkat partisipasi masyarakat dalam pemilu.
  • Evaluasi efektivitas program bantuan sosial berdasarkan kategori penerima manfaat.

Keamanan dan Kriminalitas

Di bidang keamanan dan analisis kriminal, data kategori digunakan untuk memahami pola kejahatan dan merancang strategi pencegahan. Contoh aplikasinya adalah:

  • Analisis kategori jenis kejahatan yang paling sering terjadi di suatu wilayah.
  • Studi tentang faktor demografis yang berkorelasi dengan tingkat kriminalitas.
  • Evaluasi efektivitas kebijakan penegakan hukum terhadap berbagai kategori pelanggaran.

Metode dalam Analisis data Kategori

Berbagai metode dapat digunakan dalam analisis data kategori, tergantung pada tujuan penelitian. Beberapa metode umum meliputi:

Tabel Kontingensi dan Uji Chi Square

  • Digunakan untuk menguji hubungan antara dua variabel kategori.
  • Contohnya, menguji apakah ada hubungan antara tingkat pendidikan dan status pekerjaan.

Regresi Logistik

  • Digunakan untuk memprediksi probabilitas suatu kejadian berdasarkan variabel kategori.
  • Contohnya, memprediksi apakah seorang pelanggan akan membeli produk berdasarkan kategori preferensinya.

Analisis Correspondence (CA)

  • Digunakan untuk memahami hubungan antara berbagai kategori dalam satu dataset.
  • Contohnya, analisis tentang preferensi makanan berdasarkan kelompok usia.

Decision Tree dan Random Forest

  • Metode machine learning yang sering digunakan untuk klasifikasi berbasis kategori.
  • Contohnya, mengklasifikasikan pelanggan berdasarkan tingkat loyalitas mereka.

Distribusi Probabilitas dalam Data Kategori

Distribusi Bernoulli

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

Fungsi probabilitasnya:

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

Keterangan Notasi

  • \(X\): Variabel acak biner (0 atau 1)
  • \(p\): Probabilitas sukses (X = 1)

Contoh Variabel Acak Bernoulli

  • Hasil dari lemparan koin (Kepala = 1, Ekor = 0)
  • Keberhasilan atau kegagalan dalam suatu percobaan klinis

Berikut Perhitungan dengan R

library(knitr)
## Warning: package 'knitr' was built under R version 4.3.3
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.3.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
## 
##     group_rows
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
set.seed(123)
bernoulli_sample <- rbinom(n = 10, size = 1, prob = 0.5) # 10 percobaan Bernoulli
bernoulli_sample
##  [1] 0 1 0 1 1 0 1 1 1 0

Distribusi Binomial

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

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

  • \(X\): Jumlah keberhasilan dalam \(n\) percobaan
  • \(n\): Jumlah percobaan
  • \(k\): Jumlah keberhasilan yang diamati
  • \(p\): Probabilitas keberhasilan dalam satu percobaan
  • \(\binom{n}{k}\): Kombinasi “n pilih k”, dihitung sebagai \(\frac{n!}{k!(n-k)!}\)

Contoh Variabel Acak Binomial

  • Jumlah keberhasilan dalam 10 kali lemparan koin
  • Jumlah pasien yang sembuh setelah diberikan obat tertentu dalam suatu studi klinis

Perhitungan dengan R

set.seed(123)
binomial_sample <- rbinom(n = 10, size = 5, prob = 0.5) # 10 percobaan Binomial dengan 5 ulangan
binomial_sample
##  [1] 2 3 2 4 4 1 3 4 3 2

Distribusi Multinomial

Distribusi Multinomial adalah generalisasi lebih lanjut dari distribusi Binomial, digunakan ketika setiap percobaan memiliki lebih dari dua kemungkinan hasil. Jika suatu eksperimen dilakukan 𝑛 kali, dan setiap percobaan dapat menghasilkan salah satu dari 𝑘 kategori dengan probabilitas 𝑝1, 𝑝2, …, 𝑝𝑘, maka distribusi probabilitasnya:

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

Keterangan Notasi

  • \(X_i\): Frekuensi kemunculan kategori ke-\(i\)
  • \(n\): Jumlah total percobaan
  • \(x_i\): Jumlah kejadian kategori ke-\(i\)
  • \(p_i\): Probabilitas kategori ke-\(i\)

Contoh Variabel Acak Multinomial

  • Pemilihan kandidat dalam pemilu (beberapa kandidat, satu suara per pemilih)
  • Distribusi warna permen dalam satu bungkus acak

Berikut Perhitungan dengan R

set.seed(123)
multinomial_sample <- rmultinom(n = 1, size = 10, prob = c(0.3, 0.5, 0.2))
multinomial_sample
##      [,1]
## [1,]    2
## [2,]    5
## [3,]    3

Distribusi Poisson

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

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

Keterangan Notasi

  • \(𝑋\): Jumlah kejadian dalam interval tertentu
  • \(𝜆\): Rata-rata kejadian dalam interval tersebut
  • \(𝑘\): Jumlah kejadian yang diamati

Contoh Variabel Acak Poisson

  • Jumlah panggilan telepon masuk ke pusat layanan dalam satu jam
  • Jumlah kecelakaan lalu lintas di satu jalan dalam sehari

Perhitungan dengan R

set.seed(123)
poisson_sample <- rpois(10, lambda = 3) # 10 sampel dengan rata-rata kejadian 3
poisson_sample
##  [1] 2 4 2 5 6 0 3 5 3 3

Desain Sampling dalam Analisis Data Kategori

Desain sampling sangat penting untuk validitas dan reliabilitas hasil penelitian. Pemilihan desain yang tepat bergantung pada tujuan dan jenis data yang dikumpulkan. Secara umum, desain sampling dalam analisis data kategori terbagi dalam dua pendekatan: prospective sampling dan retrospective sampling.

Prospective Sampling

Metode ini melibatkan pemilihan subjek yang diikuti dalam periode waktu tertentu untuk mengamati perkembangan variabel. Beberapa teknik yang digunakan adalah:

  1. Eksperimen: Subjek dibagi acak ke dalam kelompok perlakuan dan kontrol, menggunakan teknik seperti Simple Random Sampling (SRS), Stratified Random Sampling, dan Cluster Sampling.

  2. Studi Kohort: Mengikuti kelompok individu dari waktu ke waktu. Teknik yang digunakan termasuk Census Sampling, Systematic Sampling, dan Matched Sampling.

Retrospective Sampling

Metode ini mengumpulkan data dari peristiwa yang telah terjadi. Teknik ini sering digunakan untuk menilai hubungan antara faktor risiko dan hasil yang terjadi:

  1. Studi Kasus-Kontrol: Membandingkan kelompok dengan kondisi tertentu dengan yang tanpa kondisi tersebut. Teknik yang digunakan meliputi Purposive Sampling, Snowball Sampling, dan Incidence Density Sampling.

  2. Studi Kohort Retrospektif: Menggunakan data historis untuk mengelompokkan individu dan menganalisis hasilnya. Teknik yang digunakan antara lain Convenience Sampling, Quota Sampling, dan Case-Based Sampling.

Tabel Perbandingan Desain Sampling

Tabel 1: Perbandingan Desain Sampling

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

Tabel Kontingensi 2 x 2

Tabel kontingensi 2 × 2 adalah bentuk paling sederhana dari tabel kontingensi yang digunakan untuk menganalisis hubungan antara dua variabel kategori. Dalam banyak analisis statistik, tabel ini digunakan untuk menentukan apakah terdapat asosiasi antara dua variabel, seperti hubungan antara pengobatan dan hasil kesembuhan, atau kebiasaan merokok dan risiko kanker paru-paru. Tabel kontingensi 2 × 2 adalah bentuk paling sederhana dari tabel kontingensi yang digunakan untuk menganalisis hubungan antara dua variabel kategori. Dalam banyak analisis statistik, tabel ini digunakan untuk menentukan apakah terdapat asosiasi antara dua variabel, seperti hubungan antara pengobatan dan hasil kesembuhan, atau kebiasaan merokok dan risiko kanker paru-paru. Tabel kontingensi 2 × 2 memiliki struktur sebagai berikut:

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

Dengan Keterangan sebagai berikut.

Distribusi Peluang dalam Tabel Kontingensi 2 x 2

Peluang Bersama

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

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

Peluang Marginal

Peluang marginal adalah probabilitas kejadian suatu variabel tanpa mempertimbangkan variabel lainnya:

  • Peluang marginal baris:

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

  • Peluang marginal kolom:

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

Peluang Bersyarat

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

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

Contoh soal perhitungan manual dan R

Berikut contoh kasus dimana ada 2 variabel yaitu cara belajar dan status kelulusan ujian.

Perhitungan dengan R

# Data Observasi
data <- matrix(c(90, 30, 70, 50), nrow = 2, byrow = TRUE)
colnames(data) <- c("Lulus", "Tidak Lulus")
rownames(data) <- c("Belajar Mandiri", "Belajar Kelompok")
n <- sum(data)
# Peluang Bersama
P_joint <- data / n
# Peluang Marginal
P_marginal_rows <- rowSums(data) / n
P_marginal_cols <- colSums(data) / n
# Peluang Bersyarat
P_conditional <- data / rowSums(data)
# Hasil
list(Peluang_Bersama = P_joint, Peluang_Marginal_Baris = P_marginal_rows, Peluang_Marginal_Kolom = P_marginal_cols, Peluang_Bersyarat = P_conditional)
## $Peluang_Bersama
##                      Lulus Tidak Lulus
## Belajar Mandiri  0.3750000   0.1250000
## Belajar Kelompok 0.2916667   0.2083333
## 
## $Peluang_Marginal_Baris
##  Belajar Mandiri Belajar Kelompok 
##              0.5              0.5 
## 
## $Peluang_Marginal_Kolom
##       Lulus Tidak Lulus 
##   0.6666667   0.3333333 
## 
## $Peluang_Bersyarat
##                      Lulus Tidak Lulus
## Belajar Mandiri  0.7500000   0.2500000
## Belajar Kelompok 0.5833333   0.4166667

Interpretasi

  • Peluang Bersama menunjukkan bahwa lebih banyak mahasiswa yang memilih belajar mandiri dan lulus ujian dibandingkan mereka yang belajar dalam kelompok.

  • Peluang Marginal Baris mengindikasikan bahwa distribusi mahasiswa antara belajar mandiri dan belajar kelompok adalah seimbang, masing-masing 50%.

  • Peluang Marginal Kolom menunjukkan bahwa mayoritas mahasiswa lulus ujian (66.67%), sementara 33.33% lainnya tidak lulus.

  • Peluang Bersyarat menunjukkan bahwa mahasiswa yang belajar mandiri memiliki peluang lebih besar untuk lulus ujian (75%) dibandingkan mereka yang belajar dalam kelompok (58.33%). Namun, mahasiswa yang belajar dalam kelompok memiliki peluang lebih besar untuk tidak lulus ujian dibandingkan mereka yang belajar mandiri.

•Peluang bersama menunjukkan probabilitas gabungan dari kejadian tertentu dalam tabel.

• Peluang marginal menunjukkan probabilitas suatu kejadian tanpa mempertimbangkan variabel lain.

• Peluang bersyarat menunjukkan bagaimana probabilitas berubah ketika informasi tentang variabel lain diberikan.

• Tabel kontingensi 2 × 2 membantu dalam memahami hubungan antara dua variabel kategori.

• Peluang bersama, peluang marginal, dan peluang bersyarat adalah dasar untuk memahami hubungan statistik antar variabel dalam tabel kontingensi.

• Perhitungan manual dan implementasi dalam R membantu dalam menganalisis data secara lebih akurat dan efisien.

Ukuran Asosiasi dalam Data Kategori 2 x 2

Risk Difference (RD)

Risk Difference (RD) atau Perbedaan Risiko adalah ukuran absolut dalam epidemiologi yang menggambarkan perbedaan antara probabilitas kejadian suatu hasil dalam dua kelompok yang berbeda. RD dihitung sebagai selisih antara risiko kejadian dalam kelompok terpapar dan risiko kejadian dalam kelompok tidak terpapar.

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

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

Contoh Perhitungan dengan R:

RD <- function(n11, n12, n21, n22) {
(n11 / (n11 + n12)) - (n21 / (n21 + n22))
}
RD(50, 50, 30, 70)
## [1] 0.2

Relative Risk (RR)

Relative Risk (RR) atau Risiko Relatif adalah ukuran yang digunakan dalam epidemiologi untuk membandingkan risiko kejadian suatu peristiwa (misalnya penyakit atau kondisi tertentu) antara dua kelompok, yaitu kelompok yang terpapar dan kelompok yang tidak terpapar. RR menunjukkan seberapa besar kemungkinan kejadian tersebut terjadi pada kelompok terpapar dibandingkan dengan kelompok tidak terpapar.

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

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

Contoh Perhitungan dengan R:

RR <- function(n11, n12, n21, n22) {
(n11 / (n11 + n12)) / (n21 / (n21 + n22))
}
RR(50, 50, 30, 70)
## [1] 1.666667

Odds Ratio (OR)

Odds Ratio (OR) atau Rasio Odds adalah ukuran yang digunakan dalam epidemiologi dan statistik untuk membandingkan odds (peluang) terjadinya suatu kejadian antara dua kelompok, yaitu kelompok yang terpapar dan kelompok yang tidak terpapar. OR sering digunakan dalam studi case-control tetapi juga dapat digunakan dalam studi kohort dan eksperimental.

\[ OR = \frac{n_{11} \times n_{22}}{n_{12} \times n_{21}} \] - Jika \(OR>1\), maka peluang kejadian lebih besar di Grup 1 dibandingkan Grup 2.

  • Jika \(OR<1\), maka peluang kejadian lebih kecil di Grup 1 dibandingkan Grup 2.
  • Jika \(OR=1\), maka tidak ada perbedaan peluang kejadian antara dua kelompok.

Contoh Perhitungan R

OR <- function(n11, n12, n21, n22) {
(n11 * n22) / (n12 * n21)
}
OR(50, 50, 30, 70)
## [1] 2.333333

Perbandingan RD, RR, dan OR

Tabel Ukuran, Desain Sampling, dan Interpretasi

Ukuran Definisi Desain Sampling yang Cocok Interpretasi
Risk Difference (RD) Selisih probabilitas kejadian antara dua kelompok Studi kohort atau eksperimen acak Menunjukkan tambahan atau pengurangan risiko absolut
Relative Risk (RR) Perbandingan risiko kejadian di dua kelompok Studi kohort atau eksperimen klinis Mengukur berapa kali lebih besar risiko di satu kelompok dibandingkan kelompok lainnya
Odds Ratio (OR) Perbandingan odds antara dua kelompok Studi kasus-kontrol atau studi observasional Mengukur hubungan antara variabel eksposur dan kejadian dalam desain studi kasus-kontrol

Kondisi Tepat dalam Penggunaan RD, RR, dan OR

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

Inferensi Tabel Kontingensi Dua Arah

Estimasi

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

Estimasi Titik

Estimasi titik digunakan untuk menentukan satu nilai spesifik sebagai perkiraan terbaik dari parameter populasi.

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

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

Estimasi Interval

Estimasi interval bertujuan untuk memberikan rentang nilai yang diyakini mengandung parameter populasi dengan tingkat kepercayaan tertentu.

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

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

Uji Hipotesis

Uji Proporsi

kontingensi, terutama untuk menentukan apakah terdapat perbedaan yang signifikan dalam proporsi kejadian antara dua kelompok yang berbeda. Tabel kontingensi 2 × 2 memiliki struktur sebagai berikut:

Tabel Kejadian dan Tidak Kejadian

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

Formulasi Uji Proporsi Untuk menguji hipotesis bahwa tidak ada perbedaan proporsi antara dua kelompok, kita menggunakan uji z dua proporsi, dengan hipotesis:

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

Estimasi proporsi dalam masing-masing kelompok diberikan oleh:

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

Estimasi proporsi gabungan (pooling proportion):

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

Statistik uji untuk uji proporsi dua sampel:

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

Contoh Kasus

#Pastikan variabel data_matrix terdefinisi sebelum digunakan
set.seed(123)
data<- matrix(c(50, 30, 30, 50), nrow = 2, byrow = TRUE)
dimnames(data) <- list("Terpapar" = c("Ya", "Tidak"), "Kejadian" = c("Ya", "Tidak"))
print(data)
##         Kejadian
## Terpapar Ya Tidak
##    Ya    50    30
##    Tidak 30    50
# 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 = 9.025, df = 1, p-value = 0.002663
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.08747151 0.41252849
## sample estimates:
## prop 1 prop 2 
##  0.625  0.375

Interpretasi hasil: p-value < 0.05, maka terdapat perbedaan proporsi kejadian antara kelompok terpapar dan tidak terpapar.

Uji Asosiasi

  1. Risk Difference (RD)

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

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

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

  1. Relative Risk

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

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

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

  1. Odds Ratio

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

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

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

Hipotesis \(H_0\) : Tidak ada asosiasi antara dua variabel \(H_1\) : Terdapat asosiasi antara dua variabel

n11 <- 50; n12 <- 30; n21 <- 30; n22 <- 50
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.25
## 
## $SE_RD
## [1] 0.07654655
## 
## $Z_RD
## [1] 3.265986
## 
## $RR
## [1] 1.666667
## 
## $SE_Ln_RR
## [1] 0.1683251
## 
## $Z_RR
## [1] 3.034756
## 
## $OR
## [1] 2.777778
## 
## $SE_Ln_OR
## [1] 0.3265986
## 
## $Z_OR
## [1] 3.128155

Interpretasi Hasil

  • Semua nilai Z-score untuk RD, RR, dan OR lebih besar dari 1.96, yang menunjukkan bahwa perbedaan antara Grup 1 dan Grup 2 dalam hal kejadian positif adalah signifikan secara statistik.
  • Grup 1 lebih berisiko mengalami kejadian positif (misalnya, lulus ujian atau mengalami kejadian tertentu) dibandingkan dengan Grup 2, baik dari perspektif Risk Difference, Relative Risk, maupun Odds Ratio.

Uji Independensi

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

Uji Chi-Square

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

Dimana: - \(O\) adalah observed frequency (frekuensi yang diamati). - \(E\) adalah expected frequency (frekuensi yang diharapkan).

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

Dimana: - \(R_i\) adalah total baris untuk kategori \(i\), - \(C_j\) adalah total kolom untuk kategori \(j\), - \(N\) adalah total keseluruhan sampel.

# Contoh Data
set.seed(123)
data <- matrix(c(30, 10, 15, 45), nrow = 2, byrow = TRUE)
dimnames(data) <- list("Terpapar" = c("Ya", "Tidak"), "Kejadian" = c("Ya", "Tidak"))
print(data)
##         Kejadian
## Terpapar Ya Tidak
##    Ya    30    10
##    Tidak 15    45
# Uji Chi-Square
chisq_test <- chisq.test(data)
print(chisq_test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data
## X-squared = 22.264, df = 1, p-value = 2.376e-06

Interpretasi: Terdapat hubungan antara variabel terpapar dan kejadian

Partisi Chi-Square

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

Langkah-langkah Partisi Chi-Square

  1. Menggunakan tabel I × J dan menghitung Chi-Square keseluruhan.
  2. Memecah tabel menjadi beberapa tabel 2 × 2 sebanyak (I-1)(J-1).
  3. Menghitung statistik Chi-Square pada masing-masing tabel 2 × 2. 4. Menginterpretasikan kategori mana yang memberikan kontribusi signifikan.

Langkah 1

#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,"Reject H0", "Accept H0"))
## $Chi_Square
## X-squared 
##  12.56926 
## 
## $P_Value
## [1] 0.00186475
## 
## $Decision
## [1] "Reject H0"

Interpretasi: Hasil uji Chi-Square menunjukkan bahwa hubungan antara Gender dan Identifikasi Partai Politik signifikan secara statistik

Langkah 2

#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

Interpretasi: Gender lebih mempengaruhi pilihan antara Democrat dan Republican, tetapi tidak berpengaruh signifikan terhadap pilihan Independent.

Uji Likehood Ratio (G\(^{2}\))

\[ G^2 = 2 \sum \sum n_{ij} \ln \left( \frac{n_{ij}}{\\hat{\mu}}} \right) \] Dimana:

  • \(n_{ij}\) adalah frekuensi observasi dalam tabel kontingensi.
  • \(E_{ij}\) adalah frekuensi ekspektasi, dihitung dengan rumus:
  • \(n\) adalah total sampel,
  • \(p_i\) adalah proporsi baris \(i\),
  • \(p_j\) adalah proporsi kolom \(j\).
## CONTOH ##
# Data Observasi
data_matrix <- matrix(c(688, 650, 21, 59), nrow = 2, byrow = TRUE)
colnames(data_matrix) <- c("Cancer (+)", "Control (-)")
rownames(data_matrix) <- c("Smoker", "Non-Smoker")
# Hitung Frekuensi Ekspektasi
data_expected <- chisq.test(data_matrix)$expected
# Hitung Statistik G²
G2 <- 2 * sum(data_matrix * log(data_matrix / data_expected))
29
## [1] 29
# 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", "Failed to Reject H0"))
## $G2
## [1] 19.87802
## 
## $Critical_Value
## [1] 3.841459
## 
## $Decision
## [1] "Reject H0"

Keputusan: ada hubungan

  • Uji Likelihood Ratio (G²) digunakan sebagai alternatif dari uji Chi-Square untuk menguji independensi dalam tabel kontingensi.
  • Jika 𝐺2 lebih besar dari nilai kritis 𝜒2, maka hipotesis nol ditolak, menunjukkan adanya hubungan antara dua variabel.
  • Contoh perhitungan manual dan implementasi di R menunjukkan bahwa ada asosiasi signifikan antara merokok dan kanker paru-paru.

Uji Exact Fisher

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

Keunggulan:

• Cocok untuk ukuran sampel kecil. • Tidak memerlukan asumsi normalitas atau chi-square. • Memberikan hasil yang lebih akurat dibandingkan uji Chi-square pada data dengan frekuensi kecil.

Keterbatasan:

• Perhitungan bisa menjadi sangat berat secara komputasi jika ukuran tabel besar. • Hanya cocok untuk tabel kontingensi kecil (misalnya 2x2 atau 3x3).

Distribusi Hipergeometrik

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

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

Keterangan: - 𝑁 = total objek dalam populasi - 𝐾 = jumlah objek dalam kategori tertentu (misalnya, sukses) - 𝑛 = jumlah sampel yang diambil - 𝑥 = jumlah objek dalam kategori tertentu yang diamati dalam sampel

Contoh

# Definisi parameter
N <- 40 # Total populasi
K <- 29 # Jumlah kategori sukses (bola putih)
n <- 20 # Jumlah sampel diambil
x <- 18 # Jumlah sukses dalam sampel

# Hitung probabilitas P(X = 18)
dhyper(x, m = K, n = N - K, k = n)
## [1] 0.01380413
choose(29, 18) * choose(11, 2) / choose(40, 20)
## [1] 0.01380413
p.value<-0.00007+0.00160+0.01380+0.01380+0.00160+0.00007
p.value
## [1] 0.03094
data <- matrix(c(18, 2, 11, 9), nrow = 2, byrow = TRUE)
fisher.test(data)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  data
## p-value = 0.03095
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   1.147793 78.183838
## sample estimates:
## odds ratio 
##   6.994073

Kesimpulan Materi Sebelumnya

Tabel Uji Statistik

Uji Tujuan Output
Uji Proporsi Menguji apakah proporsi kejadian berbeda antara dua kelompok p-value (signifikansi perbedaan proporsi)
Risk Difference (RD) Mengukur perbedaan absolut risiko antar kelompok Nilai perbedaan risiko
Relative Risk (RR) Mengukur risiko relatif antar kelompok Rasio risiko
Odds Ratio (OR) Mengukur peluang kejadian antar kelompok Rasio odds
Uji Chi-Square Menguji apakah ada hubungan antara dua variabel p-value (signifikansi hubungan)
Uji Fisher’s Exact Menguji hubungan pada sampel kecil p-value (signifikansi hubungan)

Analisis Residual dalam Tabel Kontingensi

Jenis Residual

Pearson Residual

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

dimana:

  • \(O_{ij}\) adalah nilai observasi pada sel ke-\(i,j\),
  • \(E_{ij}\) adalah nilai ekspektasi pada sel ke-\(i,j\).

Standardized Residual (Adjusted Residual)

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

dimana:

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

Contoh Perhitungan Residual R

#Data Observasi
observed <- matrix(c(20, 10, 30, 20), 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)
standardized_residual <- (observed - expected) / sqrt(expected * (1 - row_sum / total_sum) * (1 - col_sum /total_sum))

# Menampilkan hasil
list(Pearson_Residual = pearson_residual,
Standardized_Residual = standardized_residual
)
## $Pearson_Residual
##            [,1]       [,2]
## [1,]  0.2886751 -0.3726780
## [2,] -0.2236068  0.2886751
## 
## $Standardized_Residual
##            [,1]       [,2]
## [1,]  0.5962848 -0.7698004
## [2,] -0.4618802  0.5962848

Deteksi Outlier dalam Analisis Data Kategori Menggunakan Residual

Outlier dalam analisis data kategori adalah sel dalam tabel kontingensi yang memiliki nilai residual yang sangat besar, baik positif maupun negatif.

Menggunakan Residual untuk Mendeteksi Outlier

# Data Observasi
observed <- matrix(c(50, 20, 30, 50), 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)
42
## [1] 42
# 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))

list(
Observed = observed,
Expected = expected,
Pearson_Residual = pearson_residual,
Standardized_Residual = standardized_residual
)
## $Observed
##        Sukses Gagal
## Grup A     50    20
## Grup B     30    50
## 
## $Expected
##          Sukses    Gagal
## Grup A 37.33333 32.66667
## Grup B 42.66667 37.33333
## 
## $Pearson_Residual
##           Sukses     Gagal
## Grup A  2.073070 -2.216205
## Grup B -1.939179  2.073070
## 
## $Standardized_Residual
##           Sukses     Gagal
## Grup A  4.155384 -4.442293
## Grup B -3.887006  4.155384
  • Deteksi outlier menggunakan residual membantu mengidentifikasi kategori dalam tabel kontingensi yang sangat menyimpang dari ekspektasi.
  • Jika |r| > 3, maka sel tersebut dapat dianggap sebagai outlier signifikan yang mungkin mempengaruhi kesimpulan analisis.
  • Residual besar menunjukkan hubungan signifikan antara variabel, sedangkan residual mendekati nol menunjukkan tidak adanya hubungan.

Tabel Kontingensi Tiga Arah

Tabel Parsial dan Marginal

  • Tabel Parsial: Tabel yang menyajikan hubungan antara 𝑋 dan 𝑌 pada setiap kategori 𝑍. Ini memungkinkan analisis hubungan bersyarat, di mana efek 𝑍 dikendalikan.
  • Tabel Marginal: Tabel yang diperoleh dengan mengabaikan 𝑍, yaitu dengan menjumlahkan semua kategori 𝑍. Ini memberikan gambaran umum hubungan antara 𝑋 dan 𝑌 tanpa mempertimbangkan 𝑍, yang bisa menyebabkan distorsi interpretasi, seperti dalam Simpson’s Paradox.

Contoh

data3 <- array(c(10, 40, 50, 5, 10, 20, 40, 10, 0, 45, 40, 30),
dim = c(2, 2, 3),
dimnames = list(
Merokok = c("Ya", "Tidak"),
Kanker = c("Ya", "Tidak"),
Usia = c("Muda", "Dewasa", "Tua")
))
# Ekstrak tabel parsial berdasarkan usia
freq_parsial_muda <- data3[, , "Muda"]
freq_parsial_dewasa <- data3[, , "Dewasa"]
freq_parsial_tua <- data3[, , "Tua"]
# Tampilkan hasil
freq_parsial_muda
##        Kanker
## Merokok Ya Tidak
##   Ya    10    50
##   Tidak 40     5
# Hitung frekuensi marginal
freq_marginal_X <- apply(data3, 1, sum)
freq_marginal_Z <- apply(data3, 3, sum)
# Tampilkan hasil
freq_marginal_X
##    Ya Tidak 
##   150   150

Kesimpulan: - Tabel frekuensi parsial memungkinkan analisis hubungan antara dua variabel dengan mempertimbangkan variabel kontrol. Ini membantu mengidentifikasi pola yang mungkin tersembunyi dalam tabel frekuensi marginal. - Tabel frekuensi marginal memberikan gambaran umum tentang distribusi kategori variabel dalam tabel kontingensi tiga arah. Ini berguna dalam memahami pola distribusi data tanpa mempertimbangkan hubungan antarvariabel.

Distribusi Peluang

  1. Peluang Bersama (Joint Probability)

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

data <- data.frame(
Usia = rep(c("Muda", "Dewasa", "Tua"), each = 2),
Merokok = rep(c("Ya", "Tidak"), 3),
Kanker = c(10, 5, 40, 10, 50, 20),
Total = c(50, 50, 50, 50, 50, 50)
)

data$P_ZXY <- data$Kanker / sum(data$Total)
data
data3 <- array(c(10, 40, 50, 5, 10, 20, 40, 10, 0, 45, 40, 30),
dim = c(2, 2, 3),
dimnames = list(
Merokok = c("Ya", "Tidak"),
Kanker = c("Ya", "Tidak"),
Usia = c("Muda", "Dewasa", "Tua")
))

# Hitung probabilitas bersama
total <- sum(data3)
joint_prob <- data3 / total
ftable(joint_prob)
##                Usia       Muda     Dewasa        Tua
## Merokok Kanker                                      
## Ya      Ya          0.03333333 0.03333333 0.00000000
##         Tidak       0.16666667 0.13333333 0.13333333
## Tidak   Ya          0.13333333 0.06666667 0.15000000
##         Tidak       0.01666667 0.03333333 0.10000000
  1. Peluang Marginal (Marginal Probability)
P_Y1 <- sum(data$Kanker) / sum(data$Total)
P_Y1
## [1] 0.45
data3 <- array(c(10, 40, 50, 5, 10, 20, 40, 10, 0, 45, 40, 30),
dim = c(2, 2, 3),
dimnames = list(
Merokok = c("Ya", "Tidak"),
Kanker = c("Ya", "Tidak"),
Usia = c("Muda", "Dewasa", "Tua")
))
# Hitung probabilitas bersama
total <- sum(data3)
joint_prob <- data3 / total
# Hitung probabilitas marginal
marginal_X <- apply(joint_prob, 1, sum)
marginal_Z <- apply(joint_prob, 3, sum)
# Tampilkan hasil
marginal_X
##    Ya Tidak 
##   0.5   0.5

Tabel Peluang Bersyarat

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

data3 <- array(c(10, 40, 50, 5, 10, 20, 40, 10, 0, 45, 40, 30),
dim = c(2, 2, 3),
dimnames = list(
Merokok = c("Ya", "Tidak"),
Kanker = c("Ya", "Tidak"),
Usia = c("Muda", "Dewasa", "Tua")
))
ftable(data3)
##                Usia Muda Dewasa Tua
## Merokok Kanker                     
## Ya      Ya            10     10   0
##         Tidak         50     40  40
## Tidak   Ya            40     20  45
##         Tidak          5     10  30

Ukuran Asosiasi

  1. Risk Difference (RD)

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

  1. RR (Risk Ratio)

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

  1. OR (Odds Ratio)

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

data <- matrix(c(10, 40, 5, 45), nrow = 2, byrow = TRUE)
rownames(data) <- c("Merokok", "Tidak")
colnames(data) <- c("Kanker", "Tidak Kanker")
data
##         Kanker Tidak Kanker
## Merokok     10           40
## Tidak        5           45
#Hitung Risk Difference
p1 <- data[1, 1] / sum(data[1, ])
p2 <- data[2, 1] / sum(data[2, ])
RD <- p1 - p2
RD
## [1] 0.1
# Hitung Relative Risk
RR <- p1 / p2
RR
## [1] 2
#Hitung Odds Ratio
odds1 <- data[1, 1] / data[1, 2]
odds2 <- data[2, 1] / data[2, 2]
OR <- odds1 / odds2
OR
## [1] 2.25

Hasil analisis menunjukkan bahwa:

• Risk Difference sebesar 0.1 (10%) mengindikasikan bahwa perokok memiliki peluang terkena kanker paru-paru yang lebih tinggi sebesar 10% dibandingkan non-perokok. • Relative Risk sebesar 2.0 menunjukkan bahwa perokok memiliki risiko 2 kali lebih tinggi terkena kanker paru-paru dibandingkan non-perokok. • Odds Ratio sebesar 2.25 menunjukkan bahwa odds terkena kanker paru-paru pada perokok adalah 2.25 kali lebih besar dibandingkan dengan non-perokok.

Struktur dalam Tabel Kontingensi 3D

data_tabel <- data.frame(
"Victims' Race" = c("White", "", "Black", "", "Total", ""),
"Defendant's Race" = c("White", "Black", "White", "Black", "White", "Black"),
"Yes" = c(53, 11, 0, 4, 53, 15),
"No" = c(414, 37, 16, 139, 430, 176),
"Percentage Yes" = c(11.3, 22.9, 0.0, 2.8, 11.0, 7.9)
)
kable(data_tabel, format = "latex", booktabs = TRUE,
caption = "Putusan Hukuman Mati berdasarkan Ras Tersangka dan Korban") %>%
kable_styling(full_width = FALSE, position = "center")

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

\[ P(X, Y | Z) = P(X | Z) P(Y | Z) \]

Inferensi Tabel Kontingensi Tiga Arah

Tabel kontingensi tiga arah digunakan untuk menganalisis hubungan antara dua variabel kategorik dengan mempertimbangkan variabel kontrol. Contohnya adalah hubungan antara kebiasaan merokok (X) dan kanker paru-paru (Y ) dengan variabel kontrol usia (Z).

Jika odds ratio parsial relatif konstan, kita dapat menghitung odds ratio bersama menggunakan estimasi Mantel-Haenszel.

Independensi Bersyarat dalam Tabel Kontingensi Tiga Arah

Independensi Bersyarat: Dua variabel, 𝑋 dan 𝑌, dikatakan independen bersyarat terhadap variabel ketiga, 𝑍, jika rasio odds mereka dalam setiap strata 𝑍 sama dengan 1.

Secara matematis, ini dapat dituliskan sebagai:

\[𝑅(𝑋, 𝑌 |𝑍) = 1\]

Artinya, setelah mengendalikan pengaruh 𝑍, tidak ada hubungan antara 𝑋 dan 𝑌 dalam setiap strata.

Pengujian Statistik untuk Independensi Bersyarat

  1. UJI CMH Metode Cochran-Mantel-Haenszel Tujuan Uji CMH Uji Cochran-Mantel-Haenszel (CMH) digunakan untuk menguji hubungan antara dua variabel kategori dengan mempertimbangkan efek dari variabel perancu (confounder). Uji ini berguna dalam:

• Menguji hubungan yang dikendalikan oleh faktor perancu dalam tabel kontingensi berlapis (stratified contingency tables). • Menguji hipotesis independensi antara dua variabel kategori dengan mempertimbangkan efek dari variabel ketiga. • Mengatasi bias akibat confounding, misalnya dalam studi epidemiologi atau eksperimen sosial.

data_cmh <- array(c(50, 30, 20, 50, # Laki-laki
40, 20, 10, 40), # Perempuan
dim = c(2, 2, 2),
dimnames = list(
Merokok = c("Perokok", "Tidak Merokok"),
Kanker = c("Kanker (+)", "Kanker (-)"),
Jenis_Kelamin = c("Laki-laki", "Perempuan")
))
data_cmh
## , , Jenis_Kelamin = Laki-laki
## 
##                Kanker
## Merokok         Kanker (+) Kanker (-)
##   Perokok               50         20
##   Tidak Merokok         30         50
## 
## , , Jenis_Kelamin = Perempuan
## 
##                Kanker
## Merokok         Kanker (+) Kanker (-)
##   Perokok               40         10
##   Tidak Merokok         20         40
cmh_base <- mantelhaen.test(data_cmh, correct = FALSE)
cmh_base
## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  data_cmh
## Mantel-Haenszel X-squared = 39.86, df = 1, p-value = 2.729e-10
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  3.132410 9.187417
## sample estimates:
## common odds ratio 
##          5.364583
#library(vcdExtra) # Untuk uji CMH

# Contoh Data: Tabel Kontingensi 2 × 2 × k
data <- array(c(20, 30, 10, 40,
15, 25, 5, 35,
18, 32, 12, 38),
dim = c(2, 2, 3),
dimnames = list(
X = c("Ya", "Tidak"),
Y = c("Ya", "Tidak"),
Z = c("Stratum 1", "Stratum 2", "Stratum 3")
))
# Melakukan uji Cochran-Mantel-Haenszel
mantelhaen.test(data)
## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  data
## Mantel-Haenszel X-squared = 10.848, df = 1, p-value = 0.000989
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.482202 4.377281
## sample estimates:
## common odds ratio 
##          2.547159

Jika nilai p lebih kecil dari 0.05, maka kita menolak hipotesis nol, yang berarti ada hubungan signifikan antara merokok dan kanker paru-paru setelah mengontrol jenis kelamin.

Untuk menguji apakah terjadi independensi bersyarat antara kebiasaan merokok dan kanker paru-paru dengan mempertimbangkan provinsi sebagai kontrol, kita dapat menggunakan uji Cochran-Mantel-Haenszel (CMH).

Odds Ratio Bersama

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

Interpretasi Hasil:

• Odds ratio bersama sebesar 2.17 menunjukkan bahwa perokok memiliki risiko lebih tinggi terkena kanker paru-paru dibandingkan non-perokok. • Standard error sebesar 0.04 menunjukkan variabilitas kecil pada estimasi odds ratio. • Interval kepercayaan 95% (1.98 hingga 2.38) menunjukkan bahwa kita cukup yakin bahwa odds ratio sebenarnya berada dalam rentang ini. • Karena interval kepercayaan tidak mencakup 1, hasil ini signifikan secara statistik, mendukung adanya hubungan antara merokok dan kanker paru-paru.

Kesimpulan:

Odds ratio bersama digunakan untuk menyimpulkan hubungan antara dua variabel setelah mempertimbangkan efek dari variabel ketiga. Statistik Mantel-Haenszel adalah metode yang andal untuk mengestimasi dan menguji hubungan ini.

Model Regresi Poisson

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

Misalkan kita memiliki data simulasi sebagai berikut:

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

Estimasi Regresi Poisson

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.21817    0.06712   3.250  0.00115 ** 
## x            0.58748    0.06288   9.343  < 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: 333.72  on 199  degrees of freedom
## Residual deviance: 244.46  on 198  degrees of freedom
## AIC: 587.36
## 
## Number of Fisher Scoring iterations: 5

Plot Hasil Prediksi

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)

Diagnostik dan Overdispersion

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

Interpretasi 1. Model Regresi Poisson: Intercept: 0.21817 (signifikan pada level 0.01). x: 0.58748 (signifikan pada level 0.001). Artinya, untuk setiap kenaikan satu unit pada x, logaritma ekspektasi dari y (yaitu log(λ)) meningkat sebesar 0.58748.

  1. Signifikansi Koefisien: Kedua koefisien (Intercept dan x) sangat signifikan

  2. Deviasi: Pengurangan deviasi menunjukkan bahwa model dengan prediktor x lebih baik dibandingkan model tanpa prediktor.

  3. Dispersi: Nilai dispersion: 1.144237. Nilai ini sedikit lebih besar dari 1, yang menunjukkan bahwa data tidak sepenuhnya mengikuti distribusi Poisson murni. Ini mengindikasikan adanya sedikit overdispersion (variabilitas yang lebih besar dari yang diharapkan).

  4. Plot:

Grafik yang dihasilkan menunjukkan hubungan antara x dan y (data aktual) dengan garis prediksi dari model regresi Poisson.

Inferensi GLM

Mencari Ekspektasi dan Varians dalam GLM

  1. Ekspektasi Estimator Ekspektasi menunjukkan apakah suatu estimator tak bias, yaitu: 𝔼[𝛽]̂ = 𝛽 Dalam GLM, MLE dari 𝛽̂ bersifat asymptotically unbiased.

  2. Varians Estimator Varians menunjukkan presisi dari estimasi parameter: Var(𝛽)̂ ≈ [X⊤WX]−1 di mana W adalah matriks bobot yang tergantung pada distribusi dan fungsi link.

Detail Metode Estimasi dan Inferensi Regresi Logistik

Estimasi MLE dengan Newton-Raphson (Manual di R)

#Data simulasi
set.seed(1)
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 # hasil estimasi akhir
##        [,1]
##   -1.516679
## x  2.155658

Inferensi Parameter

  1. Uji Wald
set.seed(123)
n <- 100
x <- rnorm(n)
log_odds <- -0.5 + 1.2 * x
p <- 1 / (1 + exp(-log_odds))
y <- rbinom(n, 1, p)
data <- data.frame(x, y)
model <- glm(y ~ x, data = data, family = binomial)
summary(model)
## 
## Call:
## glm(formula = y ~ x, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3097     0.2296  -1.349    0.177    
## x             1.2663     0.3080   4.111 3.94e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 137.99  on 99  degrees of freedom
## Residual deviance: 114.76  on 98  degrees of freedom
## AIC: 118.76
## 
## Number of Fisher Scoring iterations: 4
#Langkah 1: Ambil nilai koefisien dan SE
beta_hat <- coef(model)["x"]
se_beta <- summary(model)$coefficients["x", "Std. Error"]

#Langkah 2: Hitung Statistik Z
Z <- beta_hat / se_beta
Z
##        x 
## 4.110965
#Langkah 3: Hitung Statistik Wald
Wald_stat <- Z^2
Wald_stat
##        x 
## 16.90003
#Langkah 4: Hitung p-value
p_value <- 1 - pchisq(Wald_stat, df = 1)
p_value
##            x 
## 3.940095e-05
  1. Uji Likelihood Ratio (Chi-Square)
#Bandingkan model penuh dengan model tanpa prediktor.

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

Evaluasi Kebaikan Model 1. Akaike Information Criterion (AIC) Semakin kecil AIC, semakin baik model.

AIC(model)
## [1] 118.7598
# Bayesian Information Criterion (BIC)
BIC(model)
## [1] 123.9701

Detail Metode Estimasi dan Inferensi Regresi Poisson

#Simulasi Data

set.seed(123)
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)

#IRLS Manual Step-by-Step

# Inisialisasi
beta <- c(0, 0)
tol <- 1e-6
max_iter <- 100
for (i in 1:max_iter) {
eta <- X %*% beta
lambda <- exp(eta)
W <- diag(as.numeric(lambda))
z <- eta + (y - lambda) / lambda
beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
if (sum(abs(beta_new - beta)) < tol) {
cat("Konvergen pada iterasi ke-", i, "\n")
break
}
beta <- beta_new
}
## Konvergen pada iterasi ke- 8
beta # hasil estimasi
##        [,1]
##   0.4494951
## x 0.8600048
#Perbandingan dengan glm()
model_glm <- glm(y ~ x, family = poisson)
summary(model_glm)
## 
## Call:
## glm(formula = y ~ x, family = poisson)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.44950    0.08872   5.066 4.05e-07 ***
## x            0.86000    0.07463  11.523  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.05  on 99  degrees of freedom
## Residual deviance: 106.78  on 98  degrees of freedom
## AIC: 325.76
## 
## Number of Fisher Scoring iterations: 5
#Pengujian hipotesis Uji Wald
#Untuk menguji H0: = 0

# 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.110965 
## Chi-Square: 16.90003 
## p-value: 3.940095e-05
#Uji Likelihood Ratio (Chi-Square)
model_null <- glm(y ~ 1, family = poisson, data = data)
anova(model_null, model, test = "Chisq")
#Evaluasi Model (AIC & BIC)
AIC(model)
## [1] 118.7598
BIC(model)
## [1] 123.9701

Contoh Kasus

Data yang disajikan adalah hasil dari sebuah studi yang menganalisis hubungan antara frekuensi penggunaan alat kesehatan online dan jumlah kunjungan ke ruang gawat darurat (UGD) pada pasien dengan penyakit tertentu, seperti Penyakit Paru Obstruktif Kronis (COPD). Berikut adalah penjelasan terkait variabel yang ada dalam data:

  1. Frekuensi Penggunaan Alat Online:
    • Jarang (≤1x/minggu): Pasien yang menggunakan alat kesehatan online dengan frekuensi lebih rendah atau kurang dari satu kali per minggu.
    • Kadang (1–3x/minggu): Pasien yang menggunakan alat kesehatan online sekitar 1 hingga 3 kali per minggu.
    • Sering (≥4x/minggu): Pasien yang menggunakan alat kesehatan online lebih sering, yaitu 4 kali atau lebih dalam seminggu.
  2. Jumlah Kunjungan ke UGD:
    • 0 Kunjungan: Pasien yang tidak mengunjungi UGD dalam periode waktu yang diteliti.
    • 1–3 Kunjungan: Pasien yang mengunjungi UGD antara 1 hingga 3 kali dalam periode waktu yang diteliti.
    • >3 Kunjungan: Pasien yang mengunjungi UGD lebih dari 3 kali dalam periode waktu yang diteliti.

Struktur Data

Dalam tabel ini, kita melihat hubungan antara dua variabel kategorikal: - Frekuensi Penggunaan Alat Online di sisi baris - Jumlah Kunjungan ke UGD di sisi kolom

Setiap sel dalam tabel menggambarkan jumlah pasien yang berada dalam kombinasi tertentu dari kedua variabel tersebut. Misalnya, 12 pasien dengan frekuensi penggunaan alat kesehatan online “Jarang (≤1x/minggu)” tidak mengunjungi UGD sama sekali (0 kunjungan), sedangkan 55 pasien dalam kategori yang sama mengunjungi UGD 1–3 kali.

Tujuan Studi

Studi ini bertujuan untuk mengevaluasi apakah ada hubungan antara seringnya penggunaan alat kesehatan online dengan jumlah kunjungan ke UGD. Dengan kata lain, studi ini ingin mengetahui apakah pasien yang lebih sering menggunakan alat kesehatan online cenderung lebih sedikit atau lebih banyak mengunjungi UGD, atau jika tidak ada hubungan antara keduanya.

Tujuan Analisis

  • Analisis Tabel Kontingensi 3x3 digunakan untuk menilai apakah ada hubungan yang signifikan antara frekuensi penggunaan alat kesehatan online dan jumlah kunjungan ke UGD.
  • Hasil analisis ini dapat memberikan wawasan tentang apakah penggunaan alat kesehatan online bisa membantu mengurangi frekuensi kunjungan ke UGD atau jika tidak ada pengaruh yang signifikan.

Dengan menggunakan uji statistik seperti Chi-Square, Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR), kita dapat menguji apakah ada hubungan yang signifikan antara dua variabel ini dan seberapa besar pengaruhnya.

# Membuat data kontingensi 3x3
data <- matrix(c(12, 55, 100, 
                 21, 37, 19, 
                 105, 11, 15), 
               nrow = 3, 
               byrow = TRUE)

# Menyimpan nama kolom dan baris
colnames(data) <- c("0 Kunjungan", "1-3 Kunjungan", ">3 Kunjungan")
rownames(data) <- c("Jarang (≤1x/minggu)", "Kadang (1–3x/minggu)", "Sering (≥4x/minggu)")

# Menampilkan data
data
##                      0 Kunjungan 1-3 Kunjungan >3 Kunjungan
## Jarang (≤1x/minggu)           12            55          100
## Kadang (1–3x/minggu)          21            37           19
## Sering (≥4x/minggu)          105            11           15
# Menghitung peluang bersama (joint probabilities)
prob_joint <- data / sum(data)
prob_joint
##                      0 Kunjungan 1-3 Kunjungan >3 Kunjungan
## Jarang (≤1x/minggu)        0.032    0.14666667   0.26666667
## Kadang (1–3x/minggu)       0.056    0.09866667   0.05066667
## Sering (≥4x/minggu)        0.280    0.02933333   0.04000000
# Peluang marginal (probabilitas setiap kategori)
prob_marginal <- margin.table(data, 1) / sum(data)  # Baris
prob_marginal_column <- margin.table(data, 2) / sum(data)  # Kolom

# Peluang bersyarat (conditional probabilities)
prob_conditional <- sweep(prob_joint, 2, prob_marginal_column, "/")
prob_conditional
##                      0 Kunjungan 1-3 Kunjungan >3 Kunjungan
## Jarang (≤1x/minggu)   0.08695652     0.5339806    0.7462687
## Kadang (1–3x/minggu)  0.15217391     0.3592233    0.1417910
## Sering (≥4x/minggu)   0.76086957     0.1067961    0.1119403
# Menghitung Risk Difference (RD)
rd <- (data[2,2] / sum(data[2,])) - (data[1,2] / sum(data[1,]))
rd
## [1] 0.1511782
# Menghitung Relative Risk (RR)
rr <- (data[2,2] / sum(data[2,])) / (data[1,2] / sum(data[1,]))
rr
## [1] 1.459032
# Menghitung Odds Ratio (OR)
or <- (data[2,2] * data[1,1]) / (data[1,2] * data[2,1])
or
## [1] 0.3844156
# Uji Chi-Square
chi_square <- chisq.test(data)
chi_square$p.value
## [1] 3.005891e-40
# Hipotesis Uji Chi-Square
# H0: Tidak ada hubungan antara frekuensi penggunaan alat online dan kunjungan ke UGD.
# H1: Ada hubungan antara frekuensi penggunaan alat online dan kunjungan ke UGD.

if (chi_square$p.value < 0.05) {
  print("H0 ditolak: Ada hubungan antara penggunaan alat online dan kunjungan ke UGD.")
} else {
  print("Gagal menolak H0: Tidak ada hubungan antara penggunaan alat online dan kunjungan ke UGD.")
}
## [1] "H0 ditolak: Ada hubungan antara penggunaan alat online dan kunjungan ke UGD."
# Analisis residual
residuals <- chi_square$residuals
residuals
##                      0 Kunjungan 1-3 Kunjungan >3 Kunjungan
## Jarang (≤1x/minggu)    -6.308656      1.348159     5.220150
## Kadang (1–3x/minggu)   -1.378130      3.446666    -1.623250
## Sering (≥4x/minggu)     8.179516     -4.164635    -4.649437
# Menyusun model GLM
data <- matrix(c(12, 55, 100, 21, 37, 19, 105, 11, 15), nrow = 3, byrow = TRUE)
rownames(data) <- c("Jarang (≤1x/minggu)", "Kadang (1–3x/minggu)", "Sering (≥4x/minggu)")
colnames(data) <- c("0 Kunjungan", "1–3 Kunjungan", ">3 Kunjungan")
# Membuat data long untuk analisis multinomial
data_long <- data.frame(
  outcome = factor(rep(c("0 Kunjungan", "1–3 Kunjungan", ">3 Kunjungan"),
                       times = c(sum(data[1,]), sum(data[2,]), sum(data[3,])))),
  predictor = rep(c("Jarang", "Kadang", "Sering"), each = sum(data)),
  count = as.vector(data)
)
# Memasang paket nnet untuk regresi multinomial
#install.packages("nnet")
library(nnet)
## Warning: package 'nnet' was built under R version 4.3.3
# Model multinomial regresi untuk tabel 3x3
glm_model <- multinom(outcome ~ predictor, data = data_long)
## # weights:  12 (6 variable)
## initial  value 1235.938825 
## final  value 1184.305279 
## converged
summary(glm_model)
## Call:
## multinom(formula = outcome ~ predictor, data = data_long)
## 
## Coefficients:
##               (Intercept) predictorKadang predictorSering
## 0 Kunjungan     0.2427966    1.252644e-06    1.252648e-06
## 1–3 Kunjungan  -0.5313873   -2.811809e-06   -2.811815e-06
## 
## Std. Errors:
##               (Intercept) predictorKadang predictorSering
## 0 Kunjungan     0.1167117       0.1650553       0.1650553
## 1–3 Kunjungan   0.1435986       0.2030792       0.2030792
## 
## Residual Deviance: 2368.611 
## AIC: 2380.611

Interpretasi Hasil:

  • Peluang Bersama dan Peluang Marginal memberikan gambaran umum tentang seberapa sering kombinasi kategori terjadi dalam dataset.

  • Peluang Bersyarat memberikan gambaran bagaimana setiap kategori terkait dengan kategori lainnya, berdasarkan distribusi dalam kolom.

  • Risk Difference menunjukkan bahwa kelompok yang lebih sering menggunakan alat online (Kadang) memiliki risiko lebih tinggi dalam jumlah kunjungan dibandingkan dengan kelompok yang jarang menggunakannya.

  • Relative Risk menunjukkan bahwa kelompok yang lebih sering menggunakan alat online memiliki risiko jauh lebih tinggi untuk mengunjungi UGD dibandingkan dengan yang jarang menggunakannya.

  • Odds Ratio mengindikasikan bahwa odds kunjungan ke UGD lebih tinggi untuk kelompok yang lebih sering menggunakan alat online.

  • Uji Chi-Square menunjukkan hubungan yang signifikan antara penggunaan alat online dan jumlah kunjungan ke UGD.

Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio

Regresi logistik merupakan salah satu metode analisis statistik yang digunakan untuk memodelkan hubungan antara satu variabel respons biner (dua kategori) dengan satu atau lebih variabel prediktor. Dalam penerapannya, prediktor yang digunakan dapat memiliki skala pengukuran berbeda, yaitu:

– Sebagai nominal dengan membuat dummy untuk setiap kategori.

– Atau diperlakukan seperti rasio dengan memberi nilai tertentu (misalnya SMA = 1, Sarjana = 2, Master = 3, Doktor = 4) untuk mencerminkan urutan dan jarak antar tingkatannya.

Simulasi Data

Data terdiri dari 500 observasi yang merepresentasikan individu dengan berbagai karakteristik perjalanan. Variabel respons berupa keputusan membeli asuransi (ya atau tidak), sementara variabel prediktor terdiri atas campuran skala: nominal (jenis perjalanan: domestik atau internasional), ordinal (tingkat kepanikan terhadap risiko perjalanan: Low, Medium, High, Very High), dan rasio (durasi perjalanan dalam hari). Simulasi ini memungkinkan analisis hubungan antara karakteristik perjalanan dan keputusan pembelian asuransi dengan pendekatan statistik yang fleksibel dan dapat dikendalikan.

library(tibble)
## Warning: package 'tibble' was built under R version 4.3.3
set.seed(2025)
n <- 500

# Nominal: Jenis perjalanan
travel_type <- sample(c("Domestic", "International"), n, replace = TRUE)

# Ordinal: Tingkat kepanikan terhadap risiko
risk_level <- sample(c("Low", "Medium", "High", "Very High"), n, 
                     replace = TRUE, prob = c(0.3, 0.4, 0.2, 0.1))

# Rasio: Durasi perjalanan dalam hari
trip_duration <- rpois(n, lambda = 7) + 1  # durasi minimal 1 hari

# Konversi risk_level jadi nilai ordinal numerik
risk_numeric <- as.numeric(factor(risk_level, 
                         levels = c("Low", "Medium", "High", "Very High"), 
                         ordered = TRUE))

# Model logit
logit_p <- -2 + 0.6 * (travel_type == "International") + 
               0.5 * risk_numeric + 
               0.1 * trip_duration

# Probabilitas membeli asuransi
p <- 1 / (1 + exp(-logit_p))

# Respons: membeli atau tidak
buy_insurance <- rbinom(n, 1, p)

# Gabungkan ke dalam data frame
sim_data <- tibble::tibble(buy_insurance, travel_type, risk_level, trip_duration)
head(sim_data)

Eksplorasi Data

library(dplyr)

sim_data %>% 
  dplyr::group_by(buy_insurance) %>% 
  dplyr::summarise(Jumlah = n(), Rata2_Durasi = mean(trip_duration))

Interpretasi Jumlah individu yang membeli dan tidak membeli asuransi adalah sama (masing-masing 250). Rata-rata durasi perjalanan bagi yang membeli asuransi (8.436 hari) sedikit lebih panjang dibandingkan yang tidak membeli (7.584 hari), yang mengindikasikan bahwa semakin lama perjalanan, kemungkinan membeli asuransi cenderung meningkat.

Perlakuan Variabel Ordinal

Perlakuan sebagai Nominal

sim_data_nominal <- sim_data %>%
  mutate(
    risk_level = factor(risk_level, levels = c("Low", "Medium", "High", "Very High"))
  )

model_nominal <- glm(buy_insurance ~ travel_type + risk_level + trip_duration, 
                     data = sim_data_nominal, family = binomial)
summary(model_nominal)
## 
## Call:
## glm(formula = buy_insurance ~ travel_type + risk_level + trip_duration, 
##     family = binomial, data = sim_data_nominal)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -1.77335    0.35021  -5.064 4.11e-07 ***
## travel_typeInternational  0.77484    0.19125   4.051 5.09e-05 ***
## risk_levelMedium          0.35464    0.22307   1.590  0.11187    
## risk_levelHigh            0.81989    0.27575   2.973  0.00295 ** 
## risk_levelVery High       1.81932    0.37303   4.877 1.08e-06 ***
## trip_duration             0.11289    0.03603   3.133  0.00173 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 693.15  on 499  degrees of freedom
## Residual deviance: 633.92  on 494  degrees of freedom
## AIC: 645.92
## 
## Number of Fisher Scoring iterations: 4

Interpretasi

Intercept: Log-odds membeli asuransi untuk referensi (perjalanan domestik, risiko Low, durasi = 0), yaitu -1.773.

travel_typeInternational: Peluang membeli asuransi lebih tinggi secara signifikan untuk perjalanan internasional (OR ≈ exp(0.775) ≈ 2.17).

risk_level:

Medium: Tidak signifikan.

High: Signifikan (p < 0.01), peluang meningkat.

Very High: Sangat signifikan (p < 0.001), peluang membeli meningkat tajam.

trip_duration: Semakin lama perjalanan, peluang membeli asuransi meningkat signifikan.

Perlakuan Risk Level sebagai Numeric

sim_data_numeric <- sim_data %>%
  mutate(
    risk_numeric = case_when(
      risk_level == "Low" ~ 1,
      risk_level == "Medium" ~ 2,
      risk_level == "High" ~ 3,
      risk_level == "Very High" ~ 4
    )
  )

model_numeric <- glm(buy_insurance ~ travel_type + risk_numeric + trip_duration,
                     data = sim_data_numeric, family = binomial)
summary(model_numeric)
## 
## Call:
## glm(formula = buy_insurance ~ travel_type + risk_numeric + trip_duration, 
##     family = binomial, data = sim_data_numeric)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -2.3974     0.3873  -6.190 6.00e-10 ***
## travel_typeInternational   0.7671     0.1903   4.031 5.55e-05 ***
## risk_numeric               0.5292     0.1014   5.219 1.80e-07 ***
## trip_duration              0.1136     0.0359   3.165  0.00155 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 693.15  on 499  degrees of freedom
## Residual deviance: 635.97  on 496  degrees of freedom
## AIC: 643.97
## 
## Number of Fisher Scoring iterations: 4

Interpretasi

Perlakuan ordinal menghasilkan koefisien yang signifikan (p < 0.001) → semakin tinggi tingkat risiko (Low sampai Very High), log-odds membeli asuransi meningkat secara linier.

Hal ini lebih cocok jika tingkatan risiko memang teratur secara logika dan pengaruhnya juga bertingkat.

Perbandingan Model

nullmod <- glm(buy_insurance ~ 1, data = sim_data, family = binomial)

# AIC
AIC(model_nominal)
## [1] 645.9207
AIC(model_numeric)
## [1] 643.9686

Interpretasi

Model dengan risk level sebagai numeric memiliki nilai AIC lebih rendah, artinya model ini lebih baik secara keseluruhan (lebih fit dan lebih parsimonious) dibandingkan pendekatan nominal.

Goodness-of-Fit

1 - (logLik(model_nominal)/logLik(nullmod))
## 'log Lik.' 0.0854457 (df=6)
1 - (logLik(model_numeric)/logLik(nullmod))
## 'log Lik.' 0.08249125 (df=4)

Interpretasi

Nilai ini menyerupai pseudo R². Artinya, baik model nominal maupun numeric hanya menjelaskan sekitar 8–9% variasi dalam keputusan membeli asuransi.

Ini umum pada model logistik karena variabilitas seringkali besar dan kompleks.

Visualisasi Prediksi

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(dplyr)

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"))

Model nominal

ggplot(sim_data_nominal, aes(x = trip_duration, y = predicted, color = risk_level)) +
  geom_point(alpha = 0.5) +
  labs(title = "Prediksi Probabilitas Pembelian Asuransi (Risk Level sebagai Dummy)",
       x = "Durasi Perjalanan (hari)", y = "Probabilitas Membeli") +
  theme_minimal()

Interpretasi

Grafik menunjukkan bahwa semakin lama durasi perjalanan, probabilitas membeli asuransi meningkat untuk semua level risiko.

Perbedaan antara level risiko terlihat jelas: “Very High” menghasilkan probabilitas tertinggi di semua durasi.

Model numeric

ggplot(sim_data_numeric, aes(x = trip_duration, y = predicted, color = as.factor(risk_numeric))) +
  geom_point(alpha = 0.5) +
  labs(title = "Prediksi Probabilitas Pembelian Asuransi (Risk Level sebagai Skor)",
       x = "Durasi Perjalanan (hari)", y = "Probabilitas Membeli") +
  theme_minimal()

Interpretasi

Tren yang serupa seperti pada model nominal, tapi di sini level risiko diperlakukan sebagai skala bertingkat.

Prediksi cenderung lebih “halus” dan linier karena risk level dipetakan menjadi 1–4.

Ringkasan model nominal

library(broom)
## Warning: package 'broom' was built under R version 4.3.3
library(knitr)
library(kableExtra)

tidy(model_nominal) %>%
  mutate(estimate = round(estimate, 3),
         std.error = round(std.error, 3),
         p.value = round(p.value, 3)) %>%
  kable(format = "html", caption = "Ringkasan Koefisien Model Pembelian Asuransi (Nominal)") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
Ringkasan Koefisien Model Pembelian Asuransi (Nominal)
term estimate std.error statistic p.value
(Intercept) -1.773 0.350 -5.063636 0.000
travel_typeInternational 0.775 0.191 4.051494 0.000
risk_levelMedium 0.355 0.223 1.589841 0.112
risk_levelHigh 0.820 0.276 2.973290 0.003
risk_levelVery High 1.819 0.373 4.877162 0.000
trip_duration 0.113 0.036 3.133421 0.002

Interpretasi

Sama dengan interpretasi sebelumnya: variabel perjalanan internasional, risiko tinggi, risiko sangat tinggi, dan durasi semuanya secara signifikan meningkatkan kemungkinan membeli asuransi.

Ringkasan model numeric

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

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

Pemilihan Model Regresi Logistik dan Evaluasi

Membangun Model Regresi Logistik: Pendekatan Confrmatory dan Exploratory

Simulasi Data

Data terdiri dari 300 observasi dengan tiga prediktor: tingkat pendidikan (dalam tahun), status tempat tinggal (kota atau desa), dan tingkat kepercayaan terhadap pemerintah (skala 1–10). Variabel respons adalah kepatuhan memakai masker (1 = patuh, 0 = tidak), yang dibentuk berdasarkan kombinasi logit dari ketiga prediktor. Simulasi ini digunakan untuk mengevaluasi model regresi logistik dan membandingkan kinerjanya dengan berbagai metode seleksi dan evaluasi model.

library(knitr)
library(dplyr)
library(ggplot2)
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(caret)
## Warning: package 'caret' was built under R version 4.3.3
## Loading required package: lattice
library(pROC)
## Warning: package 'pROC' was built under R version 4.3.3
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
library(DescTools)
## Warning: package 'DescTools' was built under R version 4.3.3
## 
## Attaching package: 'DescTools'
## The following objects are masked from 'package:caret':
## 
##     MAE, RMSE
set.seed(123)

n <- 300

# x1: Tingkat edukasi (rata-rata sekitar 14 tahun pendidikan formal)
x1 <- rnorm(n, mean = 14, sd = 2)

# x2: Status tinggal di kota besar (1 = kota, 0 = desa)
x2 <- rbinom(n, 1, 0.6)

# x3: Kepercayaan terhadap pemerintah (skala 1–10)
x3 <- round(runif(n, min = 1, max = 10), 1)

# Model logit: konstruksi logika dari pengaruh faktor-faktor terhadap kepatuhan masker
lin_pred <- -2 + 0.25 * x1 + 1 * x2 + 0.3 * x3

# Probabilitas menggunakan masker
p <- 1 / (1 + exp(-lin_pred))

# y: Kepatuhan (1 = patuh, 0 = tidak)
y <- rbinom(n, 1, p)

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

Pemilihan 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)  2.88296    2.55022   1.130    0.258
## x1          -0.04047    0.17207  -0.235    0.814
## x2           0.87110    0.66178   1.316    0.188
## x3           0.12069    0.13544   0.891    0.373
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 87.687  on 299  degrees of freedom
## Residual deviance: 85.267  on 296  degrees of freedom
## AIC: 93.267
## 
## Number of Fisher Scoring iterations: 6

Interpretasi

  • Koefisien tidak signifikan (p-value semua > 0.05), artinya tidak ada prediktor yang secara signifikan mempengaruhi kepatuhan masker.
  • AIC = 93.27, digunakan untuk perbandingan model.

Metode Stepwise: Forward, Backward, dan Kedua Arah

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)

Interpretasi

  • AIC model stepwise lebih rendah (89.68) dibandingkan model penuh (93.27), artinya model stepwise lebih baik.

Evaluasi Model: ROC dan AUC

pred_prob <- predict(step_both, type = "response")
roc_obj <- roc(df$y, pred_prob)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj, main = "Kurva ROC", col = "blue")

Interpretasi

  • Kurva ROC menggambarkan hubungan antara True Positive Rate (Sensitivity) di sumbu Y dan False Positive Rate (1 - Specificity) di sumbu X untuk berbagai nilai ambang (threshold).

  • Titik (0,1) pada grafik adalah titik ideal (100% sensitif dan 100% spesifik)

auc(roc_obj)
## Area under the curve: 0.5

Interpretasi

  • AUC = 0.5, artinya model tidak memiliki kemampuan klasifikasi yang lebih baik dari tebakan acak.

Pseudo R-Squared

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

Interpretasi

  • Semua nilai = 0 → model sangat buruk dalam menjelaskan variasi variabel respons.

Tabel Klasifkasi dan Evaluasi

pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
conf_matrix <- confusionMatrix(factor(pred_class), df$y, positive = "1")
## Warning in confusionMatrix.default(factor(pred_class), df$y, positive = "1"):
## Levels are not in the same order for reference and data. Refactoring data to
## match.
conf_matrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0   0   0
##          1  10 290
##                                           
##                Accuracy : 0.9667          
##                  95% CI : (0.9396, 0.9839)
##     No Information Rate : 0.9667          
##     P-Value [Acc > NIR] : 0.583052        
##                                           
##                   Kappa : 0               
##                                           
##  Mcnemar's Test P-Value : 0.004427        
##                                           
##             Sensitivity : 1.0000          
##             Specificity : 0.0000          
##          Pos Pred Value : 0.9667          
##          Neg Pred Value :    NaN          
##              Prevalence : 0.9667          
##          Detection Rate : 0.9667          
##    Detection Prevalence : 1.0000          
##       Balanced Accuracy : 0.5000          
##                                           
##        'Positive' Class : 1               
##

Interpretasi

  • Model hanya memprediksi kelas 1 (patuh).
  • Semua data kelas 1 (290) diprediksi benar → Sensitivity = 1
  • Semua data kelas 0 (10) diprediksi salah → Specificity = 0
  • Akurasi = 96.67% terlihat tinggi, tapi Kappa = 0 menunjukkan tidak ada kesesuaian selain kebetulan.
  • Balanced Accuracy: 0.5 → Seimbang hanya secara matematis
  • Mcnemar’s Test: p = 0.004 → Perbedaan signifikan antara kesalahan klasifikasi
conf_matrix$byClass[c("Sensitivity", "Specificity")]
## Sensitivity Specificity 
##           1           0

Interpretasi

  • Sensitivitas = 1, Spesifisitas = 0 → model hanya mampu mendeteksi kelas mayoritas.

Metode Perbandingan Model dalam Regresi Logistik

library(MASS)
library(broom)
library(DescTools)

Simulasi Data

Data terdiri dari 300 responden, dengan variabel respons berupa status kepemilikan asuransi (1 = memiliki, 0 = tidak). Variabel prediktor terdiri dari usia (x1), status pekerjaan (x2), dan pendapatan bulanan (x3).

set.seed(456)

n <- 300

# x1: Usia (20–60 tahun)
x1 <- round(runif(n, 20, 60), 0)

# x2: Status pekerjaan (1 = tetap, 0 = tidak tetap)
x2 <- rbinom(n, 1, 0.6)

# x3: Pendapatan bulanan (rata-rata 5 juta, sd 2 juta)
x3 <- round(rnorm(n, mean = 5, sd = 2), 2)
x3[x3 < 0] <- 0.5  # pendapatan tidak boleh negatif

# Linear predictor
lin_pred <- -3 + 0.06 * x1 + 1.2 * x2 + 0.4 * x3

# Probabilitas memiliki asuransi
p <- 1 / (1 + exp(-lin_pred))

# Respon: y = memiliki asuransi (1) atau tidak (0)
y <- rbinom(n, 1, p)

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

Pembuatan Model

# Model 1: hanya usia
model1 <- glm(y ~ x1, data = data, family = binomial)

# Model 2: usia + status pekerjaan
model2 <- glm(y ~ x1 + x2, data = data, family = binomial)

# Model 3: semua variabel
model3 <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)

Perbandingan AIC dan Deviance

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

Interpretasi

  • Penambahan variabel menurunkan AIC dan deviance.
  • Model 3 terbaik dengan AIC = 224.53.

Likelihood-Ratio Test

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

Interpretasi

  • Setiap penambahan variabel memberikan peningkatan signifikan terhadap model (p < 0.05).

Prinsip Parsimony

Evaluasi Tabel Klasifkasi 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  13   7
##          1  37 243
##                                           
##                Accuracy : 0.8533          
##                  95% CI : (0.8082, 0.8914)
##     No Information Rate : 0.8333          
##     P-Value [Acc > NIR] : 0.1984          
##                                           
##                   Kappa : 0.3053          
##                                           
##  Mcnemar's Test P-Value : 1.232e-05       
##                                           
##             Sensitivity : 0.9720          
##             Specificity : 0.2600          
##          Pos Pred Value : 0.8679          
##          Neg Pred Value : 0.6500          
##              Prevalence : 0.8333          
##          Detection Rate : 0.8100          
##    Detection Prevalence : 0.9333          
##       Balanced Accuracy : 0.6160          
##                                           
##        'Positive' Class : 1               
##

Interpretasi

  • Akurasi: 85.33% → Cukup tinggi
  • Kappa: 0.305 → Kesesuaian sedang
  • Balanced Accuracy: 0.616 → Kinerja total terhadap kedua kelas cukup baik
  • PPV (Precision): 0.8679 → Prediksi positif sebagian besar benar
  • NPV: 0.65 → Prediksi negatif masih bisa dipercaya, tapi tidak kuat
  • Mcnemar’s p-value: < 0.001 → Ada ketidakseimbangan prediksi FN dan FP
  • Kesimpulan: Model sangat baik dalam mengenali kelas 1, tapi masih lemah dalam deteksi kelas 0.

Sensitivitas dan Spesifsitas

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

Interpretasi

  • Sensitivity: 0.972 → Model sangat baik mendeteksi yang memiliki asuransi

  • Specificity: 0.260 → Lemah dalam mengenali yang tidak memiliki

Detail ROCPenjelasan Kurva ROC (Receiver Operating Characteristic)

Visualisasi dalam R

library(pROC)
set.seed(123)
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)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj)

auc(roc_obj)
## Area under the curve: 0.8686

Interpretasi

  • AUC = 0.868 → model memiliki performa klasifikasi yang sangat baik.

Simulasi Pemilihan Threshold Optimal

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.91489362   0.5947712
## 2       0.15  0.85106383   0.6862745
## 3       0.20  0.80851064   0.7320261
## 4       0.25  0.76595745   0.7712418
## 5       0.30  0.72340426   0.8104575
## 6       0.35  0.68085106   0.8366013
## 7       0.40  0.61702128   0.8954248
## 8       0.45  0.59574468   0.9150327
## 9       0.50  0.51063830   0.9281046
## 10      0.55  0.51063830   0.9477124
## 11      0.60  0.42553191   0.9607843
## 12      0.65  0.36170213   0.9738562
## 13      0.70  0.29787234   0.9803922
## 14      0.75  0.19148936   0.9869281
## 15      0.80  0.12765957   0.9869281
## 16      0.85  0.06382979   1.0000000
## 17      0.90  0.02127660   1.0000000

Interpretasi

  • Trade-off antara sensitivitas dan spesifisitas tergantung ambang batas klasifikasi.

Precision-Recall Curve (PR Curve)

Visualisasi PR Curve di R

library(PRROC)
## Warning: package 'PRROC' was built under R version 4.3.3
## Loading required package: rlang
set.seed(123)
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)

Interpretasi

  • PR Curve cocok jika data tidak seimbang (class imbalance).

Pseudo R-squared pada Regresi Logistik

Simulasi Data

Data terdiri dari 300 pelanggan dengan tiga prediktor: kecepatan internet (x1), kualitas layanan pelanggan (x2), dan harga langganan bulanan (x3). Variabel respons berupa kepuasan pelanggan (1 = puas, 0 = tidak puas) dibentuk dari kombinasi logit terhadap ketiga prediktor. Simulasi ini digunakan untuk menghitung nilai pseudo R-squared sebagai ukuran kebaikan model regresi logistik, dengan membandingkan model penuh terhadap model null.

set.seed(789)

n <- 300

# x1: Kecepatan internet (10–100 Mbps)
x1 <- round(runif(n, 10, 100), 1)

# x2: Layanan pelanggan memuaskan (1 = ya, 0 = tidak)
x2 <- rbinom(n, 1, 0.7)

# x3: Harga langganan bulanan (200–800 ribu rupiah)
x3 <- round(runif(n, 2, 8), 1)

# Linear predictor
lin_pred <- -4 + 0.05 * x1 + 1.5 * x2 - 0.8 * x3

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

# y: Kepuasan pelanggan (1 = puas, 0 = tidak puas)
y <- rbinom(n, 1, p)

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

Model Logistik dan Null Model

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

Perhitungan Manual R-squared

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

Perhitungan Otomatis dengan Package Tambahan

Menggunakan pscl

library(pscl)
## Warning: package 'pscl' was built under R version 4.3.3
## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.
pR2(model)
## fitting null model for pseudo-r2
##         llh     llhNull          G2    McFadden        r2ML        r2CU 
## -37.1896054 -59.5545730  44.7299353   0.3755374   0.1385169   0.4227115

Menggunakan rcompanion

#if (!require(rcompanion)) 
#install.packages("rcompanion")
#library(rcompanion)
#nagelkerke(model)

Menggunakan DescTools

library(DescTools)
PseudoR2(model, which = "all")
##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.3755374       0.3083721       0.1385169       0.4227115       0.1297536 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.4565636       0.2262188       0.6818220       0.2400787      82.3792107 
##             BIC          logLik         logLik0              G2 
##      97.1943406     -37.1896054     -59.5545730      44.7299353

Interpretasi

Model regresi logistik yang digunakan menunjukkan kinerja yang cukup baik. Nilai McFadden R² sebesar 0,375 dan Nagelkerke R² sebesar 0,423 menunjukkan bahwa model memiliki kemampuan penjelasan yang cukup kuat terhadap variabel dependen. Selain itu, uji likelihood ratio signifikan (p < 0,001), yang berarti model secara signifikan lebih baik dibandingkan model tanpa prediktor (null model).

Distribusi Multinomial

Distribusi multinomial adalah perluasan dari distribusi binomial untuk lebih dari dua kategori.

Jika \(X_1, X_2, ..., X_k\) menyatakan banyaknya kejadian dalam masing-masing dari \(k\) kategori, maka:

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

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

Studi Kasus

Sebuah survei dilakukan terhadap 10 orang tentang genre musik favorit mereka:

  • Pop: 5
  • Rock: 3
  • Jazz: 2

Probabilitas teori:

  • \(p_{Pop} = 0.4\)
  • \(p_{Rock} = 0.4\)
  • \(p_{Jazz} = 0.2\)

Hitung Peluang Manual di R:

n <- 10
x <- c(5, 3, 2)
p <- c(0.4, 0.4, 0.2)

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

Multinomial Logistic Regression

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

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

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

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

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 j = 1, \dots, c-1 \]

dengan:

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

Terdapat sebanyak \((c - 1)\) fungsi logit.

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

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

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

Contoh Kasus: 3 Kategori Respon

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

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

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

Relasi antar kategori selain baseline:

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

Estimasi Parameter

Estimasi dilakukan dengan metode Maximum Likelihood menggunakan algoritma iteratif (seperti Newton-Raphson).

Log-likelihood:

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

dengan: - \(\pi_{ij} = P(Y_i = j | x_i)\) - \(y_{ij} = 1\) jika \(Y_i = j\)

Contoh Kasus

Simulasi Data

library(nnet)
set.seed(123)

n <- 150
Genre <- sample(c("Pop", "Jazz", "Rock"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 30, sd = 5))
ListeningTime <- round(pmax(rnorm(n, mean = 10, sd = 3), 0))

Platform <- sapply(Genre, function(g) {
  if (g == "Pop") {
    sample(c("Spotify", "AppleMusic", "YouTube"), size = 1, prob = c(0.5, 0.3, 0.2))
  } else if (g == "Jazz") {
    sample(c("Spotify", "AppleMusic", "YouTube"), size = 1, prob = c(0.2, 0.5, 0.3))
  } else {
    sample(c("Spotify", "AppleMusic", "YouTube"), size = 1, prob = c(0.3, 0.3, 0.4))
  }
})

df <- data.frame(Platform = factor(Platform), Age, Genre = factor(Genre), ListeningTime)
df$Platform <- relevel(df$Platform, ref = "Spotify") # baseline
head(df)

Estimasi Model

model <- multinom(Platform ~ Age + ListeningTime + Genre, data = df)
## # weights:  18 (10 variable)
## initial  value 164.791843 
## iter  10 value 155.170878
## final  value 155.008122 
## converged
summary(model)
## Call:
## multinom(formula = Platform ~ Age + ListeningTime + Genre, data = df)
## 
## Coefficients:
##            (Intercept)          Age ListeningTime  GenrePop  GenreRock
## AppleMusic   1.6287355 -0.006700859   -0.03947295 -1.601187 -0.9527845
## YouTube      0.5124881  0.041165907   -0.09094170 -2.126068 -0.7515400
## 
## Std. Errors:
##            (Intercept)        Age ListeningTime  GenrePop GenreRock
## AppleMusic    1.381306 0.03900864    0.06746306 0.5346170 0.5181676
## YouTube       1.460503 0.04227322    0.07202683 0.6125704 0.5299225
## 
## Residual Deviance: 310.0162 
## AIC: 330.0162

Nilai P-Value dan Interpretasi

# Koefisien
coef(model)
##            (Intercept)          Age ListeningTime  GenrePop  GenreRock
## AppleMusic   1.6287355 -0.006700859   -0.03947295 -1.601187 -0.9527845
## YouTube      0.5124881  0.041165907   -0.09094170 -2.126068 -0.7515400
# P-value
z <- summary(model)$coefficients / summary(model)$standard.errors
p <- 2 * (1 - pnorm(abs(z)))
p
##            (Intercept)       Age ListeningTime     GenrePop  GenreRock
## AppleMusic   0.2383475 0.8636114     0.5584773 0.0027442927 0.06595088
## YouTube      0.7256645 0.3301530     0.2067298 0.0005190413 0.15613021

Interpretasi

  • Variabel Genre Pop dan Rock berpengaruh signifikan terhadap pemilihan platform AppleMusic dibanding Spotify, karena nilai p-value untuk GenrePop adalah 0.0027 dan GenreRock sebesar 0.0659.
  • Untuk platform YouTube, genre Pop juga signifikan dibanding Spotify (p-value = 0.0005).
  • Namun, variabel Age dan ListeningTime tidak signifikan terhadap pemilihan platform karena p-value > 0.05.

Artinya, preferensi genre lebih menentukan pemilihan platform musik daripada usia atau lama waktu mendengarkan.

Prediksi dan Validasi

df$Predicted <- predict(model)
table(Predicted = df$Predicted, Actual = df$Platform)
##             Actual
## Predicted    Spotify AppleMusic YouTube
##   Spotify         28         14      10
##   AppleMusic      12         29      17
##   YouTube          8         13      19

Interpretasi

Tabel klasifikasi menunjukkan seberapa akurat model dalam memprediksi platform musik. Contoh:

  • Dari pengguna yang menggunakan Spotify, 28 diprediksi benar, sisanya salah klasifikasi.
  • Untuk AppleMusic, 29 diprediksi benar dari total 60.
  • Untuk YouTube, prediksi tepat sebanyak 19 dari 46.

Hal ini menunjukkan bahwa model cukup baik memprediksi, terutama untuk kategori AppleMusic dan YouTube, meskipun masih terdapat mis-klasifikasi yang cukup banyak.

library(ggplot2)

ggplot(df, aes(x = Age, y = ListeningTime, color = Predicted)) +
  geom_point(size = 2) +
  labs(title = "Multinomial Logistic Regression Predictions",
       x = "Age", y = "Listening Time") +
  theme_minimal()

Interpretasi

Grafik menunjukkan hasil prediksi model regresi logistik multinomial terhadap preferensi platform musik berdasarkan usia dan waktu mendengarkan. Secara umum, Apple Music lebih dipilih oleh pengguna berusia muda dengan waktu mendengarkan tinggi, Spotify dominan pada usia menengah, dan YouTube cenderung dipilih oleh pengguna berusia lebih tua. Meskipun pola klasifikasi cukup terlihat, masih terdapat tumpang tindih antar kategori, yang menunjukkan bahwa prediksi belum sepenuhnya akurat hanya dengan dua variabel ini.

Regresi logistik ordinal

Regresi logistik ordinal digunakan ketika variabel respon \(Y\) bersifat ordinal (memiliki urutan), misalnya tingkat kepuasan: Rendah, Sedang, Tinggi. Model ini berbeda dengan: - Regresi logistik biner: hanya 2 kategori - Regresi logistik multinomial: kategori lebih dari 2 tetapi tidak berurutan

Konsep Cumulative Logit Model

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

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

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

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

Interpretasi Koefisien

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

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

Odds ratio:

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

Contoh Data

set.seed(2025)
n <- 250

# Variabel prediktor
tasks_per_week <- round(runif(n, 3, 15))         # Jumlah tugas per minggu
sleep_hours <- round(rnorm(n, mean = 6, sd = 1), 1)  # Jam tidur per malam

# Kombinasi linear (logit ordinal)
linpred <- -0.4 * sleep_hours + 0.3 * tasks_per_week + rnorm(n, 0, 1)

# Konversi menjadi tingkat stres ordinal
stress <- cut(linpred,
              breaks = c(-Inf, 3, 6, Inf),
              labels = c("Rendah", "Sedang", "Tinggi"),
              ordered_result = TRUE)

# Dataset
df_stress <- data.frame(stress, tasks_per_week, sleep_hours)
head(df_stress)

Estimasi Model Ordinal

model_ord <- polr(stress ~ tasks_per_week + sleep_hours, data = df_stress, Hess = TRUE)
(ctable <- coef(summary(model_ord)))
##                     Value Std. Error   t value
## tasks_per_week  0.2809024  0.1060779  2.648078
## sleep_hours    -0.2506101  0.2808765 -0.892243
## Rendah|Sedang   4.4097390  2.0633794  2.137144
## Sedang|Tinggi  84.8715804  2.0633794 41.132319

Nilai P-Value

(ctable <- coef(summary(model_ord)))
##                     Value Std. Error   t value
## tasks_per_week  0.2809024  0.1060779  2.648078
## sleep_hours    -0.2506101  0.2808765 -0.892243
## Rendah|Sedang   4.4097390  2.0633794  2.137144
## Sedang|Tinggi  84.8715804  2.0633794 41.132319
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
## tasks_per_week  0.2809024  0.1060779  2.648078  0.0081
## sleep_hours    -0.2506101  0.2808765 -0.892243  0.3723
## Rendah|Sedang   4.4097390  2.0633794  2.137144  0.0326
## Sedang|Tinggi  84.8715804  2.0633794 41.132319  0.0000

Prediksi Probabilitas

# Prediksi probabilitas stres untuk kombinasi tugas dan jam tidur
newdata <- data.frame(
  tasks_per_week = c(5, 8, 11),
  sleep_hours = c(7, 6, 5)
)
predict(model_ord, newdata = newdata, type = "probs")
##      Rendah      Sedang Tinggi
## 1 0.9915029 0.008497093      0
## 2 0.9750638 0.024936247      0
## 3 0.9290945 0.070905460      0

Goodness-of-Fit dan Proportional Odds

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

Alternatif Model Ordinal

Selain cumulative logit, terdapat beberapa model ordinal lainnya: - Adjacent-category logit - Continuation-ratio (sequential) logit

Model alternatif digunakan jika asumsi proportional odds tidak terpenuhi.

Kesimpulan

  • Regresi ordinal efektif untuk respon berurutan.
  • Model cumulative logit menginterpretasikan efek dalam bentuk log-odds kumulatif.
  • Implementasi di R dengan fungsi polr() dari package MASS.
  • Untuk validasi lebih lanjut, gunakan uji devians atau likelihood ratio test.

Asumsi Paralelisme dalam Regresi Logistik Ordinal

Model cumulative logit mengasumsikan proportional odds atau parallel lines assumption, yaitu:

\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j + \beta x, \quad j = 1, \dots, c - 1 \]

  • Hanya \(\alpha_j\) yang berbeda-beda.
  • Koefisien \(\beta\) tetap sama untuk semua fungsi logit kumulatif.

Visualisasi: Kurva logit kumulatif dari tiap kategori terhadap prediktor akan paralel (kemiringan sama, intercept berbeda).

Konsekuensi Pelanggaran Asumsi: - Efek prediktor berbeda untuk tiap batas kategori. - Model cumulative logit tidak valid. - Gunakan model alternatif: Generalized Ordinal Logistic Regression, Partial Proportional Odds Model.

Pengujian Asumsi Paralelisme:

Gunakan Likelihood Ratio Test atau Brant Test.

library(brant)
## Warning: package 'brant' was built under R version 4.3.3
newdata <- data.frame(
  tasks_per_week = c(5, 8, 11),
  sleep_hours = c(7, 6, 5)
)
predict(model_ord, newdata = newdata, type = "probs")
##      Rendah      Sedang Tinggi
## 1 0.9915029 0.008497093      0
## 2 0.9750638 0.024936247      0
## 3 0.9290945 0.070905460      0

Interpretasi

  • Model cumulative logit menunjukkan bahwa jumlah tugas per minggu adalah prediktor signifikan terhadap tingkat stres mahasiswa.

  • Jam tidur memiliki efek negatif tetapi tidak signifikan secara statistik.

  • Prediksi menunjukkan kecenderungan stres meningkat seiring bertambahnya tugas dan berkurangnya waktu tidur, walaupun seluruh prediksi masih dominan di tingkat “Rendah”.

Log Linear Model

Ringkasan Dalam analisis data kategorik, terdapat beberapa pendekatan statistik yang umum digunakan, antara lain:

  1. Tabel Kontingensi: penyajian frekuensi gabungan dari dua atau lebih variabel kategorik.
  2. Model Loglinear: digunakan untuk memodelkan struktur asosiasi di dalam tabel kontingensi tanpa menganggap ada variabel dependen.
  3. Model Regresi Logistik: digunakan untuk memodelkan probabilitas dari kategori variabel dependen berdasarkan variabel independen.

Meskipun ketiganya dapat digunakan pada data kategorik, pendekatan dan interpretasinya sangat berbeda.

Tabel Kontingensi

Contoh tabel 2x2

tabel_pelatihan <- matrix(c(45, 15, 30, 60), nrow=2, byrow=TRUE,
       dimnames = list(Pelatihan = c("Intensif", "Biasa"),
                       Lulus = c("Ya", "Tidak")))
tabel_pelatihan
##           Lulus
## Pelatihan  Ya Tidak
##   Intensif 45    15
##   Biasa    30    60

Model Loglinear

Model log-linier dua arah (\(I \times J\)) ditulis sebagai:

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

library(MASS)
loglm(~ Pelatihan * Lulus, data = tabel_pelatihan)
## Call:
## loglm(formula = ~Pelatihan * Lulus, data = tabel_pelatihan)
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1

Model Regresi Logistik

Model regresi logistik memodelkan:

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

data_glm <- data.frame(
  Lulus = c(1, 0, 1, 0),
  Pelatihan = factor(c("Intensif", "Intensif", "Biasa", "Biasa")),
  Frek = c(45, 15, 30, 60)
)

model_logit <- glm(Lulus ~ Pelatihan, weights = Frek, family = binomial, data = data_glm)
summary(model_logit)
## 
## Call:
## glm(formula = Lulus ~ Pelatihan, family = binomial, data = data_glm, 
##     weights = Frek)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -0.6931     0.2236  -3.100  0.00194 ** 
## PelatihanIntensif   1.7918     0.3727   4.808 1.52e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 207.94  on 3  degrees of freedom
## Residual deviance: 182.05  on 2  degrees of freedom
## AIC: 186.05
## 
## Number of Fisher Scoring iterations: 4

Interpretasi

Model Saturated

library(MASS)

# Buat tabel frekuensi (kontingensi)
tabel_pelatihan <- matrix(c(45, 15, 30, 60), nrow=2, byrow=TRUE,
       dimnames = list(Pelatihan = c("Intensif", "Biasa"),
                       Lulus = c("Ya", "Tidak")))

# Fit model saturated (interaksi penuh)
model_saturated <- loglm(~ Pelatihan * Lulus, data = tabel_pelatihan)
summary(model_saturated)
## Formula:
## ~Pelatihan * Lulus
## attr(,"variables")
## list(Pelatihan, Lulus)
## attr(,"factors")
##           Pelatihan Lulus Pelatihan:Lulus
## Pelatihan         1     0               1
## Lulus             0     1               1
## attr(,"term.labels")
## [1] "Pelatihan"       "Lulus"           "Pelatihan:Lulus"
## 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

Interpretasi

Model saturated selalu cocok sempurna dengan data (pasti p = 1), jadi tidak ada deviasi, digunakan sebagai acuan untuk membandingkan model yang lebih sederhana.

Model Independent

Model independen mengasumsikan tidak ada interaksi antar variabel:

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

# Contoh data frame
data <- data.frame(
  Pelatihan = c("Intensif", "Intensif", "Biasa", "Biasa"),
  Lulus = c("Ya", "Tidak", "Ya", "Tidak"),
  Frek = c(45, 15, 30, 60)
)

tbl <- xtabs(Frek ~ Pelatihan + Lulus, data = data)

library(MASS)
model_indep <- loglm(~ Pelatihan + Lulus, data = tbl)
summary(model_indep)
## Formula:
## ~Pelatihan + Lulus
## attr(,"variables")
## list(Pelatihan, Lulus)
## attr(,"factors")
##           Pelatihan Lulus
## Pelatihan         1     0
## Lulus             0     1
## attr(,"term.labels")
## [1] "Pelatihan" "Lulus"    
## 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 25.89139  1 3.611772e-07
## Pearson          25.00000  1 5.733031e-07

Interpretasi

Model ini tidak cocok dengan data. Nilai p sangat kecil → asumsi independensi ditolak, artinya ada interaksi nyata antara jenis pelatihan dan kelulusan.

Odds Ratio dan Interpretasi

\[ \text{OR} = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}}, \quad \log(\text{OR}) = \log\left(\frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}}\right) \]

  • OR > 1: Pelatihan intensif berhubungan positif dengan kelulusan
  • OR < 1: hubungan negatif
  • OR = 1: tidak ada asosiasi

Estimasi Parameter

tabel <- matrix(c(45, 15, 30, 60), nrow = 2, byrow = TRUE,
                dimnames = list(Pelatihan = c("Intensif", "Biasa"),
                                Lulus = c("Ya", "Tidak")))

logOR <- log((tabel[1,1] * tabel[2,2]) / (tabel[1,2] * tabel[2,1]))
logOR
## [1] 1.791759

Interpretasi

Peserta pelatihan intensif 6 kali lebih besar kemungkinannya untuk lulus dibandingkan peserta pelatihan biasa. Ada asosiasi positif yang kuat.

Model Lebih Sederhana dan Perbandingan Model

anova(model_indep, model_saturated)
## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Pelatihan + Lulus 
## Model 2:
##  ~Pelatihan * Lulus 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   25.89139  1                                    
## Model 2    0.00000  0   25.89139         1              0
## Saturated  0.00000  0    0.00000         0              1

Interpretasi

Model saturated jauh lebih baik dari model independent (ΔX² = 25.89, p = 0), interaksi harus dipertimbangkan.

Studi Kasus: Kepuasan dan Sertifikasi

data_studi <- matrix(c(20, 80,
                       50, 50,
                       40, 60),
                     nrow = 3, byrow = TRUE,
                     dimnames = list(Kepuasan = c("Rendah", "Sedang", "Tinggi"),
                                     Sertifikasi = c("Tidak", "Ya")))
ftable(data_studi)
##          Sertifikasi Tidak Ya
## Kepuasan                     
## Rendah                  20 80
## Sedang                  50 50
## Tinggi                  40 60
loglm(~ Kepuasan + Sertifikasi, data = data_studi)
## Call:
## loglm(formula = ~Kepuasan + Sertifikasi, data = data_studi)
## 
## Statistics:
##                       X^2 df     P(> X^2)
## Likelihood Ratio 20.98240  2 2.777980e-05
## Pearson          20.09569  2 4.327883e-05

Interpretasi

  • Model tanpa interaksi tidak cukup baik, berarti hubungan antara kepuasan dan sertifikasi mungkin bersifat interaktif atau kompleks.

  • Perlu model dengan interaksi untuk menjelaskan data ini dengan baik.

TUGAS MODEL LOG LINEAR 2 ARAH

Model log-linear pada tabel kontingensi.

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

\[ \log(\mu_{ij}) = \lambda + \lambda_i^{(A)} + \lambda_j^{(B)} + \lambda_{ij}^{(AB)} \] ## Perbedaan utama antara model log-linear dan model 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).

Analisis Data Tabel Kontingensi 2x2

Misalkan data :

Data:

Sakit Sehat Total
Olahraga: Ya 25 35 60
Olahraga: Tidak 15 25 40
Total 40 60 100

Notasi: - n11 = 25 - n12 = 35 - n21 = 15 - n22 = 25

Model Log Linear

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

Model tanpa interaksi (independen):

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

Estimasi Parameter

Rata-rata log frekuensi sel

\[ \lambda = \frac{1}{4} (\log(25) + \log(35) + \log(15) + \log(25)) = 3.175287 \]

Efek Utama A (Olahraga)

\[ \lambda^{A}_{1} = \frac{1}{2}[(\log(25) + \log(35)) - (\log(15) + \log(25))] / 2 = 0.4236489 \]

\[ \lambda^{A}_{2} = -0.4236489 \]

Efek Utama B (Kesehatan)

\[ \lambda^{B}_{1} = \frac{1}{2}[(\log(25) + \log(35)) - (\log(15) + \log(25))] / 2 = -0.4236489 \]

\[ \lambda^{B}_{2} = 0.4236489 \]

Efek Interaksi

\[ \lambda^{AB}_{11} = \frac{1}{4}[\log(25) - \log(35) - \log(15) + \log(25)] = 0.04358835 \]

\[ \lambda^{AB}_{12} = -0.04358835, \quad \lambda^{AB}_{21} = -0.04358835, \quad \lambda^{AB}_{22} = 0.04358835 \]

Ringkasan Parameter

  • \(\lambda = 3.175287\)
  • \(\lambda^A = (0.4236489, - 0.4236489)\)
  • \(\lambda^B = (-0.4236489, 0.4236489)\)
  • \(\lambda^{AB} = \begin{bmatrix} 0.04358835 & -0.04358835 \\ -0.04358835 & 0.04358835 \end{bmatrix}\)

Odds Ratio dan Confidence Interval

Odds Ratio (OR)

\[ OR = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} = \frac{25 \times 25}{35 \times 15} = 1.190476 \]

Log Odds Ratio

\[ \log(OR) = \log(1.190476) = 0.1743534 \]

Standard Error (SE)

\[ SE = \sqrt{\frac{1}{25} + \frac{1}{35} + \frac{1}{15} + \frac{1}{25}} = \sqrt{0.2083} = 0.4186145 \]

Interval Kepercayaan 95% untuk log(OR)

\[ 0.1743534 \pm 1.96 \times 0.4186145 = (-0.6461310,0.9948378) \]

Back-transform ke OR

\[ \text{Lower CI} = \exp(0.8968) = 0.5240695, \quad \text{Upper CI} = \exp(2.6868) = 2.7042857 \]

Kesimpulan: Odds Ratio (OR) = 1.190476, dengan 95% CI = (0.5240695, 2.7042857)

Fitting Model Log-Linear dengan R

tabel <- matrix(c(25, 35, 15, 25), nrow=2, byrow=TRUE)
dimnames(tabel) <- list(
  Olahraga = c("Ya", "Tidak"),
  Kesehatan = c("Sakit", "Sehat")
)
tabel
##         Kesehatan
## Olahraga Sakit Sehat
##    Ya       25    35
##    Tidak    15    25
data <- as.data.frame(as.table(tabel))
colnames(data) <- c("Olahraga", "Kesehatan", "Freq")
data
# Model tanpa interaksi
fit_no_inter <- glm(Freq ~ Olahraga + Kesehatan, family = poisson, data = data)
summary(fit_no_inter)
## 
## Call:
## glm(formula = Freq ~ Olahraga + Kesehatan, family = poisson, 
##     data = data)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      3.1781     0.1780  17.859   <2e-16 ***
## OlahragaTidak   -0.4055     0.2041  -1.986    0.047 *  
## KesehatanSehat   0.4055     0.2041   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: 8.22829  on 3  degrees of freedom
## Residual deviance: 0.17408  on 1  degrees of freedom
## AIC: 26.256
## 
## Number of Fisher Scoring iterations: 3
# Model dengan interaksi
fit_inter <- glm(Freq ~ Olahraga * Kesehatan, family = poisson, data = data)
summary(fit_inter)
## 
## Call:
## glm(formula = Freq ~ Olahraga * Kesehatan, family = poisson, 
##     data = data)
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    3.2189     0.2000  16.094   <2e-16 ***
## OlahragaTidak                 -0.5108     0.3266  -1.564    0.118    
## KesehatanSehat                 0.3365     0.2619   1.285    0.199    
## OlahragaTidak:KesehatanSehat   0.1744     0.4186   0.417    0.677    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 8.2283e+00  on 3  degrees of freedom
## Residual deviance: 6.6613e-15  on 0  degrees of freedom
## AIC: 28.082
## 
## Number of Fisher Scoring iterations: 3

Interpretasi

  • Model tanpa interaksi sudah cukup baik.
  • Efek Olahraga dan Kesehatan signifikan secara terpisah, tapi interaksinya tidak signifikan.
  • Artinya, pengaruh olahraga terhadap kesehatan bersifat additif, bukan tergantung satu sama lain.

Analisis Data Tabel Kontingensi 2x3

Sebuah survei dilakukan untuk meneliti hubungan antara Tingkat Pendidikan (SMA, Diploma, Sarjana) dan Status Pekerjaan (Bekerja, Tidak Bekerja). Data dari hasil pengamatan disajikan dalam tabel berikut:

Bekerja Tidak Bekerja
SMA 25 15
Diploma 30 10
Sarjana 28 8

Bentuk umum model log-linear untuk tabel 3x2 (dengan sum-to-zero constraint) adalah sebagai berikut:

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

dengan:

  • \(\mu_{ij}\): ekspektasi frekuensi pada baris ke-\(i\), kolom ke-\(j\)
  • A: Tingkat Pendidikan
    (i = 1: SMA, i = 2: Diploma, i = 3: Sarjana)
  • B: Status Pekerjaan
    (j = 1: Bekerja, j = 2: Tidak Bekerja)
  • Constraint: \[ \sum_i \lambda_i^A = 0,\quad \sum_j \lambda_j^B = 0,\quad \sum_i \lambda_{ij}^{AB} = 0,\quad \sum_j \lambda_{ij}^{AB} = 0 \]

Secara eksplisit:

\[ \log(\mu_{ij}) = \lambda + \begin{cases} \lambda_1^A & \text{(SMA)} \\ \lambda_2^A & \text{(Diploma)} \\ \lambda_3^A & \text{(Sarjana)} \end{cases} + \begin{cases} \lambda_1^B & \text{(Bekerja)} \\ \lambda_2^B & \text{(Tidak Bekerja)} \end{cases} + \lambda_{ij}^{AB} \]

Fitting Model Log-Linear di R

# Membuat data frame dari tabel 3x2
tabel3x2 <- matrix(c(25, 15, 30, 10, 28, 8), nrow = 3, byrow = TRUE)
colnames(tabel3x2) <- c("Bekerja", "Tidak Bekerja")
rownames(tabel3x2) <- c("SMA", "Diploma", "Sarjana")
tabel3x2
##         Bekerja Tidak Bekerja
## SMA          25            15
## Diploma      30            10
## Sarjana      28             8
# Ubah ke format data.frame untuk analisis GLM
data3x2 <- as.data.frame(as.table(tabel3x2))
colnames(data3x2) <- c("Pendidikan", "StatusKerja", "Freq")
data3x2
# Model log-linear tanpa interaksi (asumsi independen)
fit_no_inter <- glm(Freq ~ Pendidikan + StatusKerja, family = poisson, data = data3x2)
summary(fit_no_inter)
## 
## Call:
## glm(formula = Freq ~ Pendidikan + StatusKerja, family = poisson, 
##     data = data3x2)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               3.354e+00  1.686e-01  19.893  < 2e-16 ***
## PendidikanDiploma        -4.880e-16  2.236e-01   0.000    1.000    
## PendidikanSarjana        -1.054e-01  2.297e-01  -0.459    0.647    
## StatusKerjaTidak Bekerja -9.223e-01  2.058e-01  -4.482  7.4e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 25.0386  on 5  degrees of freedom
## Residual deviance:  2.4852  on 2  degrees of freedom
## AIC: 38.622
## 
## Number of Fisher Scoring iterations: 4
# Model log-linear dengan interaksi (untuk cek asosiasi)
fit_inter <- glm(Freq ~ Pendidikan * StatusKerja, family = poisson, data = data3x2)
summary(fit_inter)
## 
## Call:
## glm(formula = Freq ~ Pendidikan * StatusKerja, family = poisson, 
##     data = data3x2)
## 
## Coefficients:
##                                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                  3.2189     0.2000  16.094   <2e-16
## PendidikanDiploma                            0.1823     0.2708   0.673    0.501
## PendidikanSarjana                            0.1133     0.2752   0.412    0.680
## StatusKerjaTidak Bekerja                    -0.5108     0.3266  -1.564    0.118
## PendidikanDiploma:StatusKerjaTidak Bekerja  -0.5878     0.4899  -1.200    0.230
## PendidikanSarjana:StatusKerjaTidak Bekerja  -0.7419     0.5171  -1.435    0.151
##                                               
## (Intercept)                                ***
## PendidikanDiploma                             
## PendidikanSarjana                             
## StatusKerjaTidak Bekerja                      
## PendidikanDiploma:StatusKerjaTidak Bekerja    
## PendidikanSarjana:StatusKerjaTidak Bekerja    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.5039e+01  on 5  degrees of freedom
## Residual deviance: 3.7748e-15  on 0  degrees of freedom
## AIC: 40.136
## 
## Number of Fisher Scoring iterations: 3

Interpretasi

  • Model tanpa interaksi (Pendidikan + StatusKerja) sudah cukup baik → tidak perlu interaksi.

  • Tidak terdapat asosiasi signifikan antara tingkat pendidikan dan status pekerjaan berdasarkan data ini.

  • Satu-satunya efek signifikan adalah: orang cenderung lebih banyak bekerja daripada tidak bekerja, terlepas dari tingkat pendidikannya.

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.

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:

Model Saturated

Model ini memuat semua kemungkinan interaksi, termasuk interaksi tiga arah (X, Y, dan Z):

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

Model ini hanya mengakomodasi interaksi dua arah antar variabel tanpa memasukkan interaksi tiga arah:

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

Model Conditional

a. Conditional pada X Memuat interaksi X dengan Y dan X dengan Z:

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

b. Conditional pada Y Memuat interaksi Y dengan X dan Y dengan Z:

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

c. Conditional pada Z Memuat interaksi Z dengan X dan Z dengan Y:

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

Model Joint Independence

a. Independensi antara X dan Y:

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

b. Independensi antara X dan Z:

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

c. Independensi antara Y dan Z:

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

5. Model Tanpa Interaksi

Model ini hanya memasukkan efek utama tanpa interaksi antar variabel:

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

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 homogen.
  • Jika perbedaan signifikan → terdapat interaksi tiga arah.

** 2. Pengujian Interaksi Dua Arah (XY, XZ, YZ)**

  • Bandingkan model homogen 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.

Studi Kasus

Tabel berikut menyajikan data hasil survei nasional yang dilakukan pada tahun 2022 mengenai jenis pekerjaan, jenis kelamin responden, dan preferensi kerja jarak jauh (remote work). 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.

Tabel Data Survei

Pekerjaan Jenis Kelamin Setuju Remote Tidak Setuju
Profesional Laki-laki 140 40
Profesional Perempuan 130 50
Administratif Laki-laki 100 70
Administratif Perempuan 120 60
Pekerja Lapangan Laki-laki 60 100
Pekerja Lapangan Perempuan 70 90

Keterangan:

  • Pekerjaan: Profesional, Administratif, Pekerja Lapangan
  • Jenis Kelamin: Laki-laki, Perempuan
  • Preferensi Kerja Jarak Jauh: Setuju (Favor), Tidak Setuju (Oppose)

Analisis Log-Linear untuk Tabel Tiga Arah

Package yang Digunakan

library(epitools)
library(DescTools)
library(lawstat)
## Warning: package 'lawstat' was built under R version 4.3.3

Input Data

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

Membentuk Tabel Kontingensi 3 Arah

table3d <- xtabs(Frekuensi ~ Fundamentalisme + Jenis_Kelamin + Sikap, data = data)
ftable(table3d)
##                               Sikap 1fav 2opp
## Fundamentalisme Jenis_Kelamin                
## 1fund           1M                   128   32
##                 2F                   123   73
## 2mod            1M                   182   56
##                 2F                   168  105
## 3lib            1M                   119   49
##                 2F                   111   70

UJI MODEL INTERAKSI TIGA ARAH (SATURATED VS HOMOGENOUS)

Penentuan Kategori Referensi

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

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

# 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(coef(model_saturated))
##                   (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

Hasil Estimasi Koefisien

Hasil Estimasi Koefisien

Berdasarkan output model saturated yang telah diperoleh:

Term Estimate Std. Error z value Pr(> z
(Intercept) 4.49981 0.10541 42.689 < 2e-16 90.00
PekerjaanAdministratif -0.40547 0.16667 -2.433 0.015 0.67
PekerjaanProfesional -0.58779 0.17638 -3.332 0.0009 0.56
Jenis_KelaminLaki-laki 0.10536 0.14530 0.725 0.468 1.11
PreferensiSetuju -0.25131 0.15936 -1.577 0.115 0.78
PekerjaanAdministratif:Jenis_KelaminLaki-laki 0.04879 0.22817 0.214 0.831 1.05
PekerjaanProfesional:Jenis_KelaminLaki-laki -0.32850 0.25712 -1.278 0.201 0.72
PekerjaanAdministratif:PreferensiSetuju 0.94446 0.22449 4.207 2.59e-05 2.57
PekerjaanProfesional:PreferensiSetuju 1.20683 0.23041 5.238 1.63e-07 3.34
Jenis_KelaminLaki-laki:PreferensiSetuju -0.25951 0.22817 -1.137 0.255 0.77
PekerjaanAdministratif:Jenis_KelaminLaki-laki:PreferensiSetuju -0.07696 0.31835 -0.242 0.809 0.93
PekerjaanProfesional:Jenis_KelaminLaki-laki:PreferensiSetuju 0.55676 0.33451 1.664 0.096 1.75

Interpretasi Koefisien

  • Intercept : Frekuensi dasar responden pada kategori referensi (Pekerja Lapangan, Perempuan, Tidak Setuju) sekitar 90 orang.
  • Pekerjaan : Responden Administratif dan Profesional lebih sedikit secara umum dibanding Pekerja Lapangan, tapi lebih cenderung setuju kerja jarak jauh (interaksi signifikan).
  • Jenis Kelamin : Tidak ada perbedaan signifikan antara laki-laki dan perempuan dalam jumlah responden.
  • Preferensi : Tidak ada perbedaan signifikan frekuensi antara setuju dan tidak setuju pada kategori dasar.
  • Interaksi Pekerjaan-Preferensi : Responden Administratif dan Profesional lebih sering setuju kerja remote dibanding Pekerja Lapangan.
  • Interaksi Jenis Kelamin-Preferensi : Tidak signifikan, jenis kelamin tidak berpengaruh pada preferensi.
  • Interaksi Tiga Arah : Tidak signifikan, model dua arah sudah cukup menjelaskan hubungan antar variabel.

Goodness-of-Fit

  • Nilai deviance residual model saturated adalah 0, menunjukkan model fit dengan sempurna data.

  • ANOVA membandingkan model saturated dengan model homogeneous menunjukkan perbedaan yang tidak signifikan (p > 0.05), mendukung bahwa model dua arah sudah cukup.

Kesimpulan

  • Asosiasi utama dalam data ini adalah hubungan antara jenis pekerjaan dan preferensi kerja remote.
  • Jenis kelamin tidak berperan signifikan dalam preferensi kerja jarak jauh.
  • Model log-linear dua arah memberikan gambaran parsimonious yang baik tanpa kehilangan informasi penting.
  • Prediksi frekuensi dari model menunjukkan bahwa pekerja profesional dan administratif lebih cenderung setuju kerja jarak jauh dibanding pekerja lapangan.

Model Homogenous

Model log-linear homogenus 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

Uji Hipotesis

Hipotesis

• H0: Tidak ada interaksi tiga arah (model homogenous sudah cukup) • H1: Ada interaksi tiga arah (model saturated diperlukan)

Hitung Selisih Deviance

Deviance.model <- model_homogenous$deviance - model_saturated$deviance
Deviance.model
## [1] 1.797977

Hitung Derajat Bebas

derajat.bebas <- model_homogenous$df.residual - model_saturated$df.residual
derajat.bebas
## [1] 2

Hitung Chi-Square Tabel

chi.tabel <- qchisq(1 - 0.05, df = derajat.bebas)
chi.tabel
## [1] 5.991465

Keputusan Uji

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

Rangkuman

  • Hipotesis
    • \(H_0: \lambda^{XYZ} = 0\) (Tidak ada interaksi tiga arah; model yang terbentuk adalah model homogenous)
    • \(H_1: \lambda^{XYZ} \ne 0\) (Ada interaksi tiga arah; model yang terbentuk adalah model saturated)
  • Tingkat Signifikansi
    • \(\alpha = 5\%\)
  • Statistik Uji
    • \(\Delta \text{Deviance} = \text{Deviance model homogenous} - \text{Deviance model saturated}\)
    • \(\Delta \text{Deviance} = 4.254 - 0 = 4.254\)
    • \(df = df_{homogenous} - df_{saturated} = 2 - 0 = 2\)
  • Kriteria Uji
    • Tolak \(H_0\) jika \(\Delta \text{Deviance} > \chi^2_{(0.05, 2)} = 5.991\)
  • Keputusan
    • Karena \(4.254 < 5.991\), maka tidak tolak \(H_0\)
  • Interpretasi
    Pada taraf nyata 5%, belum cukup bukti untuk menolak \(H_0\) sehingga tidak ada interaksi tiga arah antara jenis kelamin, pekerjaan, dan preferensi kerja jarak jauh.

Artinya, model homogenous (yang hanya memasukkan interaksi dua arah) sudah cukup menjelaskan hubungan dalam data.

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

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

\(H_0 : \chi^2_{jk} = 0\) (Tidak ada interaksi antara pendapat hukuman mati (Y) dan fundamentalisme (Z))
\(H_1 : \chi^2_{jk} \ne 0\) (Ada interaksi antara pendapat hukuman mati (Y) dan fundamentalisme (Z))

  • Tingkat Signifikansi

\(\alpha = 5\%\)

  • Statistik Uji

  • \(\Delta\)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 \(H_0\) jika \(\Delta\)Deviance \(> \chi^2_{0.05, 2} = 5.991\)

# Pengujian Hipotesis
# Deviance of Model
Deviance.model <- model_conditional_X$deviance - model_homogenous$deviance 
Deviance.model
## [1] 2.132302
  • Keputusan

Karena \(2.132 < 5.991\), maka tidak tolak \(H_0\)

  • Kesimpulan

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

Pengujian Selisih Deviance (Conditional on X vs Homogenous)

Uji Hipotesis menggunakan residual deviance

# Selisih deviance antar model
Deviance.model <- model_conditional_X$deviance - model_homogenous$deviance
Deviance.model
## [1] 2.132302

HItung Derajat Bebas

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

Nilai Chi-Square Tabel

chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 5.991465

Keputusan Uji

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

Interpretasi

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 sampaidua interaksi dengan X (conditional on X), sehingga interaksi Y*Z tidak signifikan secara statistik.

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

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

  • \(H_0 : \lambda^{XZ}_{ik} = 0\) (Tidak ada interaksi antara jenis kelamin (X) dan fundamentalisme (Z))

  • \(H_1 : \lambda^{XZ}_{ik} \ne 0\) (Ada interaksi antara jenis kelamin (X) dan fundamentalisme (Z))

  • Tingkat Signifikansi

\(\alpha = 5\%\)

  • Statistik Uji

  • \(\Delta\)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 \(H_0\) jika \(\Delta\)Deviance \(> \chi^2_{0.05,2} = 5.991\)

# Deviance of Model
Deviance.model <- model_conditional_Y$deviance - model_homogenous$deviance # model_conditional_Y: conditional Deviance.model
  • Keputusan

Karena \(1.1223 < 5.991\), maka tidak tolak \(H_0\)

  • Kesimpulan

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

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

# Selisih deviance antar model
Deviance.model <- model_conditional_X$deviance - model_homogenous$deviance
Deviance.model
## [1] 2.132302

HItung Derajat Bebas

derajat.bebas <- (4 - 2)
derajat.bebas
## [1] 2

Nilai Chi-Square Tabel

chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 5.991465

Keputusan Uji

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

Interpretasi dan Kesimpulan

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

Pada taraf nyata 5%, belum cukup bukti untuk menolak \(H_0\).
Artinya, tidak ada interaksi antara jenis kelamin (X) dan fundamentalisme (Z) yang signifikan secara statistik.
Model tanpa parameter \(\lambda^{XZ}_{ik}\) sudah cukup baik untuk data ini.

UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON Z)Uji Model Interaksi 2 arah (homogenous vs conditional on Z)

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

  • \(H_0 : \lambda^{XY}_{ij} = 0\) (Tidak ada interaksi antara jenis kelamin (X) dan pendapat tentang hukuman mati (Y))

  • \(H_1 : \lambda^{XY}_{ij} \ne 0\) (Ada interaksi antara jenis kelamin (X) dan pendapat tentang hukuman mati (Y))

  • Tingkat Signifikansi

\(\alpha = 5\%\)

  • Statistik Uji

  • \(\Delta\)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 \(H_0\) jika \(\Delta\)Deviance \(> \chi^2_{0.05,1} = 3.841\)

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

Karena \(27.931 > 3.841\), maka tolak \(H_0\)

  • 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 \(\lambda^{XY}_{ij}\).

Pengujian Hipotesis Interaksi X dan Z (Conditional on Y 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

HItung Derajat Bebas

derajat.bebas <- (3 - 2)
derajat.bebas
## [1] 1

Nilai Chi-Square Tabel

chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 3.841459

Keputusan Uji

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

Interpretasi dan Kesimpulan

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

Pada taraf nyata 5%, terdapat bukti yang cukup untuk menolak \(H_0\).
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 \(\lambda^{XY}_{ij}\).

PEMILIHAN MODEL TERBAIK

Ringkasan Model Log Linier

Ringkasan Model Log Linier
Model Parameter Deviance Jumlah.Parameter df AIC
Saturated \(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk}\) 0.0000 12 0 100.140
Homogenous \(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk}\) 1.7980 10 2 97.934
Conditional on X \(\lambda + \lambda^Y_j + \lambda^Z_k + \lambda^{YZ}_{jk}\) 3.9303 8 4 96.067
Conditional on Y \(\lambda + \lambda^X_i + \lambda^Z_k + \lambda^{XZ}_{ik}\) 2.9203 8 4 95.057
Conditional on Z \(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^{XY}_{ij}\) 29.7290 9 3 123.870

Ringkasan Pengujian Interaksi 3 Arah dan 2 Arah

Ringkasan Pengujian Interaksi
Interaksi Pengujian Delta.deviance Delta.df Chi.square.Tabel Keputusan Keterangan
XYZ Saturated vs Homogenous 1.7980 2 5.991 Tidak Tolak H₀ tidak ada interaksi
YZ Conditional on X vs Homogenous 2.1323 2 5.991 Tidak Tolak H₀ tidak ada interaksi
XZ Conditional on Y vs Homogenous 1.1223 2 5.991 Tidak Tolak H₀ tidak ada interaksi
XY Conditional on Z vs Homogenous 27.9310 1 3.841 Tolak H₀ ada interaksi

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.

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

INTERPRETASI KOEFISIEN MODEL TERBAIK

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

Interpretasi

  • \(\exp(\lambda^X_{1\text{(Male)}}) = \exp(-0.593) = 0.552 \rightarrow\) 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(\lambda^Y_{\text{favor}}) = \exp(0.483) = 1.621 \rightarrow\) nilai odds

Tanpa memperhatikan jenis kelamin dan fundamentalisme, peluang seseorang mendukung hukuman mati adalah 1.621 kali dibandingkan yang menolak.

  • \(\exp(\lambda^Z_{\text{fund}}) = \exp(0.01986) = 1.02 \rightarrow\) nilai odds

Tanpa memperhatikan jenis kelamin dan pendapat mengenai hukuman mati, peluang seseorang fundamentalis adalah 1.02 kali dibandingkan liberal.

  • \(\exp(\lambda^Z_{\text{mod}}) = \exp(0.381) = 1.464 \rightarrow\) nilai odds

Tanpa memperhatikan jenis kelamin dan pendapat mengenai hukuman mati, peluang seseorang moderate adalah 1.464 kali dibandingkan liberal.

  • \(\exp(\lambda^{XY}_{1\text{(M)},\ \text{favor}}) = \exp(0.658) = 1.932 \rightarrow\) 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.

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
)

Referensi