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.
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:
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.
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.
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.
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.
Data kategori adalah data yang mengelompokkan objek ke dalam kategori yang berbeda, yang tidak diukur secara numerik tetapi diklasifikasikan.
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.
Analisis data kategori memiliki manfaat yang luas dalam berbagai sektor:
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.
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.
Pemasaran dan Bisnis
Segmentasi pasar berdasarkan kategori seperti jenis kelamin, wilayah,
atau preferensi produk memungkinkan perusahaan mengembangkan strategi
pemasaran yang lebih tepat sasaran.
Pendidikan
Kategori seperti tingkat pendidikan, jenis sekolah, atau status siswa
digunakan untuk mengevaluasi kebijakan pendidikan dan hasil
belajar.
Kebijakan Publik dan Pemerintahan
Pemerintah menggunakan data kategori untuk perencanaan dan pengambilan
kebijakan, seperti klasifikasi wilayah menurut tingkat kemiskinan atau
pekerjaan penduduk.
Keamanan dan Kriminalitas
Data kategori digunakan untuk mengklasifikasikan jenis kejahatan, profil
pelaku, dan lokasi kejadian. Analisis ini mendukung kebijakan preventif
dan tindakan hukum.
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 (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 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 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 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.
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 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 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 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:
# 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 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:
# 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 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:
# 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 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:
# 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 adalah aspek penting dalam analisis data kategori karena menentukan bagaimana data dikumpulkan sehingga hasil analisis valid dan dapat digeneralisasi.
Desain sampling prospektif mengamati subjek dari titik awal dan mengikuti ke depan dalam waktu untuk melihat hasil yang terjadi.
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.
Studi kohort mengikuti sekelompok orang yang memiliki karakteristik tertentu selama periode waktu untuk melihat apakah mereka mengalami kejadian tertentu (misal: penyakit).
Desain sampling retrospektif melihat ke belakang pada data atau kejadian yang sudah terjadi.
Desain yang membandingkan kelompok dengan kondisi (kasus) dan tanpa kondisi (kontrol) untuk menemukan faktor risiko.
Mirip studi kohort, tetapi data dikumpulkan dari catatan masa lalu.
| 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 adalah alat dasar dalam analisis data kategori untuk melihat hubungan antara dua variabel kategori yang masing-masing memiliki dua level.
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 mencakup peluang bersama, marginal, dan bersyarat dari kejadian pada tabel.
Peluang bersama adalah probabilitas bahwa dua kejadian terjadi bersama-sama.
Contoh: \[ P(\text{Diare} \cap \text{Tidak Memenuhi}) = \frac{36}{72} = 0.5 \]
\[ 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} \]
\[ P(\text{Diare} | \text{Tidak Memenuhi}) = \frac{36}{43}, \quad P(\text{Diare} | \text{Memenuhi}) = \frac{14}{29} \]
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\) |
\[ RD = \frac{a}{a+b} - \frac{c}{c+d} = \frac{36}{43} - \frac{14}{29} \approx 0.8372 - 0.4828 = 0.3544 \]
\[ RR = \frac{\frac{36}{43}}{\frac{14}{29}} = \frac{0.8372}{0.4828} \approx 1.7345 \]
\[ OR = \frac{ad}{bc} = \frac{36 \times 15}{7 \times 14} = \frac{540}{98} \approx 5.51 \]
# 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
| 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:
Estimasi titik digunakan untuk mendapatkan nilai terbaik dari parameter populasi berdasarkan data sampel.
Estimasi interval memberikan rentang nilai yang kemungkinan besar mencakup parameter populasi.
Uji proporsi dilakukan untuk membandingkan proporsi antar kelompok dalam tabel kontingensi 2x2.
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 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 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)} \]
Misalkan data dari studi kasus air bersih dan diare:
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
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 Chi-Square digunakan untuk menguji apakah ada hubungan antara dua variabel kategorikal.
\[ \chi^2 = \sum \frac{(O - E)^2}{E} \]
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.
Digunakan untuk memeriksa kontribusi tiap sel terhadap nilai chi-square total.
\[ G^2 = 2 \sum O \log\left(\frac{O}{E}\right) \]
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.
Pearson Residual \[ R_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}} \]
Standardized Residual Digunakan untuk membandingkan residual antar sel.
Sel dengan residual besar (misalnya > 2 atau < -2) dianggap sebagai outlier.
chisq <- chisq.test(observed)
chisq$residuals
## [,1] [,2]
## [1,] 1.123406 -1.693598
## [2,] -1.367956 2.062271
Interpretasi
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.
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.
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:
Studi epidemologi: Hubungan merokok (\(X\)) dan kanker paru-paru (\(Y\)) dengan kontrol usia (\(Z\)).
Penelitian sosial: Hubungan antara ras tersangka (\(X\)) dan keputusan hukuman mati (\(Y\)), dengan mempertimbangkan ras korban (\(Z\)).
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 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\).
| 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 |
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
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.
\[ 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 31.0% berpendidikan SMP
Sekitar 35.5% berpendidikan SMA
\[ 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
Tabel ini menunjukkan probabilitas kondisi dua variabel diberikan nilai tertentu dari variabel ketiga.
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:
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"
Dua variabel dikatakan conditionally independent jika tidak ada hubungan antara keduanya setelah mengendalikan variabel ketiga.
Terkadang hubungan marginal antara \(X\) dan \(Y\) bisa berbeda dari hubungan pada setiap strata \(Z\), fenomena ini disebut Simpson’s paradox.
Untuk menguji apakah \(X\) dan \(Y\) independen secara bersyarat terhadap \(Z\).
Uji Chi-Square dapat dilakukan pada setiap strata \(Z\) dan dibandingkan.
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:
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
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()
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) 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
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 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
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 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
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).
curve(log(x / (1 - x)), from = 0.01, to = 0.99, col = "red", lwd = 2, ylab = "Nilai Link", main = "Fungsi Link")
curve(log(x), add = TRUE, col = "blue", lwd = 2)
legend("topleft", legend = c("Logit", "Log"), col = c("red", "blue"), lty = 1, lwd = 2)
Interpretasi: - Garis merah (logit) menunjukkan fungsi link dari regresi logistik: \(\log(\pi / (1 - \pi))\). Kurva ini memiliki bentuk cembung yang lebih curam di sekitar nilai probabilitas menengah (0.5). - Garis biru (log) merupakan fungsi link dari regresi Poisson: \(\log(\mu)\). Kurvanya lebih datar, dan hanya terdefinisi untuk \(x > 0\). - Perbandingan ini menunjukkan bahwa fungsi logit memperbesar perubahan di sekitar nilai tengah, sedangkan fungsi log menangkap pertumbuhan eksponensial untuk data hitung.
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
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.
Dalam kerangka GLM, ekspektasi dan variansi dari variabel respon \(Y\) tergantung pada bentuk distribusi dari keluarga eksponensial dan fungsi link yang digunakan.
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\)
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).
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
Setelah membangun model, perlu dilakukan pemeriksaan diagnostik untuk mengevaluasi: - Kesesuaian model terhadap data - Deteksi outlier atau leverage points - Validitas asumsi model
plot(model_logit) # Menampilkan 4 grafik diagnostik penting
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) \]
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
Regresi Poisson digunakan untuk data hitungan. Asumsi utama adalah: - Variansi sama dengan rata-rata (\(Var(Y) = \mu\))
summary(model_pois)$deviance
## [1] 98.65831
AIC(model_pois)
## [1] 373.3214
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")
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.
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
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
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:
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.
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.
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.
konsumsi_garam.L: 1.11791 (p < 0.001)
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).
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).
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.
Langkah seleksi stepwise membantu menentukan kombinasi prediktor terbaik untuk model regresi logistik. Kriteria seperti AIC digunakan untuk menilai kualitas model.
# 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
ROC dan AUC digunakan untuk mengevaluasi kemampuan model dalam mengklasifikasikan data biner dengan berbagai ambang batas prediksi.
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
Karena regresi logistik tidak memiliki R² seperti regresi linear, maka digunakan alternatif pseudo R-squared, seperti:
Nilainya berkisar antara 0 dan 1; semakin tinggi, semakin baik kesesuaian model.
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:
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.
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.
Model yang paling sederhana dengan akurasi terbaik lebih disukai. Jangan menambahkan variabel yang tidak signifikan atau tidak berguna secara praktis.
Akurasi dihitung dari proporsi klasifikasi benar terhadap total: \[ \text{Akurasi} = \frac{TP + TN}{TP + TN + FP + FN} \]
Digunakan untuk mengukur goodness-of-fit secara global:
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.
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.
Dalam bidang kesehatan, kita ingin memodelkan kategori tekanan darah seseorang (Normal, Prehipertensi, atau Hipertensi) berdasarkan:
Tujuannya adalah mengetahui apakah usia dan konsumsi garam berpengaruh terhadap kemungkinan seseorang termasuk dalam kategori tekanan darah tertentu.
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 dilakukan dengan metode Maximum Likelihood Estimation (MLE). Proses ini mengoptimalkan fungsi likelihood untuk mendapatkan parameter terbaik yang menjelaskan 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)
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
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:
Contoh:
usia untuk kategori Hipertensi
signifikan dan positif, maka semakin tua usia, semakin besar kemungkinan
masuk ke kategori Hipertensi dibanding Normal.konsumsi_garamTinggi signifikan dan positif pada
kategori Prehipertensi, artinya orang yang konsumsi garam tinggi lebih
berisiko mengalami Prehipertensi dibanding Normal.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.
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:
Evaluasi seperti akurasi dan signifikansi parameter memberikan gambaran apakah model bisa digunakan untuk prediksi maupun pengambilan keputusan kebijakan kesehatan.
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.
Koefisien \(\beta\) menunjukkan pengaruh variabel terhadap kemungkinan berada pada kategori yang lebih tinggi.
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)
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
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:
predict(model_ord, newdata = data.frame(usia = 60, konsumsi_garam = "Tinggi"), type = "probs")
## Normal Prehipertensi Hipertensi
## 0.97029652 0.02616978 0.00353370
Interpretasi:
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 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.
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 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 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 dalam model log-linear dapat digunakan untuk mengukur kekuatan asosiasi antara dua variabel. Nilai odds ratio = 1 menunjukkan tidak ada asosiasi.
Interpretasi:
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 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
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.
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.
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:
Tabel Kontingensi: penyajian frekuensi gabungan dari dua atau lebih variabel kategorik.
Model Loglinear: digunakan untuk memodelkan struktur asosiasi di dalam tabel kontingensi tanpa menganggap ada variabel dependen.
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 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 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 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 untuk tabel 2x2:
\\\Rumus///
Interpretasi nilai OR:
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
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
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
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///
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).
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} \]
Diberikan data:
| Vaksin | Terinfeksi | Tidak Terinfeksi |
|---|---|---|
| Vaksin A | 15 | 35 |
| Placebo | 45 | 5 |
Model LogLinear pada tabel 2x2:
dengan constraint sum-to-zero:
Misalkan:
Observasi:
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 \]
\[ \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 \]
\[ \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 \]
\[ \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 \]
\[ \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 \]
\[ \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)
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
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:
Interpretasi dalam skala asli:
2. VaksinPlacebo = 1.098
Efek pemberian vaksin placebo dibandingkan kontrol, untuk kelompok terinfeksi.
Interpretasi:
3. StatusTidak Terinfeksi = 0.8473
Efek dari tidak terinfeksi dibanding terinfeksi, pada kelompok vaksin kontrol.
Interpretasi:
4. Interaksi VaksinPlacebo:StatusTidak Terinfeksi = -3.0445
Efek gabungan dari placebo dan tidak terinfeksi, dibanding ekspektasi jika efeknya bersifat additif.
Interpretasi:
Semua koefisien memiliki p-value < 0.01, yang berarti:
✅ Efek dari vaksin, status infeksi, dan interaksi keduanya signifikan secara statistik.
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 |
# 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
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:
2. JenisKelaminPerempuan = 0.4055
Efek menjadi perempuan (dibanding laki-laki) dalam kelompok Tingkat = Jarang.
Interpretasi:
3. TingkatKadang-kadang = 0.5108
Efek dari frekuensi “kadang-kadang” dibanding “jarang”, dalam kelompok laki-laki.
Interpretasi:
4. TingkatSering = -0.4055
Efek dari frekuensi “sering” dibanding “jarang”, dalam kelompok laki-laki.
Interpretasi:
5. Interaksi JenisKelaminPerempuan:TingkatKadang-kadang = -0.2231
Efek tambahan dari kombinasi perempuan dan kadang-kadang dibanding ekspektasi jika efeknya additif.
Interpretasi:
6. Interaksi JenisKelaminPerempuan:TingkatSering = -0.1823
Efek tambahan dari kombinasi perempuan dan sering.
Interpretasi:
Signifikansi Statistik
Pada pembahasan sebelumnya, kita telah memahami bahwa salah satu tujuan utama dari penyusunan model log linear adalah untuk mengestimasi parameter-parameter yang menjelaskan hubungan di antara variabel-variabel kategorik.
Pada materi kali ini, kita akan membahas model log linear yang lebih kompleks, yaitu model log linear untuk tabel kontingensi tiga arah. Model ini melibatkan tiga variabel kategorik, sehingga kemungkinan interaksi yang dapat terjadi di dalam model pun menjadi lebih banyak. Dalam konteks ini, interaksi paling tinggi yang dapat dimodelkan adalah interaksi tiga arah, yaitu interaksi yang melibatkan ketiga variabel secara bersamaan.
Model log-linear yang melibatkan tiga variabel kategorik (misal: X, Y, dan Z) dapat dibangun dalam berbagai bentuk model, tergantung pada tingkat interaksi yang ingin dimasukkan. Berikut adalah beberapa alternatif model log-linear yang umum digunakan:
Model Saturated
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk} \]
Model Homogen
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \]
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} \]
Model Joint Independence
Independensi antara X dan Y
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} \]
Independensi antara X dan Z
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} \]
Independensi antara Y dan Z
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{YZ}_{jk} \]
Model Tanpa Interaksi
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k \]
Dalam analisis model log-linear tiga arah, pengujian interaksi dilakukan untuk mengetahui ada atau tidaknya interaksi antar variabel. Pengujian ini dilakukan secara bertahap, dimulai dari tingkat interaksi tertinggi ke yang lebih rendah. Untuk model log-linear dengan tiga peubah (X, Y, dan Z), tahapan pengujian meliputi:
Pengujian Interaksi Tiga Arah (XYZ):
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
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 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
Pengujian ini menggunakan residual deviance dari kedua model (saturated dan homogenous).
Hipotesis
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 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.
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.
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.
| 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 |
| 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 |
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 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
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
# 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
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 \]
MASS dan dokumentasi fungsi
loglm()