Kata Pengantar

Puji dan syukur penulis panjatkan ke hadirat Tuhan Yang Maha Esa, karena atas rahmat dan karunia-Nya, penulis dapat menyelesaikan penulisan eBook ini yang berjudul “Analisis Data Kategori”.
eBook ini disusun sebagai hasil dari proses pembelajaran dan pendalaman materi pada mata kuliah Analisis Data Kategori yang diasuh oleh Dr. I Gede Nyoman Mindra Jaya, dosen Program Studi Statistika di Universitas Padjadjaran. Selama mengikuti perkuliahan, penulis memperoleh banyak wawasan penting mengenai prinsip-prinsip dasar hingga lanjutan dalam analisis data kategorik, meliputi teori probabilitas diskrit, model log-linear, regresi logistik biner, multinomial, ordinal, serta berbagai uji asosiasi pada data dalam bentuk tabel kontingensi.

Tujuan dari penyusunan eBook ini adalah untuk merangkum dan menyusun kembali materi yang telah dipelajari ke dalam format yang lebih sistematis, aplikatif, dan mudah dipahami oleh mahasiswa maupun peneliti yang berminat dalam analisis data kategorik. Beberapa pendekatan analisis dijelaskan secara lengkap, mulai dari teori dasar, rumus matematis, ilustrasi kasus nyata, hingga interpretasi hasil menggunakan perangkat lunak statistik.

Penulis menyadari bahwa eBook ini masih memiliki berbagai kekurangan, baik dalam isi maupun penyajiannya. Oleh karena itu, penulis sangat mengharapkan kritik dan saran yang membangun demi penyempurnaan pada edisi-edisi berikutnya.

Akhir kata, penulis mengucapkan terima kasih yang tulus kepada Dr. I Gede Nyoman Mindra Jaya atas ilmu, bimbingan, serta inspirasi akademik yang telah diberikan. Semoga eBook ini dapat memberikan manfaat bagi siapa pun yang ingin memperdalam pemahaman tentang analisis data kategori, baik dalam konteks akademik maupun penelitian terapan.

Pendahuluan

Tujuan Analisis Data Kategori

Analisis data kategori merupakan pendekatan penting dalam statistik dan ilmu data yang berfokus pada variabel yang nilainya berbentuk kategori, bukan angka kontinu. Tujuan utama dari analisis data kategori meliputi beberapa aspek berikut:

Mengidentifikasi Pola dan Tren

Salah satu tujuan utama analisis data kategori adalah untuk menemukan pola atau kecenderungan dalam distribusi kategori. Misalnya, dengan mengetahui proporsi jenis kelamin dalam suatu populasi, atau kecenderungan preferensi konsumen terhadap suatu produk berdasarkan umur atau wilayah.

Menganalisis Hubungan Antarvariabel

Analisis data kategori juga digunakan untuk mengevaluasi apakah terdapat hubungan antara dua atau lebih variabel kategori. Contohnya, apakah ada hubungan antara tingkat pendidikan dan preferensi politik seseorang. Teknik umum yang digunakan meliputi tabel kontingensi dan uji chi-square.

Membantu dalam Pengambilan Keputusan

Data kategori yang dianalisis secara tepat dapat menjadi dasar kuat dalam proses pengambilan keputusan, terutama dalam bidang bisnis, pemerintahan, maupun penelitian sosial. Misalnya, menentukan strategi pemasaran berdasarkan segmentasi pelanggan berdasarkan jenis kelamin dan status pernikahan.

Mengembangkan Model Prediktif

Model prediktif seperti regresi logistik dikembangkan untuk memprediksi probabilitas kejadian berdasarkan variabel kategori. Contoh aplikasinya adalah memprediksi kemungkinan pasien menderita suatu penyakit berdasarkan jenis kelamin, status merokok, dan riwayat keluarga.

Definisi dan Ruang Lingkup Analisis Data Kategori

Data kategori adalah data yang mengelompokkan objek ke dalam kategori yang berbeda, yang tidak diukur secara numerik tetapi diklasifikasikan.

Nominal vs Ordinal

  • Data Nominal adalah data kategori yang tidak memiliki urutan, seperti warna, jenis kelamin, atau agama.
  • Data Ordinal adalah data kategori yang memiliki urutan atau peringkat, seperti tingkat kepuasan (puas, netral, tidak puas) atau tingkat pendidikan (SMA, S1, S2, S3).

Data Biner vs Multikategori

  • Data Biner hanya memiliki dua kategori, misalnya ya/tidak, hidup/mati, sukses/gagal.
  • Data Multikategori memiliki lebih dari dua kategori, misalnya jenis pekerjaan (pegawai negeri, swasta, wiraswasta, tidak bekerja) atau status pernikahan (lajang, menikah, cerai, duda/janda).

Perbedaan dengan Data Kuantitatif

Berbeda dengan data kuantitatif yang diukur secara numerik dan dapat digunakan dalam operasi aritmatika (penjumlahan, rata-rata, dsb), data kategori lebih difokuskan pada klasifikasi. Data kuantitatif bersifat kontinu atau diskrit dan mendukung pengukuran interval dan rasio, sedangkan data kategori digunakan untuk mengelompokkan atau memberi label.

Manfaat Analisis Data Kategori dalam Berbagai Bidang

Analisis data kategori memiliki manfaat yang luas dalam berbagai sektor:

  1. Ilmu Sosial dan Psikologi
    Dalam ilmu sosial, analisis data kategori penting untuk memahami preferensi politik, perilaku masyarakat, dan respon survei. Psikologi menggunakan skala ordinal untuk mengukur tingkat depresi, kecemasan, dan aspek emosional lainnya.

  2. Kesehatan dan Kedokteran
    Dalam epidemiologi dan kesehatan masyarakat, data kategori digunakan untuk mengelompokkan pasien berdasarkan gejala, diagnosis, atau status pengobatan. Analisis dapat membantu dalam identifikasi faktor risiko penyakit tertentu.

  3. Pemasaran dan Bisnis
    Segmentasi pasar berdasarkan kategori seperti jenis kelamin, wilayah, atau preferensi produk memungkinkan perusahaan mengembangkan strategi pemasaran yang lebih tepat sasaran.

  4. Pendidikan
    Kategori seperti tingkat pendidikan, jenis sekolah, atau status siswa digunakan untuk mengevaluasi kebijakan pendidikan dan hasil belajar.

  5. Kebijakan Publik dan Pemerintahan
    Pemerintah menggunakan data kategori untuk perencanaan dan pengambilan kebijakan, seperti klasifikasi wilayah menurut tingkat kemiskinan atau pekerjaan penduduk.

  6. Keamanan dan Kriminalitas
    Data kategori digunakan untuk mengklasifikasikan jenis kejahatan, profil pelaku, dan lokasi kejadian. Analisis ini mendukung kebijakan preventif dan tindakan hukum.

Metode dalam Analisis Data Kategori

Analisis data kategori memerlukan metode khusus yang sesuai dengan karakteristik datanya. Tidak seperti data kuantitatif, variabel kategori tidak bisa dianalisis menggunakan rata-rata atau standar deviasi. Oleh karena itu, berbagai teknik telah dikembangkan untuk mengevaluasi hubungan antar kategori dan membangun model prediktif.

Tabel Kontingensi dan Uji Chi-Square

Tabel Kontingensi

Tabel kontingensi (juga disebut tabel silang) digunakan untuk menyajikan distribusi frekuensi dari dua variabel kategori. Tabel ini membantu dalam memvisualisasikan hubungan antara dua atribut. Sebagai contoh:

# Contoh data fiktif
data <- matrix(c(30, 20, 50, 40), nrow = 2, byrow = TRUE)
dimnames(data) <- list(
  Jenis_Kelamin = c("Laki-laki", "Perempuan"),
  Preferensi = c("Produk A", "Produk B")
)
tabel_kontingensi <- as.table(data)
tabel_kontingensi
##              Preferensi
## Jenis_Kelamin Produk A Produk B
##     Laki-laki       30       20
##     Perempuan       50       40

Uji Chi-Square

Uji chi-square digunakan untuk menguji apakah ada hubungan yang signifikan antara dua variabel kategori. Hipotesis nol menyatakan bahwa tidak ada asosiasi antara variabel tersebut.

# Uji Chi-Square
chisq.test(tabel_kontingensi)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabel_kontingensi
## X-squared = 0.10954, df = 1, p-value = 0.7407

Jika nilai-p lebih kecil dari 0.05, maka kita tolak hipotesis nol, yang berarti ada hubungan yang signifikan antara kedua variabel.

Regresi Logistik

Regresi logistik adalah metode prediksi untuk variabel dependen biner (misalnya: ya/tidak, sakit/sehat). Model ini mengestimasikan peluang suatu kejadian terjadi berdasarkan variabel independen kategori maupun numerik.

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

Model ini memperkirakan logit dari probabilitas am (transmisi) berdasarkan konsumsi bahan bakar dan tenaga mesin.

Analisis Korespondensi

Analisis korespondensi adalah teknik visualisasi untuk menyajikan hubungan antara dua variabel kategori dalam bentuk grafik dua dimensi. Metode ini mirip dengan analisis komponen utama tetapi khusus untuk data kategori.

# Data fiktif 3x3
tabel <- matrix(c(10, 15, 20,
                  12, 18, 25,
                  8, 10, 22),
                nrow = 3, byrow = TRUE)

colnames(tabel) <- c("Setuju", "Netral", "Tidak Setuju")
rownames(tabel) <- c("Responden A", "Responden B", "Responden C")
tabel <- as.table(tabel)

# Analisis korespondensi
library(FactoMineR)
## Warning: package 'FactoMineR' was built under R version 4.3.3
library(factoextra)
## Warning: package 'factoextra' was built under R version 4.3.3
## Loading required package: ggplot2
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
res.ca <- CA(tabel, graph = FALSE)

# Biplot berhasil
fviz_ca_biplot(res.ca, repel = TRUE)

Analisis ini berguna dalam survei dan riset pasar untuk melihat hubungan antara preferensi dan karakteristik demografis.

Decision Tree dan Random Forest

Decision Tree

Metode pohon keputusan membagi data ke dalam kelompok berdasarkan atribut-atribut kategori dan numerik, dengan membentuk struktur pohon untuk klasifikasi. Pohon ini sangat intuitif dan mudah diinterpretasikan.

library(rpart)
library(rpart.plot)
## Warning: package 'rpart.plot' was built under R version 4.3.3
# Contoh data dan model decision tree
data(iris)
tree_model <- rpart(Species ~ Sepal.Length + Sepal.Width, data = iris, method = "class")
rpart.plot(tree_model)

Random Forest

Random forest adalah pengembangan dari pohon keputusan yang menggabungkan banyak pohon (ensemble) untuk menghasilkan prediksi yang lebih akurat dan tahan terhadap overfitting.

library(randomForest)
## Warning: package 'randomForest' was built under R version 4.3.3
## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
## The following object is masked from 'package:dplyr':
## 
##     combine
# Model Random Forest
set.seed(123)
rf_model <- randomForest(Species ~ ., data = iris, ntree = 100)
print(rf_model)
## 
## Call:
##  randomForest(formula = Species ~ ., data = iris, ntree = 100) 
##                Type of random forest: classification
##                      Number of trees: 100
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 4.67%
## Confusion matrix:
##            setosa versicolor virginica class.error
## setosa         50          0         0        0.00
## versicolor      0         47         3        0.06
## virginica       0          4        46        0.08

Random forest juga dapat digunakan untuk seleksi fitur (feature importance) dalam dataset kategori maupun campuran.

Distribusi Probabilitas dalam Data Kategori

Distribusi probabilitas digunakan untuk memodelkan kemungkinan terjadinya suatu kategori atau kejadian. Dalam konteks data kategori, distribusi-distribusi ini sangat penting untuk memahami karakteristik data dan melakukan inferensi statistik.

Distribusi Bernoulli

Distribusi Bernoulli digunakan untuk memodelkan eksperimen dengan dua hasil: sukses (1) atau gagal (0). Misalnya: apakah pelanggan membeli produk (ya/tidak).

Fungsi Distribusi Bernoulli:

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

Keterangan:

  • \(X\): variabel acak bernilai 0 atau 1
  • \(p\): probabilitas sukses (X = 1)
  • \(1 - p\): probabilitas gagal (X = 0)
# Simulasi data Bernoulli
set.seed(123)
bernoulli_data <- rbinom(n = 100, size = 1, prob = 0.3)
table(bernoulli_data)
## bernoulli_data
##  0  1 
## 71 29

Distribusi Binomial

Distribusi Binomial adalah perpanjangan dari Bernoulli untuk beberapa percobaan independen. Misalnya, dari 10 pelanggan, berapa banyak yang membeli produk jika probabilitas membeli adalah 0.3.

Fungsi Distribusi Binomial:

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

Keterangan:

  • \(X\): jumlah sukses dalam \(n\) percobaan
  • \(n\): jumlah total percobaan
  • \(k\): jumlah sukses yang diamati
  • \(p\): probabilitas sukses pada setiap percobaan
  • \(\binom{n}{k}\): kombinasi dari \(n\) ambil \(k\)
# Visualisasi Distribusi Binomial
x <- 0:10
prob <- dbinom(x, size = 10, prob = 0.3)
barplot(prob, names.arg = x, main = "Distribusi Binomial (n=10, p=0.3)", ylab = "Probabilitas")

Distribusi Multinomial

Distribusi Multinomial digunakan untuk eksperimen dengan lebih dari dua kategori. Contohnya: memilih satu dari tiga rasa es krim – coklat, vanila, atau stroberi.

Fungsi Distribusi Multinomial:

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

Keterangan:

  • \(n\): jumlah total percobaan atau uji coba
  • \(X_i\): jumlah kemunculan kategori ke-\(i\)
  • \(x_i\): realisasi atau nilai dari \(X_i\)
  • \(p_i\): probabilitas kemunculan kategori ke-\(i\)
  • \(k\): jumlah kategori yang mungkin
# Paket yang dibutuhkan
library(extraDistr)
## Warning: package 'extraDistr' was built under R version 4.3.3
# Simulasi dari distribusi multinomial
prob <- c(0.3, 0.5, 0.2)
set.seed(123)
rmultinom(n = 1, size = 10, prob = prob)
##      [,1]
## [1,]    2
## [2,]    5
## [3,]    3

Distribusi Poisson

Distribusi Poisson digunakan untuk menghitung jumlah kejadian dalam interval waktu atau ruang tertentu. Misalnya: jumlah keluhan pelanggan per jam.

Fungsi Distribusi Poisson:

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

Keterangan:

  • \(X\): jumlah kejadian dalam satu interval
  • \(k\): nilai tertentu dari jumlah kejadian (0,1,2,…)
  • \(\lambda\): rata-rata kejadian per interval
  • \(e\): konstanta Euler (\(\approx 2.718\))
# Distribusi Poisson
x <- 0:10
pois_prob <- dpois(x, lambda = 2)
barplot(pois_prob, names.arg = x, main = "Distribusi Poisson (lambda = 2)", ylab = "Probabilitas")

Distribusi-distribusi ini merupakan landasan penting dalam membangun model statistik untuk data kategori.

Desain Sampling dalam Analisis Data Kategori

Desain sampling adalah aspek penting dalam analisis data kategori karena menentukan bagaimana data dikumpulkan sehingga hasil analisis valid dan dapat digeneralisasi.

Prospective Sampling

Desain sampling prospektif mengamati subjek dari titik awal dan mengikuti ke depan dalam waktu untuk melihat hasil yang terjadi.

Eksperimen

Eksperimen adalah desain di mana peneliti mengatur kondisi atau perlakuan pada subjek dan kemudian mengamati efeknya. Contohnya uji klinis obat baru di mana subjek dibagi dalam kelompok kontrol dan perlakuan.

  • Keuntungan: Kontrol penuh terhadap variabel, dapat menentukan sebab-akibat.
  • Kekurangan: Biaya dan waktu tinggi, kadang tidak etis untuk beberapa kasus.

Studi Kohort

Studi kohort mengikuti sekelompok orang yang memiliki karakteristik tertentu selama periode waktu untuk melihat apakah mereka mengalami kejadian tertentu (misal: penyakit).

  • Keuntungan: Mengukur insidensi, waktu kejadian jelas.
  • Kekurangan: Perlu waktu lama, kemungkinan kehilangan follow-up.

Retrospective Sampling

Desain sampling retrospektif melihat ke belakang pada data atau kejadian yang sudah terjadi.

Studi Kasus Kontrol

Desain yang membandingkan kelompok dengan kondisi (kasus) dan tanpa kondisi (kontrol) untuk menemukan faktor risiko.

  • Keuntungan: Efisien waktu dan biaya.
  • Kekurangan: Bias recall, sulit menentukan sebab-akibat.

Studi Kohort Retrospektif

Mirip studi kohort, tetapi data dikumpulkan dari catatan masa lalu.

  • Keuntungan: Waktu pengumpulan data lebih cepat dibanding prospektif.
  • Kekurangan: Data bisa tidak lengkap atau bias.

Tabel Perbandingan Desain Sampling

Desain Sampling Keuntungan Kekurangan Contoh Aplikasi
Eksperimen Kontrol penuh, sebab-akibat jelas Mahal, kadang tidak etis Uji klinis obat
Studi Kohort Prospektif Dapat mengukur insidensi, waktu jelas Lama dan mahal Studi penyakit kronis
Studi Kasus Kontrol Efisien waktu dan biaya Bias recall, sebab-akibat sulit Faktor risiko kanker
Studi Kohort Retrospektif Cepat pengumpulan data Data tidak lengkap, bias Studi faktor risiko menggunakan rekam medis

Desain sampling yang tepat sangat menentukan kualitas dan validitas analisis data kategori, serta interpretasi hasil yang benar.

Tabel Kontingensi 2x2

Tabel kontingensi 2x2 adalah alat dasar dalam analisis data kategori untuk melihat hubungan antara dua variabel kategori yang masing-masing memiliki dua level.

Studi Kasus: Penyediaan Air Bersih dan Kejadian Diare

Data berikut diperoleh dari studi yang mengevaluasi pengaruh penyediaan air bersih terhadap kejadian diare pada suatu populasi:

Penyediaan Air Bersih Diare Tidak Diare Total
Tidak Memenuhi Syarat 36 7 43
Memenuhi Syarat 14 15 29
Total 50 22 72

Distribusi Peluang dalam Tabel Kontingensi 2x2

Distribusi peluang mencakup peluang bersama, marginal, dan bersyarat dari kejadian pada tabel.

Peluang Bersama

Peluang bersama adalah probabilitas bahwa dua kejadian terjadi bersama-sama.

Contoh: \[ P(\text{Diare} \cap \text{Tidak Memenuhi}) = \frac{36}{72} = 0.5 \]

Peluang Marginal

\[ P(\text{Diare}) = \frac{50}{72}, \quad P(\text{Tidak Diare}) = \frac{22}{72} \]

\[ P(\text{Tidak Memenuhi}) = \frac{43}{72}, \quad P(\text{Memenuhi}) = \frac{29}{72} \]

Peluang Bersyarat

\[ P(\text{Diare} | \text{Tidak Memenuhi}) = \frac{36}{43}, \quad P(\text{Diare} | \text{Memenuhi}) = \frac{14}{29} \]


Ukuran Asosiasi dalam Data Kategori 2x2

Ukuran asosiasi penting untuk mengukur kekuatan hubungan antara dua variabel kategori.

Misalkan tabel frekuensi (jumlah) data 2x2 adalah:

Diare (\(Y\)) Tidak Diare (\(Y^c\)) Total
Tidak Memenuhi Syarat (\(X\)) \(a=36\) \(b=7\) \(a+b=43\)
Memenuhi Syarat (\(X^c\)) \(c=14\) \(d=15\) \(c+d=29\)
Total \(a+c=50\) \(b+d=22\) \(n=72\)

Risk Difference (RD)

\[ RD = \frac{a}{a+b} - \frac{c}{c+d} = \frac{36}{43} - \frac{14}{29} \approx 0.8372 - 0.4828 = 0.3544 \]

Relative Risk (RR)

\[ RR = \frac{\frac{36}{43}}{\frac{14}{29}} = \frac{0.8372}{0.4828} \approx 1.7345 \]

Odds Ratio (OR)

\[ OR = \frac{ad}{bc} = \frac{36 \times 15}{7 \times 14} = \frac{540}{98} \approx 5.51 \]


Implementasi R

# Data kasus diare dan penyediaan air bersih
a <- 36  # Diare, Tidak Memenuhi
b <- 7   # Tidak Diare, Tidak Memenuhi
c <- 14  # Diare, Memenuhi
d <- 15  # Tidak Diare, Memenuhi

# Risk Difference
RD <- (a / (a + b)) - (c / (c + d))

# Relative Risk
RR <- (a / (a + b)) / (c / (c + d))

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

# Output
list(Risk_Difference = RD,
     Relative_Risk = RR,
     Odds_Ratio = OR)
## $Risk_Difference
## [1] 0.3544507
## 
## $Relative_Risk
## [1] 1.734219
## 
## $Odds_Ratio
## [1] 5.510204

Interpretasi

  • Risk Difference (RD = 0.354): Risiko diare pada kelompok dengan penyediaan air tidak memenuhi syarat lebih tinggi sebesar 35,4% dibandingkan kelompok dengan air memenuhi syarat.
  • Relative Risk (RR = 1.73): Individu yang tinggal di tempat dengan penyediaan air yang tidak memenuhi syarat memiliki 1,73 kali lipat risiko terkena diare dibandingkan dengan yang tinggal di tempat dengan air memenuhi syarat.
  • Odds Ratio (OR = 5.51): Peluang terkena diare di kelompok tidak memenuhi syarat sekitar 5,5 kali lebih besar dibandingkan kelompok memenuhi syarat. ————————————————————————

Perbandingan RD, RR, dan OR

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

Kesimpulan:

  • RD menunjukkan bahwa tidak tersedianya air bersih meningkatkan risiko absolut terkena diare sebesar 35.44%.
  • RR mengindikasikan bahwa risiko terkena diare pada mereka yang tidak mendapat air bersih adalah sekitar 1.73 kali lebih tinggi.
  • OR sebesar 5.51 menyatakan bahwa odds terkena diare pada kelompok tanpa air bersih sekitar 5.5 kali lipat dibanding yang mendapat air bersih.

Sumber Data: https://www.academia.edu/89066874/Hubungan_Penyediaan_Air_Bersih_Dan_Penggunaan_Jamban_Dengan_Kejadian_Diare_Pada_Balita_DI_Puskesmas_Pasar_Ikan_Kota_Bengkulu?utm_source=chatgpt.com#loswp-work-container

Inferensi Tabel Kontingensi Dua Arah

Estimasi

Estimasi Titik

Estimasi titik digunakan untuk mendapatkan nilai terbaik dari parameter populasi berdasarkan data sampel.

Estimasi Interval

Estimasi interval memberikan rentang nilai yang kemungkinan besar mencakup parameter populasi.

Uji Hipotesis

Uji Proporsi

Uji proporsi dilakukan untuk membandingkan proporsi antar kelompok dalam tabel kontingensi 2x2.

Uji Asosiasi

Risk Difference (RD)

Risk Difference mengukur perbedaan absolut dalam probabilitas kejadian dua kelompok:

\[ RD = \frac{a}{a + b} - \frac{c}{c + d} \]

Standard Error untuk RD: \[ SE(RD) = \sqrt{\frac{a(a+b-a)}{(a+b)^3} + \frac{c(c+d-c)}{(c+d)^3}} \]

Statistik uji Z: \[ Z = \frac{RD}{SE(RD)} \]

Relative Risk (RR)

Relative Risk membandingkan kemungkinan kejadian antara dua kelompok:

\[ RR = \frac{a/(a + b)}{c/(c + d)} \]

Standard Error untuk log(RR): \[ SE(\log RR) = \sqrt{\frac{1}{a} - \frac{1}{a + b} + \frac{1}{c} - \frac{1}{c + d}} \]

Statistik uji Z: \[ Z = \frac{\log(RR)}{SE(\log RR)} \]

Odds Ratio (OR)

Odds Ratio membandingkan peluang kejadian antara dua kelompok:

\[ OR = \frac{a \times d}{b \times c} \]

Standard Error untuk log(OR): \[ SE(\log OR) = \sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}} \]

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

Contoh Perhitungan Manual

Misalkan data dari studi kasus air bersih dan diare:

  • a = 36 (Diare, Tidak Memenuhi)
  • b = 7 (Tidak Diare, Tidak Memenuhi)
  • c = 14 (Diare, Memenuhi)
  • d = 15 (Tidak Diare, Memenuhi)
a <- 36
b <- 7
c <- 14
d <- 15

# RD
RD <- (a / (a + b)) - (c / (c + d))
SE_RD <- sqrt((a*(a+b-a)/(a+b)^3) + (c*(c+d-c)/(c+d)^3))
Z_RD <- RD / SE_RD

# RR
RR <- (a / (a + b)) / (c / (c + d))
SE_log_RR <- sqrt(1/a - 1/(a + b) + 1/c - 1/(c + d))
Z_RR <- log(RR) / SE_log_RR

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

list(RD = RD, SE_RD = SE_RD, Z_RD = Z_RD,
     RR = RR, SE_log_RR = SE_log_RR, Z_RR = Z_RR,
     OR = OR, SE_log_OR = SE_log_OR, Z_OR = Z_OR)
## $RD
## [1] 0.3544507
## 
## $SE_RD
## [1] 0.1085356
## 
## $Z_RD
## [1] 3.265756
## 
## $RR
## [1] 1.734219
## 
## $SE_log_RR
## [1] 0.2036364
## 
## $Z_RR
## [1] 2.703629
## 
## $OR
## [1] 5.510204
## 
## $SE_log_OR
## [1] 0.5556349
## 
## $Z_OR
## [1] 3.071444

Interpretasi

  • Risk Difference (RD = 0.354, Z = 3.27): Terdapat perbedaan risiko sebesar 35,4% antara kelompok dengan air tidak memenuhi syarat dan kelompok yang memenuhi. Nilai Z = 3.27 menunjukkan bahwa perbedaan ini signifikan secara statistik (Z > 1.96 pada α = 0.05).
  • Relative Risk (RR = 1.73, Z = 2.70): Risiko diare pada kelompok air tidak layak adalah 1,73 kali lebih tinggi dibanding kelompok air layak. Nilai Z menunjukkan bahwa RR ini signifikan secara statistik.

  • Odds Ratio (OR = 5.51, Z = 3.07): Peluang terkena diare 5,5 kali lebih besar pada kelompok dengan air tidak memenuhi syarat. Nilai Z > 1.96 menunjukkan bahwa asosiasi ini signifikan.

Uji Independensi

Uji Chi-Square

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

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

  • O = observed (frekuensi aktual)
  • E = expected (frekuensi harapan)

Contoh Perhitungan Chi-Square di R

observed <- matrix(c(36, 7, 14, 15), nrow = 2, byrow = TRUE)
chisq.test(observed, correct = FALSE)
## 
##  Pearson's Chi-squared test
## 
## data:  observed
## X-squared = 10.255, df = 1, p-value = 0.001363

Interpretasi

Tolak H₀ — artinya terdapat hubungan yang signifikan antara penyediaan air bersih dan kejadian diare. Dengan kata lain, penyediaan air bersih mempengaruhi kejadian diare secara statistik.

Partisi Chi-Square

Digunakan untuk memeriksa kontribusi tiap sel terhadap nilai chi-square total.

Uji Likelihood Ratio (G^2)

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

Uji Exact Fisher

Digunakan jika ukuran sampel kecil.

fisher.test(observed)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  observed
## p-value = 0.001892
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   1.647619 19.217286
## sample estimates:
## odds ratio 
##   5.363204

Interpretasi

Tolak H₀ — terdapat perbedaan signifikan dalam odds kejadian diare antara kedua kelompok. Odds ratio sebesar 5.36 berarti peluang terkena diare pada kelompok yang tidak mendapat air bersih sekitar 5.4 kali lebih besar dibanding kelompok yang mendapat air bersih. Karena interval kepercayaan tidak mencakup 1, maka hasilnya signifikan secara statistik.

Analisis Residual dalam Tabel Kontingensi

Jenis Residual

  1. Pearson Residual \[ R_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}} \]

  2. Standardized Residual Digunakan untuk membandingkan residual antar sel.

Deteksi Outlier dalam Analisis Data Kategori Menggunakan Residual

Apa Itu Outlier dalam Tabel Kontingensi?

Sel dengan residual besar (misalnya > 2 atau < -2) dianggap sebagai outlier.

Menggunakan Residual untuk Mendeteksi Outlier

chisq <- chisq.test(observed)
chisq$residuals
##           [,1]      [,2]
## [1,]  1.123406 -1.693598
## [2,] -1.367956  2.062271

Interpretasi

  • Sel (2,2) memiliki residual 2.06, menunjukkan bahwa jumlah aktualnya jauh lebih besar dari yang diharapkan, dan merupakan kontributor utama terhadap signifikansi uji chi-square.
  • Sebaliknya, sel (1,2) dan (2,1) memiliki residual negatif, menunjukkan bahwa observasi aktual lebih sedikit dari yang diharapkan.

  • Sel (1,1) tidak terlalu menyimpang.

Kesimpulan

Bab ini membahas teknik inferensi statistik untuk tabel kontingensi dua arah, termasuk estimasi parameter, uji hipotesis terhadap ukuran asosiasi (RD, RR, OR), dan berbagai uji independensi (Chi-square, G^2, Fisher). Analisis residual juga memberikan wawasan lebih dalam untuk mengidentifikasi sel yang menyimpang.

Sumber Data: https://www.academia.edu/89066874/Hubungan_Penyediaan_Air_Bersih_Dan_Penggunaan_Jamban_Dengan_Kejadian_Diare_Pada_Balita_DI_Puskesmas_Pasar_Ikan_Kota_Bengkulu?utm_source=chatgpt.com#loswp-work-container

Tabel Kontingensi Tiga Arah

Tabel kontingensi tiga arah adalah perpanjangan dari tabel kontingensi dua arah yang digunakan untuk menganalisis hubungan antara tiga variabel kategori secara simultan. Dalam banyak situasi, hubungan antara dua variabel (misalnya \(X\) dan \(Y\)) dapat dipengaruhi oleh variabel ketiga \(Z\), yang disebut sebagai variabel kontrol atau kovariat.

Contoh penggunaan:

Tabel ini dapat dipecah menjadi: 1. Tabel Parsial: Tabel dua arah yang dikondisikan pada satu nilai dari variabel kontrol \(Z\). 2. Tabel Marginal: Tabel dua arah yang didapat dengan menjumlahkan seluruh nilai \(Z\).

Tabel Parsial dan Marginal

Tabel parsial memungkinkan kita mengevaluasi hubungan antara \(X\) dan \(Y\) pada setiap tingkat \(Z\). Tabel marginal hanya memperlihatkan hubungan keseluruhan antara \(X\) dan \(Y\) tanpa mempertimbangkan \(Z\).

Data: Pendidikan (X), Wilayah (Z), Merokok (Y)

Pendidikan (X) Wilayah (Z) Tidak Pernah Pernah Aktif
SD Kota 12 10 8
Desa 20 14 9
Pinggiran 18 16 12
SMP Kota 15 11 9
Desa 17 10 11
Pinggiran 13 14 10
SMA Kota 20 15 12
Desa 12 18 10
Pinggiran 10 16 13

Distribusi Peluang

Peluang Bersama

Rumus peluang bersama untuk tiga variabel: \[ P(X = x, Y = y, Z = z) = \frac{n_{xyz}}{N} \] Dimana:

  • \(n_{xyz}\) adalah frekuensi dari kombinasi (\(x, y, z\))

  • \(N\) adalah total observasi

Contoh manual: \[ P(\text{SMP, Aktif, Desa}) = \frac{11}{297} \approx 0.0370 \]

Dengan R:

# Buat array 3D
array3d <- array(
  c(12,10,8, 20,14,9, 18,16,12,  # SD
    15,11,9, 17,10,11, 13,14,10, # SMP
    20,15,12, 12,18,10, 10,16,13),
  dim = c(3,3,3),
  dimnames = list(
    Merokok = c("Tidak", "Pernah", "Aktif"),
    Wilayah = c("Kota", "Desa", "Pinggiran"),
    Pendidikan = c("SD", "SMP", "SMA")
  )
)

# Peluang bersama
prop.table(array3d)
## , , Pendidikan = SD
## 
##         Wilayah
## Merokok        Kota       Desa  Pinggiran
##   Tidak  0.03380282 0.05633803 0.05070423
##   Pernah 0.02816901 0.03943662 0.04507042
##   Aktif  0.02253521 0.02535211 0.03380282
## 
## , , Pendidikan = SMP
## 
##         Wilayah
## Merokok        Kota       Desa  Pinggiran
##   Tidak  0.04225352 0.04788732 0.03661972
##   Pernah 0.03098592 0.02816901 0.03943662
##   Aktif  0.02535211 0.03098592 0.02816901
## 
## , , Pendidikan = SMA
## 
##         Wilayah
## Merokok        Kota       Desa  Pinggiran
##   Tidak  0.05633803 0.03380282 0.02816901
##   Pernah 0.04225352 0.05070423 0.04507042
##   Aktif  0.03380282 0.02816901 0.03661972

Kesimpulan Output

  • Responden terbanyak tampaknya berasal dari kategori: Tidak merokok, tinggal di Desa, pendidikan SD (5.63%) Tidak merokok, tinggal di Kota, pendidikan SMA (5.63%)
  • Merokok aktif cenderung memiliki proporsi lebih kecil di hampir semua kategori dibandingkan “Tidak” atau “Pernah”.

  • Perbedaan proporsi ini bisa jadi bahan awal analisis lebih lanjut, misalnya dengan uji chi-square untuk melihat ada tidaknya asosiasi antara variabel Merokok, Wilayah, dan Pendidikan.

Peluang Marginal

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

Manual contoh: \[ P(\text{SMA}) = \frac{20+15+12 + 12+18+10 + 10+16+13}{297} = \frac{136}{297} \approx 0.458 \]

Dengan R:

margin.table(array3d, 3) / sum(array3d)  # dim = 3 untuk Pendidikan
## Pendidikan
##        SD       SMP       SMA 
## 0.3352113 0.3098592 0.3549296

Interpretasi

  • Sekitar 33.5% dari seluruh responden berpendidikan SD
  • Sekitar 31.0% berpendidikan SMP

  • Sekitar 35.5% berpendidikan SMA

Peluang Bersyarat

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

Dengan R:

prop.table(margin.table(array3d, c(2,3)), 1)  # Kondisional pada Wilayah
##            Pendidikan
## Wilayah            SD       SMP       SMA
##   Kota      0.2678571 0.3125000 0.4196429
##   Desa      0.3553719 0.3140496 0.3305785
##   Pinggiran 0.3770492 0.3032787 0.3196721

Interpretasi

  • Untuk Wilayah Kota, sekitar:
    • 26.8% responden berpendidikan SD
    • 31.3% berpendidikan SMP
    • 41.96% berpendidikan SMA Artinya, di kota, proporsi pendidikan SMA adalah yang tertinggi.
  • Untuk Wilayah Desa, distribusi pendidikan relatif lebih merata:
    • 35.5% SD
    • 31.4% SMP
    • 33.1% SMA
  • Untuk Wilayah Pinggiran, distribusi pendidikan mirip dengan desa, tapi SD sedikit lebih tinggi:
    • 37.7% SD
    • 30.3% SMP
    • 32.0% SMA

Tabel Peluang Bersyarat

Tabel ini menunjukkan probabilitas kondisi dua variabel diberikan nilai tertentu dari variabel ketiga.

Ukuran Asosiasi

Tabel Kontingensi Parsial

Ukuran asosiasi seperti Odds Ratio (OR) dapat dihitung pada masing-masing strata dari variabel \(Z\): \[ OR = \frac{a/c}{b/d} = \frac{ad}{bc} \]

Contoh: Misalkan di Wilayah Kota untuk SD:

  • Tidak Pernah = a = 12
  • Pernah = b = 10
  • Aktif = c = 8

Ubah ke 2x2 tabel misalnya: Tidak Pernah vs Aktif: \[ OR = \frac{12*1}{10*8} = \frac{12}{80} = 0.15 \]

Dengan R:

library(epitools)
# OR untuk tiap strata Wilayah (misal untuk SD):
oddsratio(array3d[,,"SD"])
## $data
##         Wilayah
## Merokok  Kota Desa Pinggiran Total
##   Tidak    12   20        18    50
##   Pernah   10   14        16    40
##   Aktif     8    9        12    29
##   Total    30   43        46   119
## 
## $measure
##         odds ratio with 95% C.I.
## Merokok   estimate     lower    upper
##   Tidak  1.0000000        NA       NA
##   Pernah 0.8424144 0.2795678 2.545883
##   Aktif  0.6808146 0.2002842 2.308871
## 
## $p.value
##         two-sided
## Merokok  midp.exact fisher.exact chi.square
##   Tidak          NA           NA         NA
##   Pernah  0.7595681    0.8948267  0.8823789
##   Aktif   0.5344601    0.7305725  0.7279645
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"

Conditional Independence

Dua variabel dikatakan conditionally independent jika tidak ada hubungan antara keduanya setelah mengendalikan variabel ketiga.

Marginal Y dan X

Terkadang hubungan marginal antara \(X\) dan \(Y\) bisa berbeda dari hubungan pada setiap strata \(Z\), fenomena ini disebut Simpson’s paradox.

Inferensi Tabel Kontingensi Tiga Arah

Independensi Bersyarat dalam Tabel Kontingensi Tiga Arah

Untuk menguji apakah \(X\) dan \(Y\) independen secara bersyarat terhadap \(Z\).

Pengujian Statistik untuk Independensi Bersyarat

Uji Chi-Square dapat dilakukan pada setiap strata \(Z\) dan dibandingkan.

Odds Ratio Bersama

Odds Ratio dihitung pada setiap strata dan kemudian digabungkan menggunakan metode Mantel-Haenszel:

\[ OR_{MH} = \frac{\sum_i \frac{a_i d_i}{n_i}}{\sum_i \frac{b_i c_i}{n_i}} \]

Penjelasan:

  • \(a_i\), \(b_i\), \(c_i\), \(d_i\): Frekuensi dalam sel 2x2 untuk strata ke-\(i\)
  • \(n_i\): Total pada strata ke-\(i\)

Uji Homogenitas Odds Ratio dengan Statistik Breslow-Day

Digunakan untuk menguji apakah OR konsisten di semua strata:

# Buat array 3D seperti sebelumnya
array3d <- array(
  c(12,10,8, 20,14,9, 18,16,12,  # SD
    15,11,9, 17,10,11, 13,14,10, # SMP
    20,15,12, 12,18,10, 10,16,13), # SMA
  dim = c(3,3,3),
  dimnames = list(
    Merokok = c("Tidak", "Pernah", "Aktif"),
    Wilayah = c("Kota", "Desa", "Pinggiran"),
    Pendidikan = c("SD", "SMP", "SMA")
  )
)

# Ambil subset 2x2x3: hanya 2 level Merokok dan Wilayah → Pendidikan sebagai strata
sub_array <- array3d[c("Tidak", "Aktif"), c("Kota", "Desa"), ]

# Fungsi manual Breslow-Day
breslow_day_manual <- function(array_2x2xk) {
  if (!all(dim(array_2x2xk)[1:2] == c(2, 2))) {
    stop("Array harus berukuran 2x2xk")
  }

  k <- dim(array_2x2xk)[3]
  OR_vec <- numeric(k)
  var_logOR <- numeric(k)

  for (i in 1:k) {
    a <- array_2x2xk[1, 1, i]
    b <- array_2x2xk[1, 2, i]
    c <- array_2x2xk[2, 1, i]
    d <- array_2x2xk[2, 2, i]

    OR_vec[i] <- (a * d) / (b * c)
    var_logOR[i] <- (1 / a) + (1 / b) + (1 / c) + (1 / d)
  }

  MH_num <- sum((array_2x2xk[1,1,]*array_2x2xk[2,2,]) / apply(array_2x2xk, 3, sum))
  MH_den <- sum((array_2x2xk[1,2,]*array_2x2xk[2,1,]) / apply(array_2x2xk, 3, sum))
  OR_MH <- MH_num / MH_den
  logOR_MH <- log(OR_MH)

  BD_stat <- sum((log(OR_vec) - logOR_MH)^2 / var_logOR)
  df <- k - 1
  p_val <- 1 - pchisq(BD_stat, df)

  return(list(
    OR_per_stratum = OR_vec,
    logOR_MH = logOR_MH,
    BD_statistic = BD_stat,
    df = df,
    p_value = p_val
  ))
}

# Jalankan uji
bd_test <- breslow_day_manual(sub_array)
bd_test
## $OR_per_stratum
## [1] 0.675000 1.078431 1.388889
## 
## $logOR_MH
## [1] 0.02301189
## 
## $BD_statistic
## [1] 0.7716713
## 
## $df
## [1] 2
## 
## $p_value
## [1] 0.6798822

Visualisasi

library(ggplot2)
library(dplyr)
# Ubah array 3d jadi data frame long
mydata <- as.data.frame.table(array3d, responseName = "Freq")

# Cek data
head(mydata)
##   Merokok Wilayah Pendidikan Freq
## 1   Tidak    Kota         SD   12
## 2  Pernah    Kota         SD   10
## 3   Aktif    Kota         SD    8
## 4   Tidak    Desa         SD   20
## 5  Pernah    Desa         SD   14
## 6   Aktif    Desa         SD    9
#   Merokok Wilayah Pendidikan Freq
# 1  Tidak    Kota         SD   12
# 2  Pernah   Kota         SD   10
# 3  Aktif    Kota         SD    8
# ... dst

# Plot dengan ggplot2
library(ggplot2)

ggplot(mydata, aes(x = Merokok, y = Freq, fill = Pendidikan)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~ Wilayah) +
  labs(title = "Distribusi Perokok berdasarkan Pendidikan dan Wilayah",
       y = "Frekuensi") +
  theme_minimal()

Kesimpulan

Tabel kontingensi tiga arah memungkinkan kita mengevaluasi hubungan dua variabel dengan kontrol terhadap variabel ketiga. Uji seperti Breslow-Day dan Mantel-Haenszel penting untuk menentukan signifikansi dan homogenitas asosiasi di seluruh strata. Dengan data baru ini, kita dapat menggali lebih lanjut pengaruh pendidikan dan wilayah terhadap kebiasaan merokok.

Sumber Data: Dummy, karena kesulitan menemukan data dengan tabel kontingensi 3x3

Generalized Linear Model (GLM)

Generalized Linear Model (GLM) adalah kerangka pemodelan statistik yang memperluas regresi linear klasik dengan memungkinkan: - Respon yang tidak harus kontinu - Distribusi dari keluarga eksponensial (seperti Binomial, Poisson, dan Gaussian) - Penggunaan fungsi link yang menghubungkan rata-rata dari variabel respon ke prediktor linear

Exponential Family

Distribusi yang termasuk keluarga eksponensial dapat dinyatakan dalam bentuk umum:

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

Dimana: - \(\theta\): parameter kanonik - \(\phi\): parameter dispersi - \(a(.), b(.), c(.)\): fungsi tertentu tergantung pada distribusi

Contoh distribusi dari keluarga eksponensial: - Normal (Gaussian) - Binomial (untuk regresi logistik) - Poisson (untuk regresi hitung)

Model Regresi Logistik

Model regresi logistik digunakan untuk memodelkan probabilitas dari suatu kejadian biner. Model ini menggunakan fungsi link logit:

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

Dimana: - \(\pi(x)\): probabilitas kejadian - \(x\): prediktor - \(\beta_0, \beta_1\): parameter model

Contoh dan Implementasi R

set.seed(123)
x <- rnorm(100)
logit_p <- 1 / (1 + exp(-(1 + 2 * x)))
y <- rbinom(100, size = 1, prob = logit_p)
model_logit <- glm(y ~ x, family = binomial)
summary(model_logit)
## 
## Call:
## glm(formula = y ~ x, family = binomial)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.0434     0.3063   3.407 0.000657 ***
## x             2.4414     0.5241   4.658  3.2e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 126.836  on 99  degrees of freedom
## Residual deviance:  81.452  on 98  degrees of freedom
## AIC: 85.452
## 
## Number of Fisher Scoring iterations: 6
# Plot
plot(x, y, main = "Regresi Logistik")
curve(predict(model_logit, newdata = data.frame(x = x), type = "response"), 
      col = "red", lwd = 2, add = TRUE)

Interpretasi - Kurva berbentuk-S (sigmoid) adalah ciri khas dari model regresi logistik. - Probabilitas mendekati 0 ketika𝑥sangat kecil (di sisi kiri), dan mendekati 1 ketika𝑥sangat besar (di sisi kanan). - Transisi tajam terjadi di sekitar nilai𝑥=0, menunjukkan bahwa sekitar titik ini, perubahan nilai𝑥sangat berpengaruh terhadap probabilitas keberhasilan. - Model ini menunjukkan hubungan positif: semakin besar nilai𝑥, semakin besar peluang terjadinya kejadian (y = 1).

Model Regresi Poisson

Model regresi Poisson digunakan untuk data hitungan (count), misalnya jumlah kejadian per satuan waktu atau ruang. Model ini menggunakan link log:

\[ \log(\lambda(x)) = \beta_0 + \beta_1 x \]

Dimana: - \(\lambda(x)\): rata-rata jumlah kejadian - \(x\): prediktor

Contoh dan Implementasi R

set.seed(123)
x <- rnorm(100)
lambda <- exp(1 + 0.5 * x)
y <- rpois(100, lambda = lambda)
model_pois <- glm(y ~ x, family = poisson)
summary(model_pois)
## 
## Call:
## glm(formula = y ~ x, family = poisson)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.96672    0.06596  14.657   <2e-16 ***
## x            0.54847    0.06237   8.794   <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: 176.656  on 99  degrees of freedom
## Residual deviance:  98.658  on 98  degrees of freedom
## AIC: 373.32
## 
## Number of Fisher Scoring iterations: 5
# Plot
plot(x, y, main = "Regresi Poisson")
lines(sort(x), fitted(model_pois)[order(x)], col = "blue", lwd = 2)

Interpretasi: - Garis biru menunjukkan nilai prediksi dari model Poisson. - Model memprediksi bahwa rata-rata jumlah kejadian (count) meningkat eksponensial terhadap nilai x. - Sebagian besar titik data berada dekat garis prediksi, menunjukkan model cocok. - Ketidakteraturan di bagian atas mungkin menunjukkan potensi overdispersion (varian melebihi mean).

Kesimpulan

GLM memungkinkan pemodelan fleksibel untuk data dengan distribusi berbeda dan berbagai tipe respon. Dua aplikasi penting dari GLM adalah: - Regresi Logistik: memodelkan probabilitas kejadian biner - Regresi Poisson: memodelkan data hitungan

Inferensi Generalized Linear Model (GLM)

Bab ini membahas bagaimana kita melakukan inferensi statistik terhadap model GLM, termasuk penaksiran parameter, pengujian signifikansi, evaluasi kualitas model, dan metode diagnostik. Inferensi ini penting untuk memastikan bahwa model yang dibangun valid dan sesuai dengan data.

Mencari Ekspektasi dan Varians dalam GLM

Dalam kerangka GLM, ekspektasi dan variansi dari variabel respon \(Y\) tergantung pada bentuk distribusi dari keluarga eksponensial dan fungsi link yang digunakan.

  • Ekspektasi dari \(Y\): \[ E[Y] = \mu = g^{-1}(\eta) \]
  • Variansi dari \(Y\): \[ Var(Y) = \phi V(\mu) \]

Dimana: - \(\eta = X \beta\): prediktor linear - \(g^{-1}\): fungsi link invers - \(\phi\): parameter dispersi (untuk Binomial dan Poisson biasanya = 1) - \(V(\mu)\): fungsi varians spesifik distribusi

Contoh Fungsi Varians: - Binomial: \(V(\mu) = \mu(1 - \mu)\) - Poisson: \(V(\mu) = \mu\) - Gaussian: \(V(\mu) = \sigma^2\)

Metode Penaksiran Parameter

Maximum Likelihood Estimation (MLE)

GLM menggunakan metode Maximum Likelihood Estimation untuk mencari nilai parameter \(\beta\) yang memaksimalkan fungsi log-likelihood:

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

Proses estimasi diselesaikan secara numerik menggunakan Iteratively Reweighted Least Squares (IRLS).

Implementasi di R

Fungsi glm() dalam R secara otomatis menggunakan IRLS untuk melakukan estimasi MLE.

set.seed(123)
x <- rnorm(100)
logit_p <- 1 / (1 + exp(-(1 + 2 * x)))
y <- rbinom(100, size = 1, prob = logit_p)
model_logit <- glm(y ~ x, family = binomial)
summary(model_logit)
## 
## Call:
## glm(formula = y ~ x, family = binomial)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.0434     0.3063   3.407 0.000657 ***
## x             2.4414     0.5241   4.658  3.2e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 126.836  on 99  degrees of freedom
## Residual deviance:  81.452  on 98  degrees of freedom
## AIC: 85.452
## 
## Number of Fisher Scoring iterations: 6

Koefisien: - Intercept (β₀ = 1.0434): log-odds kejadian saat \(x = 0\). Probabilitas kejadian dasar dapat dihitung sebagai: \[ \pi_0 = \frac{e^{1.0434}}{1 + e^{1.0434}} \approx 0.74 \] - x (β₁ = 2.4414): setiap kenaikan 1 unit pada \(x\) meningkatkan log-odds sebesar 2.44. Dalam probabilitas, ini berarti peningkatan substansial pada peluang kejadian (karena koefisien positif dan signifikan).

Statistik: - Nilai z untuk masing-masing koefisien jauh di atas 2, dan p-value < 0.001 menunjukkan kedua parameter signifikan. - Null deviance: 126.84 (model tanpa prediktor) - Residual deviance: 81.45 (model dengan prediktor) → penurunan ini menunjukkan model membaik - AIC = 85.45: digunakan untuk membandingkan antar model

Diagnostik Model GLM

Setelah membangun model, perlu dilakukan pemeriksaan diagnostik untuk mengevaluasi: - Kesesuaian model terhadap data - Deteksi outlier atau leverage points - Validitas asumsi model

Statistik Diagnostik Umum:

  • Deviance: mengukur ketidaksesuaian antara model dan data
  • AIC (Akaike Information Criterion): membantu memilih model terbaik
  • Residuals:
    • Pearson residual: \(r_i = \frac{y_i - \hat{\mu}_i}{\sqrt{Var(y_i)}}\)
    • Deviance residual: turunan dari log-likelihood

Visualisasi Diagnostik:

plot(model_logit)  # Menampilkan 4 grafik diagnostik penting

Detail Estimasi dan Inferensi Regresi Logistik

Regresi logistik digunakan untuk memodelkan variabel biner. Untuk menguji signifikansi model, digunakan likelihood ratio test dengan statistik:

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

  • \(L_0\): log-likelihood model null (hanya intercept)
  • \(L_1\): log-likelihood model penuh

Contoh Penghitungan Manual di R:

logLik_null <- logLik(glm(y ~ 1, family = binomial))
logLik_full <- logLik(model_logit)
G2 <- -2 * (as.numeric(logLik_null) - as.numeric(logLik_full))
G2
## [1] 45.38362
pchisq(G2, df = 1, lower.tail = FALSE)  # p-value
## [1] 1.619833e-11

Interpretasi:

  • Jika nilai \(p < 0.05\), maka variabel prediktor berpengaruh signifikan.
  • G² digunakan sebagai alternatif terhadap statistik Wald untuk inferensi koefisien.

Detail Estimasi dan Inferensi Regresi Poisson

Regresi Poisson digunakan untuk data hitungan. Asumsi utama adalah: - Variansi sama dengan rata-rata (\(Var(Y) = \mu\))

Diagnostik Umum:

  • Deviance: ukuran goodness-of-fit
  • AIC: digunakan untuk membandingkan beberapa model
  • Overdispersion: jika \(Var(Y) > \mu\), gunakan model quasi-Poisson atau negatif binomial

Contoh Diagnostik di R:

summary(model_pois)$deviance
## [1] 98.65831
AIC(model_pois)
## [1] 373.3214

Residual Plot

plot(fitted(model_pois), residuals(model_pois, type = "pearson"),
     main = "Pearson Residuals vs Fitted",
     xlab = "Fitted values", ylab = "Pearson residuals")
abline(h = 0, col = "red")

Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio

Dalam regresi logistik, prediktor tidak terbatas hanya pada variabel numerik kontinu. Variabel kategorik seperti nominal (misalnya jenis kelamin), ordinal (misalnya tingkat konsumsi garam), maupun rasio (misalnya usia) dapat digunakan sebagai prediktor. Perlakuan terhadap tipe data ini akan mempengaruhi interpretasi dan hasil model.

Simulasi Data

Misalkan kita ingin memodelkan kemungkinan seseorang mengalami tekanan darah tinggi berdasarkan beberapa faktor: - Jenis Kelamin (nominal: pria, wanita) - Tingkat Konsumsi Garam (ordinal: rendah, sedang, tinggi) - Usia (rasio: dalam tahun)

set.seed(123)
n <- 300
usia <- round(rnorm(n, mean = 45, sd = 12))
jenis_kelamin <- sample(c("Pria", "Wanita"), n, replace = TRUE)
konsumsi_garam <- factor(sample(c("Rendah", "Sedang", "Tinggi"), n, replace = TRUE),
                         levels = c("Rendah", "Sedang", "Tinggi"), ordered = TRUE)

# Asumsikan probabilitas hipertensi dipengaruhi oleh usia dan tingkat konsumsi garam
logit_p <- -4 + 0.05 * usia + as.numeric(konsumsi_garam) * 0.8
prob <- 1 / (1 + exp(-logit_p))
hipertensi <- rbinom(n, 1, prob)

data_kesehatan <- data.frame(hipertensi, usia, jenis_kelamin, konsumsi_garam)
head(data_kesehatan)
##   hipertensi usia jenis_kelamin konsumsi_garam
## 1          0   38        Wanita         Rendah
## 2          0   42        Wanita         Rendah
## 3          1   64        Wanita         Sedang
## 4          0   46          Pria         Tinggi
## 5          0   47          Pria         Tinggi
## 6          1   66          Pria         Sedang

Eksplorasi Data

table(data_kesehatan$hipertensi)
## 
##   0   1 
## 159 141
table(data_kesehatan$jenis_kelamin, data_kesehatan$hipertensi)
##         
##           0  1
##   Pria   85 69
##   Wanita 74 72
table(data_kesehatan$konsumsi_garam, data_kesehatan$hipertensi)
##         
##           0  1
##   Rendah 62 30
##   Sedang 60 40
##   Tinggi 37 71
summary(data_kesehatan)
##    hipertensi        usia       jenis_kelamin      konsumsi_garam
##  Min.   :0.00   Min.   :17.00   Length:300         Rendah: 92    
##  1st Qu.:0.00   1st Qu.:38.00   Class :character   Sedang:100    
##  Median :0.00   Median :44.00   Mode  :character   Tinggi:108    
##  Mean   :0.47   Mean   :45.41                                    
##  3rd Qu.:1.00   3rd Qu.:53.00                                    
##  Max.   :1.00   Max.   :84.00

Perlakuan Variabel Ordinal

Treat sebagai Nominal (Dummy Variable)

Variabel ordinal seperti konsumsi_garam bisa dianggap sebagai nominal, yaitu menggunakan dummy variable untuk setiap kategori.

model_nominal <- glm(hipertensi ~ usia + jenis_kelamin + konsumsi_garam, 
                     data = data_kesehatan, family = binomial)
summary(model_nominal)
## 
## Call:
## glm(formula = hipertensi ~ usia + jenis_kelamin + konsumsi_garam, 
##     family = binomial, data = data_kesehatan)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -3.18553    0.59696  -5.336 9.49e-08 ***
## usia                 0.06533    0.01264   5.170 2.34e-07 ***
## jenis_kelaminWanita  0.11802    0.25784   0.458    0.647    
## konsumsi_garam.L     1.11791    0.22980   4.865 1.15e-06 ***
## konsumsi_garam.Q     0.33840    0.22120   1.530    0.126    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 414.81  on 299  degrees of freedom
## Residual deviance: 357.76  on 295  degrees of freedom
## AIC: 367.76
## 
## Number of Fisher Scoring iterations: 4

Interpretasi:

  1. Intercept: -3.18553 Ini adalah log-odds hipertensi ketika semua prediktor = 0, yaitu: usia = 0, jenis_kelamin = “Pria”, dan konsumsi_garam = “Rendah”.

    Nilainya sendiri tidak terlalu penting secara praktis karena usia = 0 biasanya tidak relevan, tapi secara statistik ini adalah titik awal model.

  2. usia: 0.06533 (p < 0.001) Setiap peningkatan 1 tahun usia dikaitkan dengan peningkatan log-odds hipertensi sebesar 0.065.

    Dalam bentuk odds ratio: exp(0.06533) ≈ 1.067, artinya:

    Setiap kenaikan 1 tahun usia meningkatkan peluang terkena hipertensi sekitar 6.7%, jika faktor lain konstan.

  3. jenis_kelaminWanita: 0.11802 (p = 0.647) Wanita memiliki log-odds hipertensi sedikit lebih tinggi dari pria, tapi perbedaannya tidak signifikan secara statistik (p > 0.05).

    Odds ratio: exp(0.11802) ≈ 1.125, artinya wanita 12.5% lebih mungkin terkena hipertensi dibanding pria, tapi hasil ini tidak bisa disimpulkan secara kuat.

  4. konsumsi_garam.L: 1.11791 (p < 0.001)

  5. konsumsi_garam.Q: 0.33840 (p = 0.126) Di sini digunakan kontrasku orthogonal polynomial coding (.L, .Q) untuk faktor ordinal:

    .L = linear trend (rendah → sedang → tinggi)

    .Q = quadratic trend

    Interpretasi konsumsi_garam.L (linear):

    Efek linier signifikan dan positif.

    exp(1.11791) ≈ 3.06, artinya:

    Semakin tinggi konsumsi garam secara linear, peluang hipertensi sekitar 3 kali lebih besar dibandingkan konsumsi rendah, dengan asumsi linearitas.

    Interpretasi konsumsi_garam.Q (quadratic):

    Tidak signifikan (p = 0.126), menunjukkan bahwa tidak ada pola kuadratik (misal: sedang lebih berisiko dari rendah & tinggi) yang signifikan.

    Kecocokan Model Null deviance: 414.81, Residual deviance: 357.76 → Model menjelaskan sebagian variabilitas.

    AIC: 367.76 → Bisa digunakan untuk membandingkan model lain, lebih kecil = lebih baik.

    Kesimpulan Utama Usia dan konsumsi garam berpengaruh signifikan terhadap kemungkinan hipertensi.

    Jenis kelamin tidak signifikan.

    Efek konsumsi garam cenderung linier: semakin tinggi, semakin besar peluang hipertensi.

    Interpretasi yang kamu tulis di akhir:

    “Setiap kategori konsumsi_garam dibandingkan terhadap kategori referensi (Rendah)” Itu tidak sepenuhnya tepat, karena model menggunakan kontras polinomial, bukan dummy coding. Jadi tidak dibandingkan langsung per kategori, tapi lebih melihat pola tren (linear dan kuadratik).

Treat sebagai Rasio (Numeric Rank)

Jika ingin memperlakukan konsumsi_garam sebagai ordinal numerik, kita konversi menjadi urutan numerik.

data_kesehatan$garam_rank <- as.numeric(data_kesehatan$konsumsi_garam)
model_ordinal <- glm(hipertensi ~ usia + jenis_kelamin + garam_rank, 
                     data = data_kesehatan, family = binomial)
summary(model_ordinal)
## 
## Call:
## glm(formula = hipertensi ~ usia + jenis_kelamin + garam_rank, 
##     family = binomial, data = data_kesehatan)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -4.75817    0.73560  -6.468 9.90e-11 ***
## usia                 0.06487    0.01256   5.165 2.40e-07 ***
## jenis_kelaminWanita  0.06521    0.25433   0.256    0.798    
## garam_rank           0.80502    0.16302   4.938 7.89e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 414.81  on 299  degrees of freedom
## Residual deviance: 360.12  on 296  degrees of freedom
## AIC: 368.12
## 
## Number of Fisher Scoring iterations: 3

Interpretasi:

  • (Intercept): 4.75817 Log-odds hipertensi untuk referensi: usia = 0, jenis kelamin = pria, garam_rank = paling rendah.

    Karena usia = 0 tidak realistis, intercept ini hanya berfungsi sebagai titik awal model, bukan interpretasi praktis.

  • usia: 0.06487 (p < 0.001) Setiap penambahan 1 tahun usia meningkatkan log-odds terkena hipertensi sebesar 0.06487.

    Dalam bentuk odds ratio: exp(0.06487) ≈ 1.067 Artinya: Setiap 1 tahun bertambah, peluang hipertensi meningkat ~6.7%.

  • jenis_kelaminWanita: 0.06521 (p = 0.798) Wanita memiliki peluang hipertensi sedikit lebih besar dibanding pria, tapi tidak signifikan secara statistik.

    exp(0.06521) ≈ 1.067 → ≈6.7% lebih besar, tapi tidak cukup bukti untuk menyimpulkan perbedaan gender.

  • garam_rank: 0.80502 (p < 0.001) Setiap kenaikan satu tingkat dalam peringkat konsumsi garam dikaitkan dengan peningkatan log-odds hipertensi sebesar 0.805.

    exp(0.80502) ≈ 2.24 Artinya:

  • Setiap naik satu tingkat konsumsi garam (misal: dari rendah → sedang atau sedang → tinggi), peluang hipertensi meningkat ~2.24 kali lipat, signifikan secara statistik.

  • Kecocokan Model Null deviance: 414.81 → tanpa prediktor

    Residual deviance: 360.12 → dengan prediktor

    Penurunan deviance signifikan, menunjukkan model lebih baik daripada hanya intercept.

    AIC: 368.12 → Untuk membandingkan antar model; lebih kecil lebih baik.

  • Kesimpulan usia dan konsumsi garam (garam_rank) adalah prediktor yang signifikan terhadap hipertensi.

    Jenis kelamin tidak signifikan dalam model ini.

    Efek garam_rank linier: semakin tinggi peringkat konsumsi garam, semakin besar peluang hipertensi.

    Model cukup baik menjelaskan variasi data (devian residual menurun, AIC cukup rendah).

Kesimpulan

Dalam regresi logistik, variabel kategorik dapat diperlakukan berbeda tergantung konteks analisis: - Perlakuan sebagai dummy (nominal) cocok jika tidak ada urutan alami. - Perlakuan sebagai ordinal (numeric rank) cocok jika ada urutan logis antar kategori.

Pemilihan perlakuan yang tepat terhadap prediktor kategorik sangat penting untuk memperoleh model yang bermakna dan interpretable.

Pemilihan Model Regresi Logistik dan Evaluasi

Membangun Model Regresi Logistik: Pendekatan Confirmatory dan Explanatory

  • Confirmatory modeling: Pendekatan ini dilakukan berdasarkan teori atau temuan terdahulu. Model dibangun secara spesifik untuk menguji hipotesis yang telah ditentukan sebelumnya.
  • Explanatory modeling: Bersifat eksploratif, bertujuan untuk menemukan pola atau hubungan terbaik antar variabel dalam data, tanpa asumsi teoritis awal yang kuat.

Metode Stepwise: Forward, Backward, dan Kedua Arah

Langkah seleksi stepwise membantu menentukan kombinasi prediktor terbaik untuk model regresi logistik. Kriteria seperti AIC digunakan untuk menilai kualitas model.

  • Forward selection: menambahkan variabel satu per satu dimulai dari model kosong.
  • Backward elimination: memulai dari model penuh lalu menghapus variabel yang tidak signifikan.
  • Stepwise (dua arah): gabungan keduanya, menambah dan menghapus variabel secara iteratif.
# Model penuh
model_full <- glm(hipertensi ~ usia + jenis_kelamin + konsumsi_garam,
                  data = data_kesehatan, family = binomial)

# Stepwise dua arah
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
model_step <- stepAIC(model_full, direction = "both")
## Start:  AIC=367.76
## hipertensi ~ usia + jenis_kelamin + konsumsi_garam
## 
##                  Df Deviance    AIC
## - jenis_kelamin   1   357.97 365.97
## <none>                357.76 367.76
## - konsumsi_garam  2   386.74 392.74
## - usia            1   389.09 397.09
## 
## Step:  AIC=365.97
## hipertensi ~ usia + konsumsi_garam
## 
##                  Df Deviance    AIC
## <none>                357.97 365.97
## + jenis_kelamin   1   357.76 367.76
## - konsumsi_garam  2   387.07 391.07
## - usia            1   389.61 395.61
summary(model_step)
## 
## Call:
## glm(formula = hipertensi ~ usia + konsumsi_garam, family = binomial, 
##     data = data_kesehatan)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -3.13827    0.58694  -5.347 8.95e-08 ***
## usia              0.06555    0.01263   5.191 2.09e-07 ***
## konsumsi_garam.L  1.12443    0.22941   4.901 9.52e-07 ***
## konsumsi_garam.Q  0.32491    0.21906   1.483    0.138    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 414.81  on 299  degrees of freedom
## Residual deviance: 357.97  on 296  degrees of freedom
## AIC: 365.97
## 
## Number of Fisher Scoring iterations: 4
  • Forward selection: Mulai dari model kosong, prediktor ditambahkan satu per satu berdasarkan peningkatan kriteria seperti AIC.
  • Backward elimination: Mulai dari model penuh, prediktor dieliminasi satu per satu.
  • Stepwise (kedua arah): Gabungan forward dan backward, secara iteratif menambahkan dan menghapus prediktor.

Evaluasi Model: ROC dan AUC

ROC dan AUC digunakan untuk mengevaluasi kemampuan model dalam mengklasifikasikan data biner dengan berbagai ambang batas prediksi.

  • ROC (Receiver Operating Characteristic): Kurva yang membandingkan True Positive Rate (TPR) vs False Positive Rate (FPR) pada berbagai threshold.
  • AUC (Area Under Curve): Luas di bawah kurva ROC. Nilai mendekati 1 menunjukkan performa klasifikasi yang sangat baik.
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
roc_obj <- roc(data_kesehatan$hipertensi, fitted(model_step))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj, main = "ROC Curve")

auc(roc_obj)
## Area under the curve: 0.74
  • ROC (Receiver Operating Characteristic): Grafik yang menunjukkan trade-off antara True Positive Rate (TPR) dan False Positive Rate (FPR).
  • AUC (Area Under Curve): Luas di bawah kurva ROC. Nilai AUC yang mendekati 1 menunjukkan model klasifikasi yang sangat baik.

Pseudo R-Squared

Karena regresi logistik tidak memiliki R² seperti regresi linear, maka digunakan alternatif pseudo R-squared, seperti:

  • Cox & Snell R²: dibatasi karena tidak bisa mencapai maksimum 1.
  • Nagelkerke R²: versi penyesuaian dari Cox & Snell agar bisa mencapai 1.

Nilainya berkisar antara 0 dan 1; semakin tinggi, semakin baik kesesuaian model.

Tabel Klasifikasi dan Evaluasi

Tabel klasifikasi membandingkan hasil prediksi model terhadap kenyataan. Ini memungkinkan kita menghitung ukuran kinerja seperti akurasi, sensitivitas, dan spesifisitas.

pred <- ifelse(fitted(model_step) > 0.5, 1, 0)
table(Prediksi = pred, Aktual = data_kesehatan$hipertensi)
##         Aktual
## Prediksi   0   1
##        0 118  57
##        1  41  84

Tabel ini membandingkan prediksi model terhadap kenyataan. Dari tabel ini bisa dihitung:

  • Akurasi: Proposi prediksi yang benar.
  • True Positive, False Positive, True Negative, False Negative.

Metode Perbandingan Model dalam Regresi Logistik

Untuk membandingkan dua model: - Gunakan AIC sebagai indikator kompleksitas dan kualitas model. - Gunakan likelihood-ratio test untuk menguji apakah model baru secara signifikan lebih baik dari model lama.

Likelihood-Ratio Test

Uji ini membandingkan dua model terstruktur, biasanya model penuh dengan model sederhana: \[ G^2 = -2 (\log L_0 - \log L_1) \] G² mengikuti distribusi chi-square.

Prinsip Parsimony

Model yang paling sederhana dengan akurasi terbaik lebih disukai. Jangan menambahkan variabel yang tidak signifikan atau tidak berguna secara praktis.

Evaluasi Tabel Klasifikasi dan Akurasi Model

Akurasi dihitung dari proporsi klasifikasi benar terhadap total: \[ \text{Akurasi} = \frac{TP + TN}{TP + TN + FP + FN} \]

Sensitivitas dan Spesifisitas

  • Sensitivitas: Kemampuan model mendeteksi kasus positif dengan benar (\(TP / (TP + FN)\))
  • Spesifisitas: Kemampuan model mengidentifikasi kasus negatif dengan benar (\(TN / (TN + FP)\))

Detail ROC, Penjelasan Kurva ROC (Receiver Operating Characteristic)

  • Titik-titik pada kurva ROC menunjukkan kinerja model pada berbagai threshold.
  • Kurva yang membentuk area luas mendekati sudut kiri atas menunjukkan performa prediksi yang sangat baik.

Precision-Recall Curve (PR Curve)

  • Precision: Proporsi prediksi positif yang benar-benar positif.
  • Recall: Sinonim dari sensitivitas.
  • PR Curve berguna untuk dataset tidak seimbang (imbalance).

Pseudo R-Squares pada Regresi Logistik

Digunakan untuk mengukur goodness-of-fit secara global:

  • Cox & Snell R²: Tidak pernah mencapai 1.
  • Nagelkerke R²: Koreksi Cox & Snell agar bisa mencapai maksimum 1.

Kesimpulan

Evaluasi regresi logistik harus melibatkan aspek statistik, prediktif, dan interpretatif. Penggunaan ROC, AUC, pseudo R², serta tabel klasifikasi dan uji perbandingan model memungkinkan analisis yang menyeluruh untuk memilih model terbaik.

Apa itu Distribusi Multinomial

Distribusi multinomial merupakan perluasan dari distribusi binomial untuk situasi di mana variabel respon memiliki lebih dari dua kategori. Distribusi ini cocok digunakan ketika kejadian memiliki beberapa kemungkinan hasil yang saling eksklusif dan mencakup semua kemungkinan.

Misalnya, jika seseorang dapat dikategorikan ke dalam salah satu dari tiga tingkat tekanan darah: Normal, Prehipertensi, atau Hipertensi, maka distribusi multinomial cocok digunakan untuk memodelkan probabilitas keanggotaan dalam masing-masing kategori tersebut.

Studi Kasus

Dalam bidang kesehatan, kita ingin memodelkan kategori tekanan darah seseorang (Normal, Prehipertensi, atau Hipertensi) berdasarkan:

  • Usia (variabel rasio)
  • Tingkat Konsumsi Garam (variabel ordinal: Rendah, Sedang, Tinggi)

Tujuannya adalah mengetahui apakah usia dan konsumsi garam berpengaruh terhadap kemungkinan seseorang termasuk dalam kategori tekanan darah tertentu.

Multinomial Logistic Regression

Baseline-category logit model

Model ini memilih salah satu kategori sebagai referensi (misalnya “Normal”), lalu membandingkan log odds setiap kategori lainnya terhadap kategori referensi.

\[ \log\left( \frac{P(Y = k)}{P(Y = \text{Normal})} \right) = \beta_{0k} + \beta_{1k}X_1 + ... + \beta_{pk}X_p \]

Setiap kategori (selain referensi) memiliki satu set koefisien regresi.

Estimasi Parameter

Estimasi dilakukan dengan metode Maximum Likelihood Estimation (MLE). Proses ini mengoptimalkan fungsi likelihood untuk mendapatkan parameter terbaik yang menjelaskan data.

Simulasi Data

set.seed(123)
n <- 400
usia <- round(rnorm(n, mean = 50, sd = 12))
konsumsi_garam <- factor(sample(c("Rendah", "Sedang", "Tinggi"), n, replace = TRUE),
                         levels = c("Rendah", "Sedang", "Tinggi"), ordered = TRUE)
logit1 <- -4 + 0.04 * usia + 0.8 * as.numeric(konsumsi_garam)     # Prehipertensi
logit2 <- -6 + 0.08 * usia + 1.3 * as.numeric(konsumsi_garam)     # Hipertensi
p1 <- exp(logit1)
p2 <- exp(logit2)
p0 <- 1
total <- p0 + p1 + p2
prob <- cbind(p0/total, p1/total, p2/total)
colnames(prob) <- c("Normal", "Prehipertensi", "Hipertensi")
tekanan <- apply(prob, 1, function(p) sample(colnames(prob), 1, prob = p))

data_multi <- data.frame(tekanan = factor(tekanan, levels = colnames(prob)), usia, konsumsi_garam)

Estimasi Model

library(nnet)
model_multi <- multinom(tekanan ~ usia + konsumsi_garam, data = data_multi)
## # weights:  15 (8 variable)
## initial  value 439.444915 
## iter  10 value 339.792980
## final  value 339.182752 
## converged
summary(model_multi)
## Call:
## multinom(formula = tekanan ~ usia + konsumsi_garam, data = data_multi)
## 
## Coefficients:
##               (Intercept)       usia konsumsi_garam.L konsumsi_garam.Q
## Prehipertensi   -2.303937 0.04001155         1.778844       0.06823459
## Hipertensi      -3.389385 0.08223902         2.121761       0.51925352
## 
## Std. Errors:
##               (Intercept)       usia konsumsi_garam.L konsumsi_garam.Q
## Prehipertensi   0.7620438 0.01567856        0.3252661        0.2766978
## Hipertensi      0.6700599 0.01359363        0.2726192        0.2343406
## 
## Residual Deviance: 678.3655 
## AIC: 694.3655

Nilai P-Value dan Interpretasi

z_val <- summary(model_multi)$coefficients / summary(model_multi)$standard.errors
p_val <- 2 * (1 - pnorm(abs(z_val)))
p_val
##                (Intercept)         usia konsumsi_garam.L konsumsi_garam.Q
## Prehipertensi 2.499802e-03 1.071091e-02     4.528677e-08       0.80521528
## Hipertensi    4.229415e-07 1.450095e-09     7.105427e-15       0.02670471

Interpretasi:

  • Untuk setiap kategori (Prehipertensi dan Hipertensi), kita memiliki satu set koefisien terhadap referensi “Normal”.
  • Jika p-value < 0.05, maka prediktor dianggap berpengaruh signifikan terhadap kategori tersebut.

Contoh:

  • Jika koefisien usia untuk kategori Hipertensi signifikan dan positif, maka semakin tua usia, semakin besar kemungkinan masuk ke kategori Hipertensi dibanding Normal.
  • Jika konsumsi_garamTinggi signifikan dan positif pada kategori Prehipertensi, artinya orang yang konsumsi garam tinggi lebih berisiko mengalami Prehipertensi dibanding Normal.

Prediksi dan Validasi

prediksi <- predict(model_multi)
table(Predicted = prediksi, Actual = data_multi$tekanan)
##                Actual
## Predicted       Normal Prehipertensi Hipertensi
##   Normal            76            19         37
##   Prehipertensi      0             0          0
##   Hipertensi        47            49        172
mean(prediksi == data_multi$tekanan)  # Akurasi prediksi
## [1] 0.62

Interpretasi: - Tabel konfusi menunjukkan bagaimana prediksi model sesuai dengan kenyataan. - Akurasi dihitung sebagai proporsi prediksi yang benar terhadap total observasi. - Jika akurasi tinggi dan proporsi antar kategori seimbang, model bisa dianggap cukup baik.

Kesimpulan

Regresi logistik multinomial adalah metode yang tepat saat variabel respon memiliki lebih dari dua kategori. Dalam contoh ini, model berhasil mengaitkan usia dan konsumsi garam terhadap kategori tekanan darah. Hasil analisis menunjukkan bahwa:

  • Usia berperan penting meningkatkan risiko dari Normal → Prehipertensi → Hipertensi.
  • Tingkat konsumsi garam secara signifikan meningkatkan kemungkinan masuk ke kategori tekanan darah lebih tinggi.

Evaluasi seperti akurasi dan signifikansi parameter memberikan gambaran apakah model bisa digunakan untuk prediksi maupun pengambilan keputusan kebijakan kesehatan.

Regresi Logistik Ordinal

Konsep Cumulative Logit Model

Regresi logistik ordinal digunakan ketika variabel respon bersifat kategorik dan memiliki urutan (misalnya: rendah, sedang, tinggi). Model yang umum digunakan adalah cumulative logit model atau proportional odds model.

Cumulative logit model menghitung probabilitas kumulatif respon berada pada atau di bawah level tertentu: \[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j - \beta^T x \]

Setiap batas kategori memiliki intercept berbeda \(\alpha_j\), tetapi koefisien regresor \(\beta\) sama untuk semua.

Interpretasi Koefisien

Koefisien \(\beta\) menunjukkan pengaruh variabel terhadap kemungkinan berada pada kategori yang lebih tinggi.

  • Koefisien positif: meningkatkan peluang respon berada pada kategori lebih rendah.
  • Koefisien negatif: meningkatkan peluang berada pada kategori lebih tinggi.

Contoh Data

Kita gunakan data simulasi tekanan darah dengan tiga kategori ordinal: Normal, Prehipertensi, dan Hipertensi.

library(MASS)
set.seed(123)
n <- 300
usia <- round(rnorm(n, 50, 10))
konsumsi_garam <- factor(sample(c("Rendah", "Sedang", "Tinggi"), n, replace = TRUE),
                         levels = c("Rendah", "Sedang", "Tinggi"), ordered = TRUE)

# Probabilitas logit kumulatif
logit1 <- -1 + 0.04 * usia + 0.8 * as.numeric(konsumsi_garam)
logit2 <- 2 + 0.04 * usia + 0.8 * as.numeric(konsumsi_garam)
p1 <- 1 / (1 + exp(-logit1))
p2 <- 1 / (1 + exp(-logit2))

kategori <- mapply(function(p1, p2) {
  r <- runif(1)
  if (r < p1) "Normal" else if (r < p2) "Prehipertensi" else "Hipertensi"
}, p1, p2)

tekanan_ord <- factor(kategori, levels = c("Normal", "Prehipertensi", "Hipertensi"), ordered = TRUE)
data_ord <- data.frame(tekanan_ord, usia, konsumsi_garam)

Estimasi Model Ordinal

model_ord <- polr(tekanan_ord ~ usia + konsumsi_garam, data = data_ord, Hess = TRUE)
summary(model_ord)
## Call:
## polr(formula = tekanan_ord ~ usia + konsumsi_garam, data = data_ord, 
##     Hess = TRUE)
## 
## Coefficients:
##                     Value Std. Error t value
## usia             -0.01942    0.02877 -0.6752
## konsumsi_garam.L -0.80213    0.43137 -1.8595
## konsumsi_garam.Q  0.90033    0.63554  1.4166
## 
## Intercepts:
##                          Value   Std. Error t value
## Normal|Prehipertensi      2.1213  1.4409     1.4722
## Prehipertensi|Hipertensi  4.2768  1.5819     2.7035
## 
## Residual Deviance: 129.1768 
## AIC: 139.1768

Nilai P-Value

ctable <- coef(summary(model_ord))
p_val <- 2 * (1 - pnorm(abs(ctable[, "t value"])))
cbind(ctable, "p value" = p_val)
##                                Value Std. Error   t value    p value
## usia                     -0.01942403 0.02876948 -0.675161 0.49957350
## konsumsi_garam.L         -0.80212644 0.43136627 -1.859502 0.06295602
## konsumsi_garam.Q          0.90033257 0.63554145  1.416639 0.15658861
## Normal|Prehipertensi      2.12126585 1.44092726  1.472153 0.14097948
## Prehipertensi|Hipertensi  4.27679831 1.58191993  2.703549 0.00686033

Interpretasi:

  • p-value < 0.05 menunjukkan bahwa variabel signifikan terhadap peningkatan level tekanan darah.

Prediksi Probabilitas

predict(model_ord, newdata = data.frame(usia = 60, konsumsi_garam = "Tinggi"), type = "probs")
##        Normal Prehipertensi    Hipertensi 
##    0.97029652    0.02616978    0.00353370

Interpretasi:

  • Probabilitas untuk setiap level tekanan darah dapat dihitung untuk kombinasi prediktor tertentu.

Goodness-of-Fit dan Proportional Odds

  • Model ini mengasumsikan bahwa pengaruh prediktor konstan di semua level kategori (proportional odds).
  • Jika asumsi tidak terpenuhi, bisa gunakan model alternatif.

Alternatif Model Ordinal

  • Partial proportional odds model: memungkinkan koefisien bervariasi antar cutpoints
  • Continuation ratio model atau adjacent category model untuk kasus khusus

Kesimpulan

Model regresi logistik ordinal sangat cocok untuk data dengan respon terurut. Cumulative logit model menyediakan pendekatan efisien dengan interpretasi jelas, namun bergantung pada asumsi proportional odds.

Asumsi Paralelisme dalam Regresi Logistik Ordinal

Asumsi ini menyatakan bahwa efek prediktor sama pada setiap transisi antar kategori. Jika asumsi ini dilanggar, hasil estimasi bisa bias. Uji Brant atau perbandingan model nested bisa digunakan untuk mengevaluasi asumsi ini.

Log Linear Model

Tabel Kontingensi dan Model Log Linear

Model log-linear digunakan untuk menganalisis hubungan antara dua atau lebih variabel kategorikal yang disusun dalam tabel kontingensi. Model ini tidak memisahkan antara variabel respon dan prediktor, melainkan memperlakukan semua variabel secara simetris.

Tujuan utama dari model log-linear adalah untuk menjelaskan frekuensi dalam sel-sel tabel kontingensi berdasarkan efek dari masing-masing variabel dan interaksi di antara mereka.

Model Saturated

Model saturated adalah model yang mencakup semua kemungkinan efek utama dan interaksi antara variabel. Model ini akan memiliki fit sempurna karena menyertakan semua parameter yang mungkin.

Contoh untuk tabel 2x2x2 (A, B, C)

\[\log(\mu{ijk}) = \lambda + \lambda^A_i + \lambda^B_j + \lambda^C_k + \lambda^{AB}{ij} + \lambda^{AC}{ik} + \lambda^{BC}{jk} + \lambda^{ABC}_{ijk}\]

Model Independent

Model ini mengasumsikan bahwa variabel-variabel tidak saling berinteraksi. Hanya efek utama yang disertakan.

\[ \log(\mu_{ijk}) = \lambda + \lambda^A_i + \lambda^B_j + \lambda^C_k \]

Jika model ini cocok dengan baik, maka tidak ada asosiasi atau ketergantungan antara ketiga variabel.

Odds Ratio dan Interpretasi

Odds ratio dalam model log-linear dapat digunakan untuk mengukur kekuatan asosiasi antara dua variabel. Nilai odds ratio = 1 menunjukkan tidak ada asosiasi.

Interpretasi:

  • OR > 1: hubungan positif
  • OR < 1: hubungan negatif

Estimasi Parameter

Estimasi parameter dilakukan dengan metode maximum likelihood, dan biasanya dihitung menggunakan fungsi loglm() dari paket MASS.

library(MASS)
data(HairEyeColor)
tab <- margin.table(HairEyeColor, c(1, 2))  # Tabel 2D
model_log <- loglm(~ Hair + Eye, data = tab)
summary(model_log)
## Formula:
## ~Hair + Eye
## attr(,"variables")
## list(Hair, Eye)
## attr(,"factors")
##      Hair Eye
## Hair    1   0
## Eye     0   1
## attr(,"term.labels")
## [1] "Hair" "Eye" 
## 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 146.4436  9        0
## Pearson          138.2898  9        0

Model Lebih Sederhana dan Perbandingan Model

Model nested (satu model di dalam model lain) dapat dibandingkan menggunakan likelihood ratio test (G²). Ini berguna untuk memilih model yang lebih sederhana tetapi tetap cukup baik.

model_full <- loglm(~ Hair * Eye, data = tab)
model_reduced <- loglm(~ Hair + Eye, data = tab)
anova(model_reduced, model_full)
## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Hair + Eye 
## Model 2:
##  ~Hair * Eye 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   146.4436  9                                    
## Model 2     0.0000  0   146.4436         9              0
## Saturated   0.0000  0     0.0000         0              1

Studi Kasus

Misalkan kita memiliki data kebiasaan merokok (Rokok: Ya/Tidak), status hipertensi (Hipertensi: Ya/Tidak), dan konsumsi garam (Tinggi/Rendah). Kita dapat membangun model log-linear untuk melihat apakah terdapat interaksi antara ketiga faktor tersebut.

data_case <- expand.grid(Rokok = c("Ya", "Tidak"),
                         Hipertensi = c("Ya", "Tidak"),
                         Garam = c("Tinggi", "Rendah"))
data_case$Freq <- c(30, 10, 25, 15, 20, 5, 10, 5)

model_case <- loglm(Freq ~ Rokok * Hipertensi * Garam, data = data_case)
summary(model_case)
## Formula:
## Freq ~ Rokok * Hipertensi * Garam
## attr(,"variables")
## list(Freq, Rokok, Hipertensi, Garam)
## attr(,"factors")
##            Rokok Hipertensi Garam Rokok:Hipertensi Rokok:Garam Hipertensi:Garam
## Freq           0          0     0                0           0                0
## Rokok          1          0     0                1           1                0
## Hipertensi     0          1     0                1           0                1
## Garam          0          0     1                0           1                1
##            Rokok:Hipertensi:Garam
## Freq                            0
## Rokok                           1
## Hipertensi                      1
## Garam                           1
## attr(,"term.labels")
## [1] "Rokok"                  "Hipertensi"             "Garam"                 
## [4] "Rokok:Hipertensi"       "Rokok:Garam"            "Hipertensi:Garam"      
## [7] "Rokok:Hipertensi:Garam"
## attr(,"order")
## [1] 1 1 1 2 2 2 3
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(Freq, Rokok, Hipertensi, Garam)
## attr(,"dataClasses")
##       Freq      Rokok Hipertensi      Garam 
##  "numeric"   "factor"   "factor"   "factor" 
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1

Interpretasi output akan menunjukkan apakah interaksi tiga arah signifikan, atau apakah model dapat disederhanakan.

Kesimpulan

Model log-linear sangat bermanfaat untuk mengidentifikasi struktur asosiasi dalam data kategorikal. Dengan membandingkan model saturated dan model yang lebih sederhana, kita bisa memahami hubungan antar variabel serta menentukan model yang paling efisien dan informatif.

Model Log-Linear

Analisis data kategorik merupakan bagian penting dalam statistika terapan karena banyak fenomena di dunia nyata yang menghasilkan data dalam bentuk kategori, seperti jenis kelamin, status pekerjaan, tingkat pendidikan, preferensi konsumen, atau hasil diagnosis medis. Data kategori ini umumnya dianalisis menggunakan tabel kontingensi, model log-linier, dan model regresi logistik. Masing-masing pendekatan memiliki kekuatan dan kelemahan tergantung pada tujuan analisis dan struktur data.

Tabel kontingensi digunakan sebagai langkah awal eksplorasi untuk melihat hubungan antara dua atau lebih variabel kategorik. Misalnya, dalam studi tentang efek obat terhadap serangan jantung, tabel kontingensi dapat menyajikan jumlah pasien yang mengalami atau tidak mengalami serangan jantung, berdasarkan jenis obat yang dikonsumsi. Tabel ini membantu mengidentifikasi pola awal dan menghitung ukuran asosiasi seperti odds ratio, risk ratio, dan chi-square statistic untuk menguji independensi antar variabel.

Namun, ketika ingin membangun model statistik yang dapat mengendalikan efek dari banyak variabel dan interaksinya secara simultan, maka model log-linier menjadi sangat berguna. Model log-linier adalah bentuk khusus dari Generalized Linear Model (GLM) yang digunakan pada frekuensi sel dalam tabel kontingensi dan mengasumsikan distribusi Poisson. Tidak seperti regresi logistik, model log-linier tidak menetapkan variabel mana yang dependen dan mana yang independen, karena seluruh variabel diperlakukan secara simetris. Model ini lebih cocok ketika tujuan analisis adalah untuk memahami struktur asosiasi atau independensi antar variabel, bukan untuk prediksi.

Struktur model log-linier ditentukan berdasarkan efek utama dari masing-masing variabel serta interaksi di antara variabel-variabel tersebut. Misalnya, dalam tabel kontingensi tiga arah (misalnya: jenis kelamin, status merokok, dan penyakit paru), model log-linier dapat membedakan apakah interaksi dua variabel cukup menjelaskan data, ataukah diperlukan interaksi tiga arah untuk menangkap struktur asosiasinya. Penyesuaian model dapat dilakukan menggunakan metode likelihood ratio test untuk membandingkan model sederhana dengan model lebih kompleks.

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

Dengan demikian, meskipun ketiganya beroperasi pada data kategorik, tabel kontingensi bersifat deskriptif, model log-linier bersifat eksploratif terhadap hubungan simetris, sedangkan regresi logistik bersifat prediktif terhadap outcome kategorik. Pemilihan metode tergantung pada apakah fokus utama analisis adalah deskripsi, eksplorasi struktur, atau prediksi hasil berdasarkan variabel penjelas. Kombinasi dari ketiganya sering digunakan dalam praktik untuk memperoleh pemahaman komprehensif dari data kategorik yang dianalisis.

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: Tabel kontingensi menyajikan jumlah frekuensi dari kombinasi kategori antar variabel.

Contoh tabel 2x2:

Sebuah studi dilakukan untuk mengevaluasi efektivitas Vaksin A dalam mencegah infeksi virus flu. Penelitian ini melibatkan dua kelompok:

Setelah periode observasi:

table_data <- matrix(c(15, 35, 45, 5), nrow=2,
       dimnames = list(Vaksin = c("Vaksin A", "Placebo"),
                       Infeksi = c("Ya", "Tidak")))
table_data
##           Infeksi
## Vaksin     Ya Tidak
##   Vaksin A 15    45
##   Placebo  35     5

Tabel kontingensi bersifat deskriptif dan tidak melibatkan pemodelan probabilitas.

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

\\\Rumus///

Contoh:

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

Model Regresi Logistik

Model regresi logistik biner:

\\\Rumus///

Contoh:

data_glm <- data.frame(
  Infeksi = c(1, 0, 1, 0),
  Vaksin = factor(c("Vaksin A", "Vaksin A", "Placebo", "Placebo")),
  Frek = c(15, 35, 45, 5)
)
model_logit <- glm(Infeksi ~ Vaksin, weights = Frek, family = binomial, data = data_glm)
summary(model_logit)
## 
## Call:
## glm(formula = Infeksi ~ Vaksin, family = binomial, data = data_glm, 
##     weights = Frek)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      2.1972     0.4714   4.661 3.15e-06 ***
## VaksinVaksin A  -3.0445     0.5634  -5.404 6.53e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 134.602  on 3  degrees of freedom
## Residual deviance:  93.595  on 2  degrees of freedom
## AIC: 97.595
## 
## Number of Fisher Scoring iterations: 5

Perbandingan Ketiga Pendekatan

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

Tabel Kontingensi dan Model LogLinear

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

table_data <- matrix(c(15, 35, 45, 5), nrow=2,
       dimnames = list(Vaksin = c("Vaksin A", "Placebo"),
                       Infeksi = c("Ya", "Tidak")))
table_data
##           Infeksi
## Vaksin     Ya Tidak
##   Vaksin A 15    45
##   Placebo  35     5

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

\\\Rumus///

Model Saturated

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

  • Cocok sempurna terhadap data

  • Tidak mengasumsikan independensi antar variabel

Contoh formulasi untuk tabel 2x2:

# Data efektivitas vaksin
datal <- matrix(c(15, 35, 45, 5), nrow=2, byrow=TRUE)
dimnames(datal) <- list(Vaksin = c("Vaksin A", "Placebo"),
                       Infeksi = c("Ya", "Tidak"))
ftable(datal)
##          Infeksi Ya Tidak
## Vaksin                   
## Vaksin A         15    35
## Placebo          45     5

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

model_saturated <- loglm(~ Vaksin * Infeksi, data = datal)
summary(model_saturated)
## Formula:
## ~Vaksin * Infeksi
## attr(,"variables")
## list(Vaksin, Infeksi)
## attr(,"factors")
##         Vaksin Infeksi Vaksin:Infeksi
## Vaksin       1       0              1
## Infeksi      0       1              1
## attr(,"term.labels")
## [1] "Vaksin"         "Infeksi"        "Vaksin:Infeksi"
## attr(,"order")
## [1] 1 1 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1

Model Independen

Model independen mengasumsikan bahwa tidak ada interaksi antara variabel:

\\\Rumus///

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

model_indep <- loglm(~ Vaksin + Infeksi, data = datal)
summary(model_indep)
## Formula:
## ~Vaksin + Infeksi
## attr(,"variables")
## list(Vaksin, Infeksi)
## attr(,"factors")
##         Vaksin Infeksi
## Vaksin       1       0
## Infeksi      0       1
## attr(,"term.labels")
## [1] "Vaksin"  "Infeksi"
## 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 41.00761  1 1.516379e-10
## Pearson          37.50000  1 9.141299e-10

Odds Ratio dan Interpretasi

Odds Ratio untuk tabel 2x2:

\\\Rumus///

Interpretasi nilai OR:

  • OR = 1: Tidak ada asosiasi
  • OR > 1: Asosiasi positif
  • OR < 1: Asosiasi negatif

Estimasi Parameter

Dalam model saturated:

  • Estimasi dilakukan dengan pembatasan seperti sum-to-zero

  • Estimasi parameter dilakukan dengan iterative proportional fitting (IPF)

logOR <- log((datal[1,1] * datal[2,2]) / (datal[1,2] * datal[2,1]))
logOR
## [1] -3.044522

Model lebih Sederhana dan Perbandingan Model

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

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

Studi Kasus: Hubungan Antara Kepuasan Pelayanan dan Loyalitas Pelanggan Sebuah Restoran

data_survey <- matrix(c(32, 190,
                        113, 611,
                        51, 326),
                      nrow = 3, byrow = TRUE,
                      dimnames = list(Kepuasan = c("Tidak Puas", "Cukup Puas", "Sangat Puas"),
                                      Loyalitas = c("Tidak Akan Kembali", "Akan Kembali")))

ftable(data_survey)
##             Loyalitas Tidak Akan Kembali Akan Kembali
## Kepuasan                                             
## Tidak Puas                            32          190
## Cukup Puas                           113          611
## Sangat Puas                           51          326
loglm(~ Kepuasan + Loyalitas, data = data_survey)
## Call:
## loglm(formula = ~Kepuasan + Loyalitas, data = data_survey)
## 
## Statistics:
##                        X^2 df  P(> X^2)
## Likelihood Ratio 0.8911136  2 0.6404675
## Pearson          0.8836760  2 0.6428538

Rangkuman Model Log Linear 2 Arah

Model LogLinear Pada Tabel kontingensi

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

\\\Rumus///

Perbedaan Utaman Antara Model LogLinear 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).

Estimasi Parameter LogLinear 2 Arah

Sistem Persamaan Model LogLinear

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

Constraint Sum-to-Zero

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

Rumus Estimasi Parameter dengan Sum-to-Zero Constraint

\[ \begin{aligned} \lambda^A_1 &= \frac{1}{2} \left[ \log \mu_{11} + \log \mu_{12} - (\log \mu_{21} + \log \mu_{22}) \right] \\ \lambda^B_1 &= \frac{1}{2} \left[ \log \mu_{11} + \log \mu_{21} - (\log \mu_{12} + \log \mu_{22}) \right] \\ \lambda^{AB}_{12} &= \frac{1}{4} \left[ \log \mu_{12} - \log \mu_{11} - \log \mu_{22} + \log \mu_{21} \right] \end{aligned} \]

Analisis Data Tabel Kontingensi 2x2

Diberikan data:

Vaksin Terinfeksi Tidak Terinfeksi
Vaksin A 15 35
Placebo 45 5

Bentuk Model LogLinear

Model LogLinear pada tabel 2x2:

dengan constraint sum-to-zero:

Estimasi Parameter Model (Manual, Sum-to-zero)

Misalkan:

  • A1 = Vaksin A, A2 = Placebo
  • B1 = Terinfeksi, B2 = Tidak Terinfeksi

Observasi:

  • \(n_{11} = 15, \quad n_{12} = 35\)
  • \(n_{21} = 45, \quad n_{22} = 5\)

Model log-linear:

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

Dengan constraint sum-to-zero:

\[ \lambda^A_1 + \lambda^A_2 = 0, \quad \lambda^B_1 + \lambda^B_2 = 0, \quad \sum_{i=1}^{2} \sum_{j=1}^{2} \lambda^{AB}_{ij} = 0 \]


Langkah 1: Rata-rata log frekuensi sel

\[ \lambda = \frac{1}{4} \left[ \log(15) + \log(35) + \log(45) + \log(5) \right] = \frac{1}{4} (2.7081 + 3.5553 + 3.8067 + 1.6094) = \frac{1}{4} (11.6795) = 2.9199 \]


Langkah 2: Efek utama A (Jenis Vaksin)

\[ \lambda^A_1 = \frac{1}{2} \left[ (\log 15 + \log 35) - (\log 45 + \log 5) \right] \\ = \frac{1}{2} (2.7081 + 3.5553 - 3.8067 - 1.6094) \\ = \frac{1}{2} (6.2634 - 5.4161) = \frac{1}{2} (0.8473) = 0.4237 \]

\[ \lambda^A_2 = -0.4237 \]


Langkah 3: Efek utama B (Infeksi)

\[ \lambda^B_1 = \frac{1}{2} \left[ (\log 15 + \log 45) - (\log 35 + \log 5) \right] \\ = \frac{1}{2} (2.7081 + 3.8067 - 3.5553 - 1.6094) \\ = \frac{1}{2} (6.5148 - 5.1647) = \frac{1}{2} (1.3501) = 0.6750 \]

\[ \lambda^B_2 = -0.6750 \]


Langkah 4: Efek interaksi

\[ \lambda^{AB}_{11} = \frac{1}{4} \left[ \log(15) - \log(35) - \log(45) + \log(5) \right] \\ = \frac{1}{4} (2.7081 - 3.5553 - 3.8067 + 1.6094) = \frac{1}{4} (-3.0445) = -0.7611 \]

\[ \lambda^{AB}_{12} = -\lambda^{AB}_{11} = 0.7611 \\ \lambda^{AB}_{21} = 0.7611 \\ \lambda^{AB}_{22} = -0.7611 \]


Ringkasan Parameter

  • \(\lambda = 2.9199\)
  • \(\lambda^A_1 = 0.4237, \quad \lambda^A_2 = -0.4237\)
  • \(\lambda^B_1 = 0.6750, \quad \lambda^B_2 = -0.6750\)
  • \(\lambda^{AB}_{11} = -0.7611, \quad \lambda^{AB}_{12} = 0.7611, \quad \lambda^{AB}_{21} = 0.7611, \quad \lambda^{AB}_{22} = -0.7611\)

Hitung Odds Ratio dan Interval Kepercayaan

\[ \text{OR} = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} = \frac{15 \cdot 5}{35 \cdot 45} = \frac{75}{1575} = 0.0476 \]

Log odds ratio:

\[ \log(\text{OR}) = \log(0.0476) = -3.048 \]

Standard Error (SE):

\[ SE = \sqrt{ \frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} } \\ = \sqrt{ \frac{1}{15} + \frac{1}{35} + \frac{1}{45} + \frac{1}{5} } = \sqrt{0.3175} = 0.5634 \]

Confidence Interval untuk \(\log(\text{OR})\)

\[ \log(\text{OR}) \pm 1.96 \cdot SE = -3.048 \pm 1.96 \cdot 0.5634 \\ = (-3.048 - 1.104, -3.048 + 1.104) = (-4.152, -1.944) \]

Back transform to get Cl for OR:

\[ Lower=exp(-4.152)=0.0157 \\ Upper=exp(-1.944)=0.1431 \]

Jadi OR=0.0476 (95% Cl: 0.0157-0.1431)

Fitting Model Log-Linear dengan R

tabel <- matrix(c(15, 35, 45, 5), nrow = 2, byrow = TRUE)
colnames(tabel) <- c("Terinfeksi", "Tidak Terinfeksi")
rownames(tabel) <- c("Vaksin A", "Placebo")
tabel
##          Terinfeksi Tidak Terinfeksi
## Vaksin A         15               35
## Placebo          45                5
datam <- as.data.frame(as.table(tabel))
colnames(datam) <- c("Vaksin", "Status", "Freq")
data
##              Preferensi
## Jenis_Kelamin Produk A Produk B
##     Laki-laki       30       20
##     Perempuan       50       40
# Model tanpa interaksi
fit_no_inter <- glm(Freq ~ Vaksin + Status, family = poisson, data = datam)
summary(fit_no_inter)
## 
## Call:
## glm(formula = Freq ~ Vaksin + Status, family = poisson, data = datam)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             3.401e+00  1.633e-01  20.828   <2e-16 ***
## VaksinPlacebo           5.633e-12  2.000e-01   0.000    1.000    
## StatusTidak Terinfeksi -4.055e-01  2.041e-01  -1.986    0.047 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 45.035  on 3  degrees of freedom
## Residual deviance: 41.008  on 1  degrees of freedom
## AIC: 66.091
## 
## Number of Fisher Scoring iterations: 5
# Model dengan interaksi
fit_inter <- glm(Freq ~ Vaksin * Status, family = poisson, data = datam)
summary(fit_inter)
## 
## Call:
## glm(formula = Freq ~ Vaksin * Status, family = poisson, data = datam)
## 
## Coefficients:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                            2.7081     0.2582  10.488  < 2e-16 ***
## VaksinPlacebo                          1.0986     0.2981   3.685 0.000229 ***
## StatusTidak Terinfeksi                 0.8473     0.3086   2.746 0.006041 ** 
## VaksinPlacebo:StatusTidak Terinfeksi  -3.0445     0.5634  -5.403 6.54e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4.5035e+01  on 3  degrees of freedom
## Residual deviance: 7.9936e-15  on 0  degrees of freedom
## AIC: 27.084
## 
## Number of Fisher Scoring iterations: 3

Interpretasi Parameter

Model: glm(Freq ~ Vaksin * Status, family = poisson, data = datam)

Artinya, kita memodelkan frekuensi (Freq) sebagai variabel respons, dengan prediktor: - Vaksin - Status - Interaksi antara Vaksin dan Status

Menggunakan distribusi Poisson, yang cocok untuk data berupa jumlah kejadian.

1. (Intercept) = 2.7081

Ini adalah log-rata-rata kejadian (log(Freq)) untuk kelompok referensi:

  • Vaksin = Kontrol
  • Status = Terinfeksi

Interpretasi dalam skala asli:

  • Frekuensi = exp(2.7081) ≈ 15 orang
  • Artinya, rata-rata jumlah kejadian pada kelompok kontrol yang terinfeksi adalah sekitar 15 orang.

2. VaksinPlacebo = 1.098

Efek pemberian vaksin placebo dibandingkan kontrol, untuk kelompok terinfeksi.

Interpretasi:

  • Log-frekuensi meningkat 1.0986
  • Frekuensi = exp(1.0986) ≈ 3 kali lebih besar
  • Artinya, jumlah infeksi di kelompok placebo 3 kali lebih tinggi dibanding kelompok kontrol (di antara yang terinfeksi).

3. StatusTidak Terinfeksi = 0.8473

Efek dari tidak terinfeksi dibanding terinfeksi, pada kelompok vaksin kontrol.

Interpretasi:

  • Log-frekuensi meningkat 0.8473
  • rekuensi = exp(0.8473) ≈ 2.33 kali lebih besar
  • Artinya, dalam kelompok kontrol, mereka yang tidak terinfeksi memiliki frekuensi 2.33 kali lebih tinggi dibanding yang terinfeksi.

4. Interaksi VaksinPlacebo:StatusTidak Terinfeksi = -3.0445

Efek gabungan dari placebo dan tidak terinfeksi, dibanding ekspektasi jika efeknya bersifat additif.

Interpretasi:

  • Log-frekuensi menurun 3.0445
  • Frekuensi = exp(-3.0445) ≈ 0.0476 (hanya 4.76%)
  • Artinya, kombinasi placebo dan tidak terinfeksi menunjukkan penurunan besar, jauh lebih kecil dari yang diperkirakan tanpa interaksi.
  • Menunjukkan adanya interaksi negatif yang kuat.

Signifikansi Statistik

Semua koefisien memiliki p-value < 0.01, yang berarti:

✅ Efek dari vaksin, status infeksi, dan interaksi keduanya signifikan secara statistik.


Goodness-of-Fit

  • Residual deviance ≈ 0 → Model fit dengan sangat baik
  • AIC = 27.084 → Sangat rendah, menunjukkan model yang efisien
  • Fisher scoring iterations = 3 → Model konvergen dengan cepat

Kesimpulan

  • Vaksin placebo meningkatkan risiko infeksi secara signifikan.
  • Tidak terinfeksi meningkatkan frekuensi hanya dalam kelompok kontrol.
  • Kombinasi placebo dan tidak terinfeksi justru menyebabkan penurunan frekuensi yang sangat besar, menunjukkan interaksi negatif yang kuat.
  • Model sangat cocok dengan data, dan semua efek signifikan.

Analisis Data Tabel Kontingensi 2x3

Suatu survei dilakukan untuk mengetahui hubungan antara Jenis Kelamin (Laki-laki/Perempuan) dan Tingkat Kebiasaan Mengonsumsi Makanan Asin:


Jarang
Kadang-kadang Sering
Laki-laki 12 20 8
Perempuan 18 24 10

Bentuk Model LogLinear Untuk Tabel 2x3

Fitting Model LogLinear di R

# Membuat data frame dari tabel
tabel2x3 <- matrix(c(12, 20, 8, 18, 24, 10), nrow = 2, byrow = TRUE)
colnames(tabel2x3) <- c("Jarang", "Kadang-kadang", "Sering")
rownames(tabel2x3) <- c("Laki-laki", "Perempuan")
tabel2x3
##           Jarang Kadang-kadang Sering
## Laki-laki     12            20      8
## Perempuan     18            24     10
# Ubah menjadi data.frame untuk glm
data2x3 <- as.data.frame(as.table(tabel2x3))
colnames(data2x3) <- c("JenisKelamin", "Tingkat", "Freq")
data2x3
##   JenisKelamin       Tingkat Freq
## 1    Laki-laki        Jarang   12
## 2    Perempuan        Jarang   18
## 3    Laki-laki Kadang-kadang   20
## 4    Perempuan Kadang-kadang   24
## 5    Laki-laki        Sering    8
## 6    Perempuan        Sering   10
# Model log-linear tanpa interaksi (asumsi independen)
fit_no_inter <- glm(Freq ~ JenisKelamin + Tingkat, family = poisson, data = data2x3)
summary(fit_no_inter)
## 
## Call:
## glm(formula = Freq ~ JenisKelamin + Tingkat, family = poisson, 
##     data = data2x3)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             2.5683     0.2179  11.789   <2e-16 ***
## JenisKelaminPerempuan   0.2624     0.2103   1.248   0.2122    
## TingkatKadang-kadang    0.3830     0.2368   1.618   0.1058    
## TingkatSering          -0.5108     0.2981  -1.713   0.0866 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 13.06443  on 5  degrees of freedom
## Residual deviance:  0.22527  on 2  degrees of freedom
## AIC: 35.26
## 
## Number of Fisher Scoring iterations: 3
# Model log-linear dengan interaksi (untuk cek asosiasi)
fit_inter <- glm(Freq ~ JenisKelamin * Tingkat, family = poisson, data = data2x3)
summary(fit_inter)
## 
## Call:
## glm(formula = Freq ~ JenisKelamin * Tingkat, family = poisson, 
##     data = data2x3)
## 
## Coefficients:
##                                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                  2.4849     0.2887   8.608   <2e-16
## JenisKelaminPerempuan                        0.4055     0.3727   1.088    0.277
## TingkatKadang-kadang                         0.5108     0.3651   1.399    0.162
## TingkatSering                               -0.4055     0.4564  -0.888    0.374
## JenisKelaminPerempuan:TingkatKadang-kadang  -0.2231     0.4802  -0.465    0.642
## JenisKelaminPerempuan:TingkatSering         -0.1823     0.6032  -0.302    0.762
##                                               
## (Intercept)                                ***
## JenisKelaminPerempuan                         
## TingkatKadang-kadang                          
## TingkatSering                                 
## JenisKelaminPerempuan:TingkatKadang-kadang    
## JenisKelaminPerempuan:TingkatSering           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  1.3064e+01  on 5  degrees of freedom
## Residual deviance: -9.0719e-30  on 0  degrees of freedom
## AIC: 39.034
## 
## Number of Fisher Scoring iterations: 3

Interpretasi

Interpretasi Model Regresi Poisson

Model: glm(Freq ~ JenisKelamin * Tingkat, family = poisson, data = data2x3)

Model ini memodelkan jumlah kejadian (Freq) berdasarkan: - JenisKelamin (Laki-laki, Perempuan) - Tingkat (frekuensi perilaku: Jarang, Kadang-kadang, Sering) - Dan interaksi antara keduanya

Distribusi: Poisson (cocok untuk data hitungan).

1. (Intercept) = 2.4849

Ini adalah log-rata-rata frekuensi kejadian untuk kelompok referensi: - JenisKelamin = Laki-laki - Tingkat = Jarang

Interpretasi:

  • Frekuensi = exp(2.4849) ≈ 11.99
  • Artinya, laki-laki dengan tingkat perilaku jarang rata-rata memiliki sekitar 12 kejadian.

2. JenisKelaminPerempuan = 0.4055

Efek menjadi perempuan (dibanding laki-laki) dalam kelompok Tingkat = Jarang.

Interpretasi:

  • exp(0.4055) ≈ 1.50
  • Artinya, perempuan dengan tingkat “jarang” memiliki 50% lebih banyak kejadian daripada laki-laki, namun tidak signifikan secara statistik (p = 0.277).

3. TingkatKadang-kadang = 0.5108

Efek dari frekuensi “kadang-kadang” dibanding “jarang”, dalam kelompok laki-laki.

Interpretasi:

  • exp(0.5108) ≈ 1.67
  • Artinya, laki-laki dengan tingkat “kadang-kadang” memiliki jumlah kejadian 1.67 kali lebih banyak dibanding “jarang”. Namun tidak signifikan (p = 0.162).

4. TingkatSering = -0.4055

Efek dari frekuensi “sering” dibanding “jarang”, dalam kelompok laki-laki.

Interpretasi:

  • exp(-0.4055) ≈ 0.67
  • Artinya, laki-laki dengan tingkat “sering” justru memiliki jumlah kejadian lebih rendah, yaitu sekitar 67% dari yang “jarang”. Tidak signifikan (p = 0.374).

5. Interaksi JenisKelaminPerempuan:TingkatKadang-kadang = -0.2231

Efek tambahan dari kombinasi perempuan dan kadang-kadang dibanding ekspektasi jika efeknya additif.

Interpretasi:

  • exp(-0.2231) ≈ 0.80
  • Artinya, ada penurunan sebesar 20%, namun efek ini tidak signifikan (p = 0.642).

6. Interaksi JenisKelaminPerempuan:TingkatSering = -0.1823

Efek tambahan dari kombinasi perempuan dan sering.

Interpretasi:

  • exp(-0.1823) ≈ 0.83
  • Artinya, penurunan sekitar 17%, namun tidak signifikan (p = 0.762).

Signifikansi Statistik

  • Hanya intercept yang signifikan (p < 0.001)
  • Semua prediktor dan interaksi memiliki p-value > 0.1
  • Artinya, tidak ada bukti kuat bahwa JenisKelamin, Tingkat, atau interaksinya mempengaruhi f

Model LogLinear Tuga 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 LogLinear Untuk Tabel Tiga Arah

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

  1. Model Saturated

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

  2. Model Homogen

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

  3. Model Consitional

    • Conditional Pada X

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

    • Conditional Pada Y

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

    • Conditional Pada Z

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

  4. Model Joint Independence

    • Independensi antara X dan Y

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

    • Independensi antara X dan Z

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

    • Independensi antara Y dan Z

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

  5. Model Tanpa Interaksi

    \[ \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 homogenous.
  2. Pengujian Interaksi Dua Arah (XY, XZ, YZ):

    • Bandingkan model homogenous dengan model conditional.

    • Bandingkan model conditional dengan model joint independence.

    • Bandingkan model joint independence dengan model tanpa interaksi.

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

library("epitools")
library("DescTools")
## Warning: package 'DescTools' was built under R version 4.3.3
library("lawstat")
## Warning: package 'lawstat' was built under R version 4.3.3
# Input data sesuai tabel praktikum
z.fund <- factor(rep(c("1fund", "2mod", "3lib"), each = 4))
x.sex  <- factor(rep(c("1M", "2F"), each = 2, times = 3))
y.fav  <- factor(rep(c("1fav", "2opp"), times = 6))
counts <- c(128, 32, 123, 73, 182, 56, 168, 105, 119, 49, 111, 70)

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

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

##=============================##
# Penentuan kategori reference
##=============================##
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:

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

Interpretasi output model saturated

Ringkasan model

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

Hasil estimasi koefisien

Parameter Estimate Std. Error z value Pr(> z
(Intercept) 4.25 0.12 35.55 <2e-16*** 70.00
x.sex1M -0.36 0.19 -1.92 0.055 . 0.70
y.fav1fav 0.46 0.15 3.02 0.0025 ** 1.59
z.fund1fund 0.04 0.17 0.25 0.80 1.04
z.fund2mod 0.41 0.15 2.63 0.0086 ** 1.50
x.sex1M:y.fav1fav 0.43 0.23 1.87 0.062 . 1.53
x.sex1M:z.fund1fund -0.47 0.28 -1.66 0.097 . 0.63
x.sex1M:z.fund2mod -0.27 0.25 -1.09 0.28 0.76
y.fav1fav:z.fund1fund 0.06 0.21 0.29 0.78 1.06
y.fav1fav:z.fund2mod 0.01 0.20 0.05 0.96 1.01
x.sex1M:y.fav1fav:z.fund1fund 0.44 0.34 1.30 0.19 1.55
x.sex1M:y.fav1fav:z.fund2mod 0.28 0.30 0.94 0.35 1.33

Goodness-of-fit

Model Homogenous

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

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

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

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

Pengujian ini menggunakan residual deviance dari kedua model (saturated dan homogenous).

Langkah-langkah pangujian

Hipotesis

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

Hitung selisih deviance

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

Hitung derajat bebas

# Derajat bebas = db model homogenous - db model saturated
derajat.bebas <- (model_homogenous$df.residual - model_saturated$df.residual)
derajat.bebas
## [1] 2

Chi-square tabel (alpha=0.05)

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"

Interpretasi Pada taraf nyata 5%, belum cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi tiga arah antara jenis kelamin, fundamentalisme, dan pendapat mengenai hukuman mati.

Catatan:

  • Model pengurang adalah model yang lebih lengkap (lebih banyak parameter, df lebih kecil), yaitu model saturated.
  • Derajat bebas dihitung dari selisih derajat bebas model homogenous dan saturated.
  • Keputusan berdasarkan perbandingan deviance model dengan chi-square tabel.

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 hipotesis

# Deviance of Model
Deviance.model <- model_conditional_X$deviance - model_homogenous$deviance   # model_conditional_X: conditional on X, model_homogenous: homogenous
Deviance.model
## [1] 2.132302
# Selisih deviance antar model
Deviance.model <- model_conditional_X$deviance - model_homogenous$deviance
Deviance.model
## [1] 2.132302
# Chi Square tabel dengan alpha = 0.05
derajat.bebas <- (4 - 2)
derajat.bebas
## [1] 2
chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 5.991465
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan
## [1] "Terima"

Dengan taraf nyata 5%, belum cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi antara pendapat tentang hukuman mati dan fundamentalisme. 

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
# Deviance of Model
Deviance.model <- model_conditional_Y$deviance - model_homogenous$deviance   # model_conditional_Y: conditional on Y, model_homogenous: homogenous
Deviance.model
## [1] 1.122315
derajat.bebas <- (4 - 2)
derajat.bebas
## [1] 2
chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 5.991465
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan
## [1] "Terima"

Dengan taraf nyata 5%, belum cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi antara jenis kelamin dan fundamentalisme. 

4.5 UJI MODEL INTERAKSI DUA 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
# Deviance of Model
Deviance.model <- model_conditional_Z$deviance - model_homogenous$deviance   # model_conditional_Z: conditional on Z, model_homogenous: homogenous
Deviance.model
## [1] 27.93095
derajat.bebas <- (3 - 2)
derajat.bebas
## [1] 1
chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 3.841459
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan
## [1] "Tolak"

Dengan taraf nyata 5%, ada interaksi antara jenis kelamin dan pendapat tentang hukuman mati.

Pemilihan Model Terbaik

Ringkasan Model LogLinear

Model Parameter Deviance Jumlah Parameter df AIC
Saturated \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk} \] 0.0000 12 0 100.14
Homogenous \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \] 1.798 10 2 97.934
Conditional on X \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} \] 3.9303 8 4 96.067
Conditional on Y \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk} \] 2.9203 8 4 95.057
Conditional on Z \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \] 29.729 9 3 123.87

Ringkasan Pengujian Interaksi 3 Arah dan 2 Arah

Interaksi Pengujian Δ deviance Δ df Chi-square Tabel Keputusan Keterangan
XYZ Saturated vs Homogenous 1.798 2 5.991 Tidak Tolak H0 tidak ada interaksi
YZ Conditional on X vs Homogenous 2.1323 2 5.991 Tidak Tolak H0 tidak ada interaksi
XZ Conditional on Y vs Homogenous 1.1223 2 5.991 Tidak Tolak H0 tidak ada interaksi
XY Conditional on Z vs Homogenous 27.931 1 3.841 Tolak H0 ada interaksi

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

Dari summary model diatas terlihat bahwa best model memiliki AIC yang lebih rendah dibandingkan saturated, homogeneous, dan conditional model.

Interpretasi Koefisien Model Terbaik

# Interpretasi koefisien model terbaik
data.frame(
  koef = bestmodel$coefficients,
  exp_koef = exp(bestmodel$coefficients)
)
##                          koef   exp_koef
## (Intercept)        4.26517861 71.1776316
## x.sex1M           -0.59344782  0.5524194
## y.fav1fav          0.48302334  1.6209677
## z.fund1fund        0.01985881  1.0200573
## z.fund2mod         0.38129767  1.4641834
## x.sex1M:y.fav1fav  0.65845265  1.9318008

Interpretasi Koefisien Model Terbaik

  • \(\exp(\lambda^X_{1M}) = \exp(-0.593) = 0.552\)
    Tanpa memperhatikan fundamentalisme dan pendapat mengenai hukuman mati, peluang seseorang berjenis kelamin laki-laki adalah 0,552 kali dibandingkan perempuan.
    Atau, peluang seseorang berjenis kelamin perempuan adalah \(\frac{1}{0.552} = 1.81\) kali dibandingkan laki-laki.
  • \(\exp(\lambda^Y_{1,\text{fav}}) = \exp(0.483) = 1.621\)
    Tanpa memperhatikan jenis kelamin dan fundamentalisme, peluang seseorang mendukung hukuman mati adalah 1,621 kali dibandingkan yang menolak.
  • \(\exp(\lambda^Z_{1,\text{fund}}) = \exp(0.01986) = 1.02\)
    Tanpa memperhatikan jenis kelamin dan pendapat mengenai hukuman mati, peluang seseorang fundamentalist adalah 1,02 kali dibandingkan liberal.
  • \(\exp(\lambda^Z_{2,\text{mod}}) = \exp(0.381) = 1.464\)
    Tanpa memperhatikan jenis kelamin dan pendapat mengenai hukuman mati, peluang seseorang moderate adalah 1,464 kali dibandingkan liberal.
  • \(\exp(\lambda^{XY}_{1M,1\text{fav}}) = \exp(0.658) = 1.932\)
    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
)
##     Fund sex favor counts    fitted
## 1  1fund  1M  1fav    128 125.59539
## 2  1fund  1M  2opp     32  40.10855
## 3  1fund  2F  1fav    123 117.69079
## 4  1fund  2F  2opp     73  72.60526
## 5   2mod  1M  1fav    182 180.27878
## 6   2mod  1M  2opp     56  57.57155
## 7   2mod  2F  1fav    168 168.93257
## 8   2mod  2F  2opp    105 104.21711
## 9   3lib  1M  1fav    119 123.12582
## 10  3lib  1M  2opp     49  39.31990
## 11  3lib  2F  1fav    111 115.37664
## 12  3lib  2F  2opp     70  71.17763

Manual

Estimasi Nilai Harapan \(\hat{\mu}_{ijk}\)

\[ \hat{\mu}_{111} = \exp(\lambda + \lambda^X_{1M} + \lambda^Y_{1fav} + \lambda^Z_{1fund} + \lambda^{XY}_{1M,1fav}) \\ = \exp(4.265 - 0.593 + 0.483 + 0.01986 + 0.658) \\ = \exp(4.833) = 125.595 \]

\[ \hat{\mu}_{112} = \exp(\lambda + \lambda^X_{1M} + \lambda^Y_{1fav} + \lambda^Z_{2mod} + \lambda^{XY}_{1M,1fav}) \\ = \exp(4.265 - 0.593 + 0.483 + 0.381 + 0.658) \\ = \exp(5.195) = 180.279 \]

\[ \hat{\mu}_{113} = \exp(\lambda + \lambda^X_{1M} + \lambda^Y_{1fav} + \lambda^Z_{3lib} + \lambda^{XY}_{1M,1fav}) \\ = \exp(4.265 - 0.593 + 0.483 + 0 + 0.658) \\ = \exp(4.813) = 123.126 \]

\[ \hat{\mu}_{121} = \exp(\lambda + \lambda^X_{1M} + \lambda^Y_{2opp} + \lambda^Z_{1fund} + \lambda^{XY}_{1M,2opp}) \\ = \exp(4.265 - 0.593 + 0 + 0.01986 + 0) \\ = \exp(3.692) = 40.109 \]

\[ \hat{\mu}_{122} = \exp(\lambda + \lambda^X_{1M} + \lambda^Y_{2opp} + \lambda^Z_{2mod} + \lambda^{XY}_{1M,2opp}) \\ = \exp(4.265 - 0.593 + 0 + 0.381 + 0) \\ = \exp(4.053) = 57.572 \]

\[ \hat{\mu}_{123} = \exp(\lambda + \lambda^X_{1M} + \lambda^Y_{2opp} + \lambda^Z_{3lib} + \lambda^{XY}_{1M,2opp}) \\ = \exp(4.265 - 0.593 + 0 + 0 + 0) \\ = \exp(3.672) = 39.320 \]

\[ \hat{\mu}_{211} = \exp(\lambda + \lambda^X_{2F} + \lambda^Y_{1fav} + \lambda^Z_{1fund} + \lambda^{XY}_{2F,1fav}) \\ = \exp(4.265 + 0 + 0.483 + 0.01986 + 0) \\ = \exp(4.768) = 117.691 \]

\[ \hat{\mu}_{212} = \exp(\lambda + \lambda^X_{2F} + \lambda^Y_{1fav} + \lambda^Z_{2mod} + \lambda^{XY}_{2F,1fav}) \\ = \exp(4.265 + 0 + 0.483 + 0.381 + 0) \\ = \exp(5.129) = 168.933 \]

\[ \hat{\mu}_{213} = \exp(\lambda + \lambda^X_{2F} + \lambda^Y_{1fav} + \lambda^Z_{3lib} + \lambda^{XY}_{2F,1fav}) \\ = \exp(4.265 + 0 + 0.483 + 0 + 0) \\ = \exp(4.748) = 115.377 \]

\[ \hat{\mu}_{221} = \exp(\lambda + \lambda^X_{2F} + \lambda^Y_{2opp} + \lambda^Z_{1fund} + \lambda^{XY}_{2F,2opp}) \\ = \exp(4.265 + 0 + 0 + 0.01986 + 0) \\ = \exp(4.285) = 72.605 \]

\[ \hat{\mu}_{222} = \exp(\lambda + \lambda^X_{2F} + \lambda^Y_{2opp} + \lambda^Z_{2mod} + \lambda^{XY}_{2F,2opp}) \\ = \exp(4.265 + 0 + 0 + 0.381 + 0) \\ = \exp(4.646) = 104.236 \]

\[ \hat{\mu}_{223} = \exp(\lambda + \lambda^X_{2F} + \lambda^Y_{2opp} + \lambda^Z_{3lib} + \lambda^{XY}_{2F,2opp}) \\ = \exp(4.265 + 0 + 0 + 0 + 0) \\ = \exp(4.265) = 71.178 \]

Referensi