EBOOK MENGGUNAKAN R MARKDOWN
Mona Yola Lumban Raja (140610230081)
2025-04-14
Bab 1 Pendahuluan
Data kategori merupakan jenis data yang diklasifikasikan ke dalam kelompok atau kategori tanpa nilai numerik kontinu. Contohnya seperti jenis kelamin, tingkat pendidikan, atau status pekerjaan. Analisis terhadap data jenis ini sangat penting untuk memahami pola dan hubungan antar kelompok.
1.1 Tujuan Analisis Data Kategori
- Menemukan Pola dan Tren: Untuk mengamati kecenderungan tertentu berdasarkan kategori.
- Menganalisis Relasi antar Variabel: Misalnya hubungan antara kebiasaan merokok dan kejadian penyakit.
- Mendukung Pengambilan Keputusan: Seperti penyusunan kebijakan berbasis hasil survei.
- Membangun Model Prediksi: Termasuk penggunaan regresi logistik untuk memprediksi kejadian biner.
1.2 Ruang Lingkup
- Data Nominal: Tidak berurutan (contoh: jenis kelamin).
- Data Ordinal: Memiliki tingkatan (contoh: tingkat kepuasan).
- Data Biner vs Multikategori: Dua kategori (ya/tidak) atau lebih dari dua kategori (misal: warna favorit).
1.3 Manfaat di Berbagai Bidang
- Kesehatan: Menilai risiko penyakit berdasarkan gaya hidup.
- Pendidikan: Evaluasi metode pembelajaran.
- Sosial: Survei opini publik.
- Bisnis: Segmentasi pasar.
- Kriminalitas: Pola jenis kejahatan.
Bab 2 Metode dalam Analisis Data Kategori
Berikut beberapa metode utama yang umum digunakan:
2.1 Tabel Kontingensi dan Uji Chi-Square
Digunakan untuk melihat apakah terdapat hubungan antara dua variabel kategori.
2.2 Regresi Logistik
Digunakan untuk memprediksi probabilitas suatu kejadian (misalnya: beli produk atau tidak) berdasarkan faktor-faktor kategori.
2.3 Correspondence Analysis (CA)
Berguna untuk mengeksplorasi hubungan antar kategori secara grafis.
2.4 Decision Tree dan Random Forest
Metode machine learning berbasis kategori, sangat berguna untuk klasifikasi dan prediksi.
Analisis data kategori merupakan salah satu pilar utama dalam statistik terapan, khususnya ketika data yang dianalisis berbentuk nominal atau ordinal. Metode ini sangat relevan dalam berbagai konteks kehidupan nyata, mulai dari riset kesehatan masyarakat, survei sosial, hingga pengambilan keputusan di bidang pemasaran, pendidikan, dan pemerintahan.
Melalui teknik seperti tabel kontingensi, uji chi-square, regresi logistik, dan pendekatan berbasis machine learning, analisis data kategori memberikan landasan kuat untuk mengeksplorasi dan menguji hubungan antar variabel kategori. Teknik-teknik ini membantu peneliti untuk tidak hanya menggambarkan distribusi data, tetapi juga mengidentifikasi pola keterkaitan dan menguji hipotesis secara sistematis.
Seiring dengan perkembangan teknologi dan ketersediaan data yang semakin melimpah, kemampuan untuk menganalisis data kategori secara efisien menjadi semakin penting. Penggunaan software statistik seperti R, serta integrasi dengan metode modern seperti GLM, model log-linear, atau pendekatan Bayesian, memperluas cakupan dan ketajaman analisis.
Dengan penguasaan metode analisis data kategori, seorang analis tidak hanya mampu memahami fenomena sosial, ekonomi, dan kesehatan secara mendalam, tetapi juga mampu memberikan rekomendasi berbasis data yang lebih akurat dan strategis. Oleh karena itu, pemahaman yang baik terhadap analisis data kategori menjadi aset penting bagi para peneliti, pengambil kebijakan, maupun profesional lintas bidang di era data saat ini.
3 Distribusi Probabilitas dalam Data Kategori
Distribusi probabilitas dalam konteks data kategori merupakan fondasi penting dalam statistik, terutama ketika variabel yang diamati bersifat diskrit dan terbagi dalam beberapa kategori. Distribusi ini memungkinkan kita memahami kemungkinan terjadinya suatu kejadian tertentu dalam kelompok kategori yang ada. Penerapannya meliputi berbagai bidang, mulai dari ilmu kesehatan, psikologi, hingga ekonomi dan ilmu sosial. Variabel acak kategori memiliki distribusi peluang diskrit. Berikut beberapa distribusi penting:
3.1 Distribusi Bernoulli
Distribusi untuk percobaan biner (0 atau 1), dengan rumus:
\[ P(X = x) = p^x (1 - p)^{1 - x},\quad x \in \{0,1\} \] Keterangan: - \(p\): probabilitas sukses (misalnya mendapatkan kepala) - \(1 - p\): probabilitas gagal (misalnya mendapatkan ekor)
# Library
library(knitr)## Warning: package 'knitr' was built under R version 4.3.3library(kableExtra)## Warning: package 'kableExtra' was built under R version 4.3.3library(dplyr)## Warning: package 'dplyr' was built under R version 4.3.3##
## Attaching package: 'dplyr'## The following object is masked from 'package:kableExtra':
##
## group_rows## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# Set seed
set.seed(2025)
# ===== 3.1 Distribusi Bernoulli =====
# Simulasi 10 data Bernoulli (0/1)
bernoulli_sample <- rbinom(n = 10, size = 1, prob = 0.7)
bernoulli_sample## [1] 0 1 1 1 0 1 0 1 1 13.2 Distribusi Binomial
Distribusi Binomial merupakan generalisasi dari distribusi Bernoulli untuk \(n\) percobaan yang independen dan identik. Distribusi ini menggambarkan jumlah keberhasilan dalam sejumlah percobaan tetap dengan peluang yang konstan. Misalnya, menghitung berapa banyak pasien yang sembuh dari sepuluh orang setelah diberikan suatu obat tertentu. Rumus probabilitasnya adalah:
\[ P(Y = y) = \binom{n}{y} p^y (1 - p)^{n - y} \]
Keterangan: - \(n\): jumlah percobaan - \(y\): jumlah sukses - \(p\): probabilitas sukses tiap percobaan
Distribusi binomial sangat sering digunakan dalam pengujian hipotesis proporsi dan estimasi parameter dalam data biner.
# Simulasi 10 sampel Binomial, masing-masing 7 percobaan
binomial_sample <- rbinom(n = 10, size = 7, prob = 0.6)
binomial_sample## [1] 4 2 4 7 4 3 5 5 3 63.3 Distribusi Multinomial
Ketika hasil yang mungkin lebih dari dua kategori, kita menggunakan distribusi multinomial. Distribusi ini merupakan perluasan dari distribusi binomial ke lebih dari dua kategori. Contohnya, ketika kita meneliti warna permen dalam satu bungkus, hasilnya bisa terdiri dari merah, hijau, biru, dan lain-lain. Dalam kasus ini, kita ingin mengetahui berapa banyak dari tiap warna yang muncul dalam sejumlah total permen.
Rumus distribusi multinomial adalah:
\[ P(Y_1 = y_1, ..., Y_k = y_k) = \frac{n!}{y_1! y_2! \dots y_k!} p_1^{y_1} p_2^{y_2} \dots p_k^{y_k} \]
Keterangan: - \(y_i\): jumlah observasi pada kategori ke-\(i\) - \(p_i\): probabilitas kategori ke-\(i\) - \(\sum y_i = n\), \(\sum p_i = 1\)
Distribusi ini sangat berguna untuk memodelkan hasil dari survei, preferensi konsumen, dan klasifikasi.
# Simulasi 1 percobaan Multinomial dengan 3 kategori
set.seed(123)
multinomial_sample <- rmultinom(n = 1, size = 12, prob = c(0.4, 0.35, 0.25))
rownames(multinomial_sample) <- c("Merah", "Biru", "Hijau")
multinomial_sample## [,1]
## Merah 4
## Biru 4
## Hijau 43.4 Distribusi Poisson
Distribusi Poisson digunakan untuk memodelkan jumlah kejadian dalam interval waktu atau ruang tertentu ketika kejadian tersebut bersifat jarang (rare events). Distribusi ini memiliki satu parameter \(\lambda\), yaitu rata-rata jumlah kejadian dalam interval tersebut. Contohnya termasuk jumlah kendaraan yang melewati jalan tol dalam satu jam, atau jumlah pelanggan yang datang ke toko dalam satu hari.
Fungsi probabilitasnya adalah:
\[ P(Y = y) = \frac{e^{-\lambda} \lambda^y}{y!} \quad \text{untuk } y = 0, 1, 2, ... \]
Distribusi ini banyak digunakan dalam epidemiologi, logistik, dan bidang lainnya yang berkaitan dengan kejadian diskrit dalam waktu atau ruang.
# Simulasi 10 kejadian Poisson dengan lambda = 4
poisson_sample <- rpois(10, lambda = 4)
poisson_sample## [1] 3 6 7 1 4 7 4 4 8 44 Desain Sampling dalam Analisis Data Kategori
Dalam penelitian statistik, cara pengambilan sampel memiliki peran besar terhadap kualitas dan validitas hasil analisis, termasuk dalam konteks data kategori. Dua pendekatan utama yang sering digunakan adalah sampling prospektif dan sampling retrospektif, tergantung pada arah waktu pengumpulan data dan tujuan studi.
4.1 Prospective Sampling
Sampling prospektif dilakukan dengan memilih responden terlebih dahulu, lalu mengamati mereka dalam jangka waktu tertentu untuk melihat apakah suatu kejadian akan muncul. Pendekatan ini memberikan keuntungan karena peneliti dapat mengendalikan paparan variabel sebelum hasil diamati. Oleh karena itu, metode ini banyak digunakan dalam penelitian eksperimental dan studi jangka panjang.
4.1.1 Eksperimen
Dalam studi eksperimental, partisipan secara acak dibagi ke dalam kelompok perlakuan dan kontrol. Teknik sampling yang sering digunakan dalam konteks ini meliputi:
Simple Random Sampling (SRS): Semua anggota populasi memiliki peluang yang sama untuk dipilih.
Stratified Random Sampling: Populasi dibagi menjadi beberapa lapisan (strata) berdasarkan karakteristik tertentu, dan sampel diambil secara acak dari masing-masing strata.
Cluster Sampling: Populasi dikelompokkan ke dalam beberapa cluster, lalu beberapa cluster dipilih secara acak untuk diobservasi.
4.1.2 Studi Kohort
Studi kohort bersifat observasional, di mana sekelompok individu dengan karakteristik tertentu diikuti dari waktu ke waktu untuk melihat perkembangan kejadian. Jenis teknik sampling yang bisa digunakan antara lain:
Census Sampling: Seluruh populasi target dijadikan sebagai sampel.
Systematic Sampling: Pemilihan subjek berdasarkan interval tertentu dari daftar populasi.
Matched Sampling: Setiap individu dari satu kelompok dipasangkan dengan individu yang memiliki karakteristik serupa dari kelompok lain.
4.2 Retrospective Sampling
Pada pendekatan ini, data dikumpulkan setelah kejadian sudah terjadi. Responden dikelompokkan berdasarkan outcome (misalnya: sakit atau tidak), kemudian informasi mengenai paparan sebelumnya ditelusuri.
4.2.1 Studi Kasus-Kontrol
Penelitian ini membandingkan kelompok yang mengalami suatu kondisi (kasus) dengan kelompok yang tidak mengalaminya (kontrol). Teknik sampling yang sering digunakan:
Purposive Sampling: Subjek dipilih berdasarkan kriteria yang relevan dengan penelitian.
Snowball Sampling: Responden awal membantu merekomendasikan subjek lain yang memenuhi syarat.
Incidence Density Sampling: Pemilihan kasus dan kontrol berdasarkan waktu kemunculan kejadian.
Stratified dan Cluster Sampling: Bisa diterapkan untuk memastikan distribusi sampel yang representatif.
5 Tabel Kontingensi 2 × 2
Tabel kontingensi 2 × 2 merupakan salah satu alat penting dalam analisis data kategori yang digunakan untuk menampilkan distribusi frekuensi dari dua variabel kategorik. Analisis ini sangat umum dalam studi epidemiologi, survei sosial, dan studi kebijakan kesehatan.
Pada bab ini, akan dianalisis hubungan antara tingkat pendidikan ibu (Tinggi vs Rendah) dengan status imunisasi dasar lengkap anak usia 12–23 bulan (Lengkap vs Tidak Lengkap) berdasarkan data SDKI 2017.
5.1 Data Observasi dan Struktur Tabel
Data SDKI 2017 menunjukkan hasil berikut:
data <- matrix(c(1200, 800, 1000, 1500), nrow = 2, byrow = TRUE)
colnames(data) <- c("Imunisasi Lengkap (+)", "Imunisasi Tidak Lengkap (-)")
rownames(data) <- c("Pendidikan Tinggi", "Pendidikan Rendah")
data## Imunisasi Lengkap (+) Imunisasi Tidak Lengkap (-)
## Pendidikan Tinggi 1200 800
## Pendidikan Rendah 1000 1500Tabel ini menunjukkan bahwa: - Dari ibu berpendidikan tinggi, 1200 anak mendapatkan imunisasi lengkap dan 800 tidak. - Dari ibu berpendidikan rendah, hanya 1000 anak yang mendapatkan imunisasi lengkap, sementara 1500 tidak.
Total anak dalam sampel = 4500 anak usia 12–23 bulan.
5.2 Distribusi Probabilitas
Untuk memahami distribusi data, dihitung tiga jenis peluang: bersama, marginal, dan bersyarat.
5.2.1 Peluang Bersama (Joint Probability)
Peluang bersama menunjukkan probabilitas bahwa dua kejadian terjadi secara bersamaan. Misalnya, peluang seorang ibu berpendidikan tinggi dan anaknya mendapatkan imunisasi lengkap.
\[ P(X = x, Y = y) = \frac{n_{xy}}{n} \]
Keterangan: - \(n_{xy}\): jumlah kasus pada sel dengan kategori \(x\) dan \(y\) - \(n\): total seluruh observasi
total <- sum(data)
P_joint <- data / total
P_joint## Imunisasi Lengkap (+) Imunisasi Tidak Lengkap (-)
## Pendidikan Tinggi 0.2666667 0.1777778
## Pendidikan Rendah 0.2222222 0.33333335.2.2 Peluang Marginal
Peluang marginal menunjukkan probabilitas terjadinya suatu kejadian tanpa memperhatikan kejadian lainnya (di sepanjang baris atau kolom).
Peluang marginal baris (untuk \(X = x\)):
\[ P(X = x) = \frac{n_{x+}}{n} \]
Peluang marginal kolom (untuk \(Y = y\)):
\[ P(Y = y) = \frac{n_{+y}}{n} \] Keterangan: - \(n_{x+}\): jumlah total di baris \(x\) - \(n_{+y}\): jumlah total di kolom \(y\)
tunjukkan distribusi probabilitas masing-masing variabel secara terpisah:
P_marginal_rows <- rowSums(data) / total # Pendidikan
P_marginal_cols <- colSums(data) / total # Status Imunisasi
P_marginal_rows## Pendidikan Tinggi Pendidikan Rendah
## 0.4444444 0.5555556P_marginal_cols## Imunisasi Lengkap (+) Imunisasi Tidak Lengkap (-)
## 0.4888889 0.51111115.2.3 Peluang Bersyarat (Conditional Probability)
Peluang bersyarat mengukur probabilitas suatu kejadian terjadi dengan syarat kejadian lain sudah terjadi.
Peluang \(Y = y\) dengan syarat \(X = x\):
\[ P(Y = y \mid X = x) = \frac{P(X = x, Y = y)}{P(X = x)} = \frac{n_{xy}}{n_{x+}} \]
Peluang \(X = x\) dengan syarat \(Y = y\):
\[ P(X = x \mid Y = y) = \frac{n_{xy}}{n_{+y}} \]
P_conditional <- data / rowSums(data)
P_conditional## Imunisasi Lengkap (+) Imunisasi Tidak Lengkap (-)
## Pendidikan Tinggi 0.6 0.4
## Pendidikan Rendah 0.4 0.65.3 Ukuran Asosiasi
Untuk menilai hubungan antara pendidikan ibu dan status imunisasi, digunakan tiga ukuran asosiasi berikut:
5.3.1 Risk Difference (RD)
Mengukur selisih proporsi risiko:
\[ RD = P(Y|X=1) - P(Y|X=0) \]
RD <- function(n11, n12, n21, n22) {
(n11 / (n11 + n12)) - (n21 / (n21 + n22))
}
RD(1200, 800, 1000, 1500)## [1] 0.25.3.2 Relative Risk (RR)
Membandingkan probabilitas antara dua kelompok:
\[ RR = \frac{P(Y|X=1)}{P(Y|X=0)} \]
RR <- function(n11, n12, n21, n22) {
(n11 / (n11 + n12)) / (n21 / (n21 + n22))
}
RR(1200, 800, 1000, 1500)## [1] 1.55.3.3 Odds Ratio (OR)
Perbandingan odds dua kelompok: \[ OR = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} \]
OR <- function(n11, n12, n21, n22) {
(n11 * n22) / (n12 * n21)
}
OR(1200, 800, 1000, 1500)## [1] 2.255.4 Interpretasi
Dari hasil perhitungan: - Risk Difference (RD) menunjukkan bahwa anak dari ibu berpendidikan tinggi memiliki kelebihan probabilitas imunisasi lengkap sebesar nilai RD dibanding anak dari ibu berpendidikan rendah. - Relative Risk (RR) memberikan rasio risiko imunisasi lengkap antara dua kelompok. - Odds Ratio (OR) menunjukkan bahwa kemungkinan anak imunisasi lengkap lebih besar jika ibunya berpendidikan tinggi dibanding rendah.
Analisis ini memperlihatkan adanya perbedaan mencolok dalam cakupan imunisasi berdasarkan tingkat pendidikan ibu, yang menjadi dasar penting dalam perumusan kebijakan kesehatan berbasis pendidikan.
6 Inferensi Tabel Kontingensi Dua Arah
Bab ini membahas langkah-langkah inferensi dalam tabel kontingensi dua arah untuk mengetahui apakah ada hubungan yang signifikan secara statistik antara dua variabel kategorik. Contoh kasus diambil dari SDKI 2017 mengenai hubungan antara pendidikan ibu dan status imunisasi dasar lengkap anak usia 12–23 bulan.
6.1 Estimasi Parameter
Estimasi digunakan untuk menaksir parameter populasi berdasarkan data sampel. Dalam konteks data kategori (2x2), kita ingin mengestimasi proporsi, dan dari situ menghitung ukuran asosiasi seperti Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR).
6.1.1 Estimasi Titik
Estimasi titik merupakan nilai tunggal yang digunakan untuk memperkirakan parameter populasi berdasarkan data sampel.
Proporsi dasar: \[ \hat{p} = \frac{x}{n} \]
Keterangan: - \(\hat{p}\): estimasi titik proporsi - \(x\): jumlah individu pada kategori yang diinginkan - \(n\): total individu dalam sampel
# Proporsi anak imunisasi lengkap berdasarkan pendidikan ibu
# Estimasi proporsi anak yang imunisasi lengkap berdasarkan pendidikan ibu
p_tinggi <- 1200 / (1200 + 800) # proporsi imunisasi pada ibu berpendidikan tinggi
p_rendah <- 1000 / (1000 + 1500) # proporsi imunisasi pada ibu berpendidikan rendah
p_tinggi## [1] 0.6p_rendah## [1] 0.46.1.2 Estimasi Interval
Estimasi interval memberikan rentang nilai yang kemungkinan besar mengandung nilai parameter populasi sesungguhnya, dengan tingkat kepercayaan tertentu.
CI untuk Proporsi Tunggal \[ \hat{p} \pm Z_{\alpha/2} \cdot \sqrt{ \frac{ \hat{p}(1 - \hat{p}) }{n} } \]
Keterangan: - \(\hat{p}\): proporsi sampel - \(Z_{\alpha/2}\): nilai kritis dari distribusi normal standar - \(n\): ukuran sampel
CI untuk Perbedaan Proporsi (Risk Difference)
\[ RD \pm Z_{\alpha/2} \cdot SE(RD) \]
Keterangan: - \(RD = \hat{p}_1 - \hat{p}_2\) - \(SE(RD) = \sqrt{ \frac{ \hat{p}_1(1 - \hat{p}_1) }{n_1} + \frac{ \hat{p}_2(1 - \hat{p}_2) }{n_2} }\)
CI untuk Log Relative Risk (log RR)
\[ \log(RR) \pm Z_{\alpha/2} \cdot SE\left( \log(RR) \right) \]
Keterangan: - \(RR = \frac{\hat{p}_1}{\hat{p}_2}\) - \(SE\left( \log(RR) \right) = \sqrt{ \frac{1 - \hat{p}_1}{\hat{p}_1 n_1} + \frac{1 - \hat{p}_2}{\hat{p}_2 n_2} }\)
CI untuk Log Odds Ratio (log OR) \[ \log(OR) \pm Z_{\alpha/2} \cdot SE\left( \log(OR) \right) \]
Keterangan: - \(OR = \frac{a \cdot d}{b \cdot c}\) - \(SE\left( \log(OR) \right) = \sqrt{ \frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d} }\)
# Fungsi menghitung CI untuk Risk Difference
prop_diff <- function(a, b, c, d, alpha = 0.05) {
ph <- a / (a + c)
pi <- b / (b + d)
nh <- a + c
ni <- b + d
se_bp <- sqrt((ph * (1 - ph) / nh) + (pi * (1 - pi) / ni))
z_alpha <- qnorm(1 - alpha / 2)
ci_lower <- (ph - pi) - z_alpha * se_bp
ci_upper <- (ph - pi) + z_alpha * se_bp
list(estimate = ph - pi, ci = c(ci_lower, ci_upper))
}
# Contoh estimasi interval dari data SDKI
prop_diff(a = 1200, b = 1000, c = 800, d = 1500)## $estimate
## [1] 0.2
##
## $ci
## [1] 0.1711945 0.22880556.2 Uji Hipotesis
6.2.1 Uji Proporsi Dua Kelompok
Hipotesis yang diuji:
Hipotesis nol:
\[ H_0: p_1 = p_2 \]
Hipotesis alternatif:
\[ H_1: p_1 \ne p_2 \]
Uji proporsi dua kelompok dapat dilakukan menggunakan fungsi
prop.test() dalam R.
# Menggunakan data SDKI
prop.test(x = c(1200, 1000), n = c(2000, 2500)) # n = jumlah total anak untuk masing-masing pendidikan##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(1200, 1000) out of c(2000, 2500)
## X-squared = 177.07, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.1707445 0.2292555
## sample estimates:
## prop 1 prop 2
## 0.6 0.46.2.2 Uji Asosiasi
Uji asosiasi dalam tabel kontingensi 2 × 2 digunakan untuk mengukur hubungan antara dua variabel kategori, terutama dalam konteks data observasional. Tiga ukuran asosiasi yang umum digunakan adalah Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR).
Hipotesis:
- \(H_0\): Tidak ada asosiasi antara dua variabel
- \(H_1\): Terdapat asosiasi antara dua variabel
Risk Difference (RD) \[ RD = \frac{n_{11}}{n_{1.}} - \frac{n_{21}}{n_{2.}} \]
Standard Error (SE):
\[ SE(RD) = \sqrt{ \frac{\hat{p}_1(1 - \hat{p}_1)}{n_{1.}} + \frac{\hat{p}_2(1 - \hat{p}_2)}{n_{2.}} } \]
Statistik uji:
\[ Z_{RD} = \frac{RD}{SE(RD)} \] Relative Risk (RR)
\[ RR = \frac{n_{11} / n_{1.}}{n_{21} / n_{2.}} = \frac{\hat{p}_1}{\hat{p}_2} \]
Standard Error (SE log RR):
\[ SE(\log RR) = \sqrt{ \frac{1}{n_{11}} - \frac{1}{n_{1.}} + \frac{1}{n_{21}} - \frac{1}{n_{2.}} } \]
Statistik uji:
\[ Z_{RR} = \frac{\log(RR)}{SE(\log RR)} \] Odds Ratio (OR) \[ OR = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} \]
Standard Error (SE log OR):
\[ SE(\log OR) = \sqrt{ \frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} } \]
Statistik uji:
\[ Z_{OR} = \frac{\log(OR)}{SE(\log OR)} \] Menggunakan R
n11 <- 1200; n12 <- 800
n21 <- 1000; n22 <- 1500
n1. <- n11 + n12
n2. <- n21 + n22
# Risk Difference
p1 <- n11 / n1.
p2 <- n21 / n2.
rd <- p1 - p2
se_rd <- sqrt((p1 * (1 - p1) / n1.) + (p2 * (1 - p2) / n2.))
z_rd <- rd / se_rd
# Relative Risk
rr <- p1 / p2
se_ln_rr <- sqrt((1 - p1) / n11 + (1 - p2) / n21)
z_rr <- log(rr) / se_ln_rr
# Odds Ratio
or <- (n11 * n22) / (n12 * n21)
se_ln_or <- sqrt(1 / n11 + 1 / n12 + 1 / n21 + 1 / n22)
z_or <- log(or) / se_ln_or
list(RD = rd, Z_RD = z_rd, RR = rr, Z_RR = z_rr, OR = or, Z_OR = z_or)## $RD
## [1] 0.2
##
## $Z_RD
## [1] 13.60828
##
## $RR
## [1] 1.5
##
## $Z_RR
## [1] 13.27196
##
## $OR
## [1] 2.25
##
## $Z_OR
## [1] 13.242436.2.3 Uji Independensi
Uji independensi bertujuan untuk mengetahui apakah terdapat hubungan antara dua variabel kategori, atau apakah keduanya saling bebas secara statistik.
6.2.3.1 Uji Chi-Square
Uji Chi-Square (𝜒²) merupakan salah satu metode paling umum yang digunakan untuk menguji asosiasi antara dua variabel kategorikal dalam tabel kontingensi.
Rumus Statistik Uji: \[ \chi^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \]
Keterangan: - \(O_{ij}\): frekuensi yang diamati pada sel ke-\(i,j\). - \(E_{ij}\): frekuensi harapan, dihitung sebagai:
\[ E_{ij} = \frac{R_i \cdot C_j}{N} \]
dengan: - \(R_i\): total baris ke-\(i\), - \(C_j\): total kolom ke-\(j\), - \(N\): total jumlah seluruh observasi dalam tabel.
Uji ini dilakukan dengan membandingkan nilai statistik 𝜒² yang dihitung dengan nilai kritis dari distribusi Chi-Square dengan derajat kebebasan:
\[ df = (r - 1)(c - 1) \]
Jika nilai p < 0.05, maka kita menolak hipotesis nol dan menyimpulkan bahwa terdapat hubungan antara kedua variabel.
data <- matrix(c(1200, 800, 1000, 1500), nrow = 2, byrow = TRUE)
chisq.test(data, correct = FALSE)##
## Pearson's Chi-squared test
##
## data: data
## X-squared = 177.87, df = 1, p-value < 2.2e-16Uji Likelihood Ratio (G²) \[ G^2 = 2 \sum O_{ij} \log\left(\frac{O_{ij}}{E_{ij}}\right) \]
expected <- chisq.test(data)$expected
G2 <- 2 * sum(data * log(data / expected))
G2## [1] 178.9972qchisq(0.95, df = 1)## [1] 3.8414596.2.3.2 Uji Fisher
6.2.3.2.1 Distribusi Hipergeometrik
Distribusi hipergeometrik digunakan dalam uji Fisher ketika ukuran sampel kecil dan asumsi dari uji Chi-Square tidak terpenuhi. Distribusi ini merepresentasikan probabilitas memilih sejumlah elemen tertentu dari populasi terbatas tanpa pengembalian.
Distribusi hipergeometrik digunakan untuk menghitung probabilitas persis dari konfigurasi tertentu dalam tabel kontingensi 2 × 2.
Rumus Distribusi Hipergeometrik: \[ P(X = x) = \frac{\binom{K}{x} \binom{N - K}{n - x}}{\binom{N}{n}} \]
Keterangan: - \(N\): total populasi - \(K\): total “sukses” dalam populasi - \(n\): ukuran sampel - \(x\): jumlah sukses dalam sampel
Contoh Kasus (Data SDKI)
Tabel kontingensi:
| Imunisasi Lengkap | Tidak Lengkap | |
|---|---|---|
| Pendidikan Tinggi | 1200 | 800 |
| Pendidikan Rendah | 1000 | 1500 |
Uji Fisher dengan Distribusi Hipergeometrik di R
# Membuat tabel data
data <- matrix(c(1200, 800, 1000, 1500), nrow = 2, byrow = TRUE)
# Uji Fisher
fisher.test(data)##
## Fisher's Exact Test for Count Data
##
## data: data
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.991853 2.541810
## sample estimates:
## odds ratio
## 2.2495866.3 Analisis Residual dalam Tabel Kontingensi
6.3.1 Jenis Residual
Pearson Residual: \[ R_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}} \]
Standardized Residual: \[ SR_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}(1 - \frac{r_i}{n})(1 - \frac{c_j}{n})}} \]
Residual digunakan untuk mendeteksi sel mana yang menyumbang paling besar terhadap ketidaksesuaian.
expected <- chisq.test(data)$expected
pearson_resid <- (data - expected) / sqrt(expected)
row_sum <- rowSums(data)
col_sum <- colSums(data)
total_sum <- sum(data)
standardized_resid <- (data - expected) /
sqrt(expected * (1 - row_sum / total_sum) * (1 - col_sum / total_sum))
list(Pearson = pearson_resid, Standardized = standardized_resid)## $Pearson
## [,1] [,2]
## [1,] 7.106691 -6.950480
## [2,] -6.356417 6.216699
##
## $Standardized
## [,1] [,2]
## [1,] 13.33663 -13.04348
## [2,] -13.63636 13.336636.3.2 Deteksi Outlier
Residual besar (positif/negatif) dapat menandakan adanya outlier atau deviasi lokal dari asumsi independensi dalam tabel kontingensi.
# Data
tabel <- matrix(c(1200, 800, 1000, 1500), nrow = 2, byrow = TRUE)
rownames(tabel) <- c("Tinggi", "Rendah")
colnames(tabel) <- c("Lengkap", "Tidak")
# Uji Chi-Square
chi_result <- chisq.test(tabel)
# Tampilkan residual Pearson
chi_result$residuals## Lengkap Tidak
## Tinggi 7.106691 -6.950480
## Rendah -6.356417 6.216699# Menghitung residual terstandarisasi
observed <- tabel
expected <- chi_result$expected
row_sum <- rowSums(observed)
col_sum <- colSums(observed)
total_sum <- sum(observed)
# Matriks residual terstandarisasi
standardized_residual <- (observed - expected) /
sqrt(expected * (1 - row_sum / total_sum) %*% t(1 - col_sum / total_sum))
# Menampilkan hasil lengkap
list(
Observed = observed,
Expected = expected,
Pearson_Residual = chi_result$residuals,
Standardized_Residual = standardized_residual
)## $Observed
## Lengkap Tidak
## Tinggi 1200 800
## Rendah 1000 1500
##
## $Expected
## Lengkap Tidak
## Tinggi 977.7778 1022.222
## Rendah 1222.2222 1277.778
##
## $Pearson_Residual
## Lengkap Tidak
## Tinggi 7.106691 -6.950480
## Rendah -6.356417 6.216699
##
## $Standardized_Residual
## Lengkap Tidak
## Tinggi 13.33663 -13.33663
## Rendah -13.33663 13.336637 Tabel Kontingensi Tiga Arah
7.1 Tabel Parsial dan Marginal
Tabel parsial menunjukkan hubungan antara dua variabel (X dan Y) pada setiap strata dari variabel ketiga (Z). Tabel marginal menyatukan seluruh level Z.
data3 <- array(c(25, 75, 40, 60, 15, 85, 35, 65, 30, 70, 45, 55),
dim = c(2, 2, 3),
dimnames = list(
Gizi = c("Kurang", "Baik"),
Diare = c("Ya", "Tidak"),
Pendidikan = c("Rendah", "Menengah", "Tinggi")
))
# Tabel Parsial
ftable(data3[, , "Rendah"])## Diare Ya Tidak
## Gizi
## Kurang 25 40
## Baik 75 60ftable(data3[, , "Menengah"])## Diare Ya Tidak
## Gizi
## Kurang 15 35
## Baik 85 65ftable(data3[, , "Tinggi"])## Diare Ya Tidak
## Gizi
## Kurang 30 45
## Baik 70 55# Tabel Marginal (Gabungan semua strata Pendidikan)
apply(data3, c(1, 2), sum)## Diare
## Gizi Ya Tidak
## Kurang 70 120
## Baik 230 1807.2 Distribusi Peluang
Peluang Bersama: \[ P(X = x, Y = y, Z = z) = \frac{n_{xyz}}{N} \]
Peluang Marginal: \[ P(X = x) = \sum_{y,z} P(x, y, z) \]
total <- sum(data3)
joint_prob <- data3 / total
ftable(joint_prob)## Pendidikan Rendah Menengah Tinggi
## Gizi Diare
## Kurang Ya 0.04166667 0.02500000 0.05000000
## Tidak 0.06666667 0.05833333 0.07500000
## Baik Ya 0.12500000 0.14166667 0.11666667
## Tidak 0.10000000 0.10833333 0.09166667# Peluang marginal
apply(joint_prob, 1, sum) # Gizi## Kurang Baik
## 0.3166667 0.6833333apply(joint_prob, 2, sum) # Diare## Ya Tidak
## 0.5 0.5apply(joint_prob, 3, sum) # Pendidikan## Rendah Menengah Tinggi
## 0.3333333 0.3333333 0.33333337.3 Tabel Peluang Bersyarat
\[ P(Y = y | X = x, Z = z) = \frac{P(X = x, Y = y, Z = z)}{\sum_{y} P(X = x, Y = y, Z = z)} \]
conditional_probs <- lapply(dimnames(data3)$Pendidikan, function(z) {
prop.table(data3[, , z], 1)
})
names(conditional_probs) <- dimnames(data3)$Pendidikan
conditional_probs## $Rendah
## Diare
## Gizi Ya Tidak
## Kurang 0.3846154 0.6153846
## Baik 0.5555556 0.4444444
##
## $Menengah
## Diare
## Gizi Ya Tidak
## Kurang 0.3000000 0.7000000
## Baik 0.5666667 0.4333333
##
## $Tinggi
## Diare
## Gizi Ya Tidak
## Kurang 0.40 0.60
## Baik 0.56 0.447.4 Ukuran Asosiasi
7.4.1 Rumus:
- \(RD = P(Y=1|X=1) - P(Y=1|X=0)\)
- \(RR = \frac{P(Y=1|X=1)}{P(Y=1|X=0)}\)
- \(OR = \frac{a \cdot d}{b \cdot c}\)
asosiasi_2x2 <- function(a, b, c, d) {
p1 <- a / (a + b)
p2 <- c / (c + d)
rd <- p1 - p2
rr <- p1 / p2
or <- (a * d) / (b * c)
return(data.frame(RD = rd, RR = rr, OR = or))
}
# Pendidikan Rendah
asosiasi_2x2(25, 75, 40, 60)## RD RR OR
## 1 -0.15 0.625 0.5# Menengah
asosiasi_2x2(15, 85, 35, 65)## RD RR OR
## 1 -0.2 0.4285714 0.3277311# Tinggi
asosiasi_2x2(30, 70, 45, 55)## RD RR OR
## 1 -0.15 0.6666667 0.52380957.5 Conditional Independence
Jika dua variabel (X dan Y) tidak berhubungan setelah dikendalikan oleh variabel ketiga (Z), maka keduanya dikatakan independen secara kondisional terhadap Z. Hipotesis: - \(H_0\): X dan Y independen secara kondisional terhadap Z
7.6 Marginal Y dan X
Menjumlahkan seluruh level Z: \[ P(X = x, Y = y) = \sum_z P(X = x, Y = y, Z = z) \]
apply(data3, 1, sum) # Gizi## Kurang Baik
## 190 410apply(data3, 2, sum) # Diare## Ya Tidak
## 300 3007.7 Inferensi Tabel Kontingensi Tiga Arah
7.7.1 Independensi Bersyarat dalam Tabel Kontingensi
Hipotesis: - \(H_0\): Gizi dan Diare independen setelah dikontrol oleh Pendidikan - \(H_1\): Terdapat asosiasi setelah dikontrol oleh Pendidikan
7.7.2 Pengujian Statistik untuk Independensi Bersyarat
mantelhaen.test(data3)##
## Mantel-Haenszel chi-squared test with continuity correction
##
## data: data3
## Mantel-Haenszel X-squared = 18.628, df = 1, p-value = 1.589e-05
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.3166001 0.6454627
## sample estimates:
## common odds ratio
## 0.45205487.7.3 Odds Ratio Bersama
Rataan tertimbang dari OR di semua strata:
\[ OR_{MH} = \frac{\sum_k \frac{a_k d_k}{n_k}}{\sum_k \frac{b_k c_k}{n_k}} \]
# Fungsi custom untuk hitung RD, RR, OR
asosiasi_2x2 <- function(a, b, c, d) {
p1 <- a / (a + b)
p2 <- c / (c + d)
rd <- p1 - p2
rr <- p1 / p2
or <- (a * d) / (b * c)
return(data.frame(RD = rd, RR = rr, OR = or))
}
# Simpan hasil dari tiap strata pendidikan
hasil_odds <- lapply(dimnames(data3)$Pendidikan, function(p) {
tab <- data3[, , p]
a <- tab["Kurang", "Ya"]
b <- tab["Kurang", "Tidak"]
c <- tab["Baik", "Ya"]
d <- tab["Baik", "Tidak"]
hasil <- asosiasi_2x2(a, b, c, d)
hasil$Pendidikan <- p
return(hasil)
})
# Gabungkan semua jadi satu data frame
hasil_odds_df <- do.call(rbind, hasil_odds)
# Tampilkan
print(hasil_odds_df)## RD RR OR Pendidikan
## 1 -0.1709402 0.6923077 0.5000000 Rendah
## 2 -0.2666667 0.5294118 0.3277311 Menengah
## 3 -0.1600000 0.7142857 0.5238095 Tinggi7.7.4 Uji Homogenitas Odds Ratio dengan Statistik Breslow-Day
Hipotesis: - \(H_0\): OR antar strata homogen (tidak ada interaksi) - \(H_1\): OR berbeda (ada efek interaksi)
Statistik uji:
\[ Q_{BD} = \sum_k \frac{(OR_k - OR_{MH})^2}{Var(OR_k)} \]
# install.packages("metafor")
library(metafor)## Warning: package 'metafor' was built under R version 4.3.3## Loading required package: Matrix## Loading required package: metadat## Warning: package 'metadat' was built under R version 4.3.3## Loading required package: numDeriv##
## Loading the 'metafor' package (version 4.8-0). For an
## introduction to the package please type: help(metafor)dat <- escalc(measure = "OR",
ai = c(25, 15, 30),
bi = c(75, 85, 70),
ci = c(40, 35, 45),
di = c(60, 65, 55))
res <- rma(yi, vi, data = dat, method = "FE")
res$QE # Statistik Q (Cochran's Q)## [1] 1.198289res$QEp # p-value## [1] 0.5492814Dengan menggunakan pendekatan ini, kita dapat mengevaluasi apakah hubungan dua variabel tetap konsisten di berbagai level variabel kontrol, atau justru berubah (interaksi). Ini sangat berguna dalam penelitian untuk mendeteksi pengaruh perancu atau efek interaksi yang tersembunyi dalam data agregat. Hal ini memungkinkan kita mengidentifikasi adanya interaksi antar variabel atau efek perancu yang sebelumnya tidak terdeteksi pada analisis dua arah sederhana.
8 Generalized Linear Model (GLM)
8.2 Model Regresi Logistik
Model regresi logistik digunakan untuk kasus di mana variabel respons bersifat biner (0 atau 1). Model ini menghubungkan log-odds (logit) terhadap kombinasi linear prediktor.
Fungsi link logit:
\[ \log\left( \frac{\mu}{1 - \mu} \right) = X\beta \]
Simulasi dan Visualisasi Kurva Logistik
set.seed(2024)
n <- 100
x <- seq(-4, 4, length.out = n)
log_odds <- -1 + 2 * x
prob <- 1 / (1 + exp(-log_odds))
y <- rbinom(n, 1, prob)
data <- data.frame(x = x, y = y, prob = prob)
plot(x, y, pch = 16, col = "gray60",
xlab = "X", ylab = "Y / Probabilitas",
main = "Simulasi Regresi Logistik dengan Kurva Sigmoid")
lines(x, prob, col = "blue", lwd = 2)
abline(h = 0.5, col = "red", lty = 2)
legend("topleft", legend = c("Data Biner", "Kurva Logistik", "Ambang 0.5"),
col = c("gray60", "blue", "red"), pch = c(16, NA, NA), lty = c(NA, 1, 2), bty = "n")
Estimasi Model Regresi Logistik (Simulasi)
set.seed(123)
n <- 150
x <- rnorm(n)
log_odds <- -0.3 + 2 * x
prob <- 1 / (1 + exp(-log_odds))
y <- rbinom(n, 1, prob)
data <- data.frame(x = x, y = y)
model <- glm(y ~ x, data = data, family = binomial)
summary(model)##
## Call:
## glm(formula = y ~ x, family = binomial, data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4437 0.2024 -2.192 0.0284 *
## x 1.6742 0.2940 5.694 1.24e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 203.41 on 149 degrees of freedom
## Residual deviance: 150.08 on 148 degrees of freedom
## AIC: 154.08
##
## Number of Fisher Scoring iterations: 5Interpretasi Koefisien - Intercept menunjukkan log-odds ketika prediktor x = 0. - Koefisien x mengindikasikan perubahan log-odds setiap kenaikan 1 unit x. - Untuk interpretasi yang lebih mudah, transformasi ke odds ratio dapat dilakukan:
exp(coef(model))## (Intercept) x
## 0.6416859 5.3347261Visualisasi Kurva Logit
# Membuat prediksi dari model
data$pred <- predict(model, type = "response")
# Membuat plot
plot(data$x, data$y,
col = "gray40", pch = 16,
xlab = "X (Prediktor)",
ylab = "Probabilitas / Respons",
main = "Kurva Logit pada Regresi Logistik",
cex = 0.7, alpha = 0.5)## Warning in plot.window(...): "alpha" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "alpha" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "alpha" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "alpha" is not a
## graphical parameter## Warning in box(...): "alpha" is not a graphical parameter## Warning in title(...): "alpha" is not a graphical parameter# Menambahkan garis prediksi
lines(data$x, data$pred, col = "blue", lwd = 2)
Evaluasi Model
data$pred_class <- ifelse(data$pred > 0.5, 1, 0)
table(Predicted = data$pred_class, Actual = data$y)## Actual
## Predicted 0 1
## 0 72 24
## 1 16 38Contoh Data Nyata: mtcars
# Menggunakan dataset mtcars
data(mtcars)
df <- mtcars
df$am <- factor(df$am, labels = c("Automatic", "Manual"))
# Model regresi logistik
model_real <- glm(am ~ mpg + hp, data = df, family = binomial)
summary(model_real)##
## Call:
## glm(formula = am ~ mpg + hp, family = binomial, data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -33.60517 15.07672 -2.229 0.0258 *
## mpg 1.25961 0.56747 2.220 0.0264 *
## hp 0.05504 0.02692 2.045 0.0409 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 43.230 on 31 degrees of freedom
## Residual deviance: 19.233 on 29 degrees of freedom
## AIC: 25.233
##
## Number of Fisher Scoring iterations: 7# Odds Ratio
exp(coef(model_real))## (Intercept) mpg hp
## 2.543664e-15 3.524063e+00 1.056588e+00# Prediksi
df$prob <- predict(model_real, type = "response")
df$pred_class <- ifelse(df$prob > 0.5, "Manual", "Automatic")
# Visualisasi menggunakan Base R
# Menentukan warna dan bentuk
colors <- ifelse(df$pred_class == "Manual", "blue", "red")
shapes <- ifelse(df$am == "Manual", 16, 17)
# Plot titik-titik
plot## function (x, y, ...)
## UseMethod("plot")
## <bytecode: 0x0000029afb9addc0>
## <environment: namespace:base>8.3 Model Regresi Poisson
Model regresi Poisson digunakan untuk memodelkan variabel respons berupa data hitung (count data), seperti jumlah kejadian. Fungsi link yang digunakan adalah logaritma:
\[ \log(\mu) = X\beta \]
Simulasi dan Estimasi Model
set.seed(42)
n <- 200
x <- rnorm(n)
lambda <- exp(0.5 + 0.7 * x)
y <- rpois(n, lambda)
data <- data.frame(y = y, x = x)
poisson_model <- glm(y ~ x, data = data, family = poisson)
summary(poisson_model)##
## Call:
## glm(formula = y ~ x, family = poisson, data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.44933 0.06027 7.455 8.99e-14 ***
## x 0.64063 0.05540 11.564 < 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: 381.48 on 199 degrees of freedom
## Residual deviance: 244.45 on 198 degrees of freedom
## AIC: 641.2
##
## Number of Fisher Scoring iterations: 5Visualisasi Hasil Prediksi
plot(data$x, data$y, pch = 16, col = "darkgray", main = "Prediksi Regresi Poisson")
newdata <- data.frame(x = seq(min(x), max(x), length.out = 100))
pred <- predict(poisson_model, newdata = newdata, type = "response")
lines(newdata$x, pred, col = "blue", lwd = 2)
Evaluasi Overdispersion
dispersion <- sum(residuals(poisson_model, type = "pearson")^2) / poisson_model$df.residual
dispersion## [1] 1.101595Jika nilai dispersion > 1, ini mengindikasikan adanya overdispersion, dan model alternatif seperti quasi-Poisson atau negative binomial perlu dipertimbangkan.
Kesimpulan - GLM memperluas model regresi klasik untuk jenis data yang lebih kompleks.
Regresi logistik digunakan untuk data biner, sedangkan regresi Poisson untuk data hitung.
Evaluasi model dapat dilakukan melalui nilai AIC, BIC, dan residual.
Visualisasi kurva logit dan prediksi Poisson membantu interpretasi.
Periksa overdispersion untuk memastikan ketepatan model Poisson.
9 Inferensi GLM
9.1 Ekspektasi dan Varians dalam GLM
Bagian ini menjelaskan bagaimana ekspektasi \(\mathbb{E}(Y)\) dan varians \(\text{Var}(Y)\) bergantung pada distribusi dari GLM. Untuk distribusi:
- Poisson: \(\text{Var}(Y_i) = \mu_i\)
- Binomial: \(\text{Var}(Y_i) = \mu_i(1 - \mu_i)\)
Estimasi parameter \(\hat{\beta}\)
akan bersifat asimptotik normal: \[
\hat{\beta} \sim \mathcal{N}(\beta, \text{Var}(\hat{\beta}))
\]
9.2 Penaksiran Parameter: Newton-Raphson Manual
Metode Newton-Raphson digunakan untuk mencari estimasi parameter \(\hat{\beta}\) yang memaksimumkan fungsi likelihood. Pada GLM, pendekatan ini dilakukan melalui IRLS.
# Simulasi data biner untuk regresi logistik
set.seed(2025)
n <- 120
x <- rnorm(n)
X <- cbind(1, x)
beta_true <- c(-0.8, 1.8)
eta <- X %*% beta_true
p <- 1 / (1 + exp(-eta))
y <- rbinom(n, 1, p)
# Newton-Raphson
beta <- c(0, 0)
tol <- 1e-6
for (i in 1:100) {
eta <- X %*% beta
pi_hat <- 1 / (1 + exp(-eta))
W <- diag(as.numeric(pi_hat * (1 - pi_hat)))
z <- eta + solve(W) %*% (y - pi_hat)
beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
if (sum(abs(beta_new - beta)) < tol) break
beta <- beta_new
}
beta## [,1]
## -0.5825176
## x 1.52238399.3 Diagnostik Model GLM
Diagnostik model diperlukan untuk mengevaluasi kesesuaian model: - Deviance membandingkan model dengan model saturasi. - Residual plot digunakan untuk mendeteksi pola sistematis atau outlier. - Chi-square Pearson dapat digunakan sebagai ukuran goodness-of-fit.
set.seed(1)
x <- rnorm(100)
mu <- exp(0.3 + 0.6 * x)
y <- rpois(100, mu)
model <- glm(y ~ x, family = poisson)
# Residual plot
plot(residuals(model), main = "Residual Plot", ylab = "Residual", pch = 19, col = "darkblue")
abline(h = 0, col = "red", lty = 2)
9.4 Inferensi Regresi Logistik
Model regresi logistik menggunakan fungsi logit untuk memodelkan hubungan antara variabel prediktor dan probabilitas kejadian.
Uji Wald
Uji Wald digunakan untuk menilai apakah koefisien regresi berbeda signifikan dari nol. Rumus: Rumus Uji Wald: \[ Z = \frac{\hat{\beta}}{SE(\hat{\beta})}, \quad \chi^2 = Z^2 \]
set.seed(42)
x <- rnorm(100)
log_odds <- -0.4 + 1.1 * x
p <- 1 / (1 + exp(-log_odds))
y <- rbinom(100, 1, p)
data <- data.frame(x = x, y = y)
model <- glm(y ~ x, data = data, family = binomial)
## Wald Test
beta_hat <- coef(model)["x"]
se_beta <- summary(model)$coefficients["x", "Std. Error"]
z_val <- beta_hat / se_beta
wald_stat <- z_val^2
p_val <- 1 - pchisq(wald_stat, df = 1)
z_val; wald_stat; p_val## x
## 4.412698## x
## 19.4719## x
## 1.020906e-05Uji Likelihood Ratio
Membandingkan model penuh dengan model nol (tanpa prediktor) melalui statistik deviance:
\[ G^2 = 2 (\ell_{full} - \ell_{null}) \]
model_null <- glm(y ~ 1, data = data, family = binomial)
anova(model_null, model, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ 1
## Model 2: y ~ x
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 99 135.37
## 2 98 104.49 1 30.884 2.74e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC dan BIC
Digunakan untuk mengevaluasi efisiensi dan kompleksitas model:
\[ AIC = -2 \ell + 2p \quad \text{dan} \quad BIC = -2 \ell + p \log(n) \]
AIC(model)## [1] 108.4881BIC(model)## [1] 113.69849.5 Inferensi Regresi Poisson dan Estimasi IRLS Manual
Model regresi Poisson diestimasi secara manual menggunakan IRLS.
Estimasi ini kemudian dibandingkan dengan hasil fungsi
glm().
set.seed(2025)
n <- 100
x <- rnorm(n)
X <- cbind(1, x)
beta_true <- c(0.4, 0.9)
eta <- X %*% beta_true
lambda <- exp(eta)
y <- rpois(n, lambda)
# Estimasi IRLS manual
beta <- c(0, 0)
for (i in 1:100) {
eta <- X %*% beta
lambda <- exp(eta)
W <- diag(as.numeric(lambda))
z <- eta + (y - lambda) / lambda
beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
if (sum(abs(beta_new - beta)) < 1e-6) break
beta <- beta_new
}
beta## [,1]
## 0.2340082
## x 1.0372395Perbandingan dengan glm()
data <- data.frame(y = y, x = x)
model_glm <- glm(y ~ x, data = data, family = poisson)
summary(model_glm)##
## Call:
## glm(formula = y ~ x, family = poisson, data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.23401 0.09811 2.385 0.0171 *
## x 1.03724 0.06525 15.898 <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: 364.13 on 99 degrees of freedom
## Residual deviance: 109.19 on 98 degrees of freedom
## AIC: 295.76
##
## Number of Fisher Scoring iterations: 5Uji Wald Regresi Poisson
Setelah model glm() dijalankan, kita hitung nilai Wald
secara manual untuk menguji signifikansi koefisien. Jika p-value <
0.05, maka prediktor signifikan.
coef_val <- coef(model_glm)[2]
se_val <- summary(model_glm)$coefficients[2, 2]
wald_z <- coef_val / se_val
wald_chisq <- wald_z^2
p_value <- 1 - pchisq(wald_chisq, df = 1)
wald_z; wald_chisq; p_value## x
## 15.8976## x
## 252.7337## x
## 0Uji Likelihood Ratio dan AIC/BIC
Perbandingan dilakukan antara model penuh dan model nol. Nilai AIC dan BIC yang lebih kecil menunjukkan model yang lebih baik. Ini juga sejalan dengan prinsip informasi dalam pemodelan statistik.
model_null <- glm(y ~ 1, family = poisson, data = data)
anova(model_null, model_glm, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ 1
## Model 2: y ~ x
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 99 364.13
## 2 98 109.19 1 254.94 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(model_glm)## [1] 295.7613BIC(model_glm)## [1] 300.9717Kesimpulan
- Estimasi parameter pada GLM menggunakan MLE dengan pendekatan numerik.
- Uji Wald dan Likelihood Ratio digunakan untuk pengujian hipotesis.
- Evaluasi model dapat dilakukan melalui AIC, BIC, dan deviance.
- Visualisasi residual membantu deteksi outlier dan penyimpangan.
- Pendekatan manual seperti Newton-Raphson dan IRLS memberikan pemahaman mendalam terhadap proses estimasi.
10 Studi Kasus:
- Tabel Data
# Buat tabel kontingensi
tabel <- matrix(c(647, 622, 2, 27), nrow = 2, byrow = TRUE)
rownames(tabel) <- c("Perokok", "Non-perokok")
colnames(tabel) <- c("Kanker+", "Kanker-")
tabel## Kanker+ Kanker-
## Perokok 647 622
## Non-perokok 2 27Peluang
# Total observasi
total <- sum(tabel)
# Peluang bersama
pj <- tabel / total
# Peluang marginal
marginal_baris <- rowSums(tabel) / total
marginal_kolom <- colSums(tabel) / total
# Peluang bersyarat: P(Kanker | Merokok), P(Kanker | Non-merokok)
p_kanker_perokok <- tabel[1,1] / sum(tabel[1,])
p_kanker_non <- tabel[2,1] / sum(tabel[2,])
# Tampilkan
pj## Kanker+ Kanker-
## Perokok 0.498459168 0.47919877
## Non-perokok 0.001540832 0.02080123marginal_baris## Perokok Non-perokok
## 0.97765794 0.02234206marginal_kolom## Kanker+ Kanker-
## 0.5 0.5p_kanker_perokok## [1] 0.5098503p_kanker_non## [1] 0.06896552ukuran asosiasi
# Risk Difference (RD)
RD <- p_kanker_perokok - p_kanker_non
# Relative Risk (RR)
RR <- p_kanker_perokok / p_kanker_non
# Odds Ratio (OR)
OR <- (tabel[1,1] * tabel[2,2]) / (tabel[1,2] * tabel[2,1])
# Output
cat("Risk Difference (RD):", round(RD, 4), "\n")## Risk Difference (RD): 0.4409cat("Relative Risk (RR):", round(RR, 2), "\n")## Relative Risk (RR): 7.39cat("Odds Ratio (OR):", round(OR, 2), "\n")## Odds Ratio (OR): 14.0411. Pemilihan Model Regresi Logistik dan Evaluasi
11.1 Membangun Model Regresi Logistik: Pendekatan Confirmatory dan Exploratory
Ada dua strategi umum dalam membangun model regresi logistik:
- Confirmatory (berdasarkan teori yang sudah ada).
- Exploratory (berbasis eksplorasi data).
Confirmatory:
• Menggunakan teori atau penelitian sebelumnya sebagai dasar model. • Fokus pada pengujian efek prediktor tertentu (uji signifikansi). • Biasanya dilakukan Likelihood Ratio Test antar model.
Exploratory:
• Digunakan saat belum ada teori pasti. • Tujuan utamanya adalah memilih model terbaik berdasarkan kriteria statistik. • Teknik seleksi: Forward, Backward, dan Stepwise. • Kriteria: AIC, deviance, dan likelihood.
Simulasi Data Baru
{install.packages("DescTools")}
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(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, varlibrary(DescTools)## Warning: package 'DescTools' was built under R version 4.3.3set.seed(456)
n <- 250
x1 <- rnorm(n, mean = 2)
x2 <- rbinom(n, 1, 0.6)
x3 <- rnorm(n, mean = -1, sd = 1.2)
lin_pred <- -0.8 + 0.9 * x1 - 1.1 * x2 + 0.7 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
df <- data.frame(y = as.factor(y), x1, x2, x3)
head(df)## y x1 x2 x3
## 1 0 0.6564786 1 -1.2300001
## 2 1 2.6217756 1 0.2539341
## 3 0 2.8008747 1 -1.5866382
## 4 0 0.6111076 0 -0.7035646
## 5 0 1.2856431 1 -0.4038599
## 6 0 1.6759389 0 -2.519360311.2 Pemilihan Model Stepwise
# Membuat model penuh (full model)
model_full <- glm(y ~ x1 + x2 + x3, data = df, family = binomial)
summary(model_full)##
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.01248 0.37751 -0.033 0.974
## x1 0.67364 0.16093 4.186 2.84e-05 ***
## x2 -1.34397 0.30526 -4.403 1.07e-05 ***
## x3 0.79472 0.14734 5.394 6.90e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 346.56 on 249 degrees of freedom
## Residual deviance: 278.27 on 246 degrees of freedom
## AIC: 286.27
##
## Number of Fisher Scoring iterations: 4# Membuat model null
null_model <- glm(y ~ 1, data = df, family = binomial)
# Stepwise selection (forward, backward, both)
step_forward <- step(null_model, direction = "forward", scope = formula(model_full), trace = FALSE)
step_backward <- step(model_full, direction = "backward", trace = FALSE)
step_both <- step(null_model, direction = "both", scope = formula(model_full), trace = FALSE)
# Menghitung AIC untuk semua model
AIC(model_full, step_forward, step_backward, step_both)## df AIC
## model_full 4 286.2657
## step_forward 4 286.2657
## step_backward 4 286.2657
## step_both 4 286.265711.3 Evaluasi Model: ROC dan AUC
# Mendapatkan probabilitas prediksi dari model stepwise (both)
pred_prob <- predict(step_both, type = "response")
# Membuat objek ROC
roc_obj <- roc(df$y, pred_prob)## Setting levels: control = 0, case = 1## Setting direction: controls < cases# Memplot kurva ROC
plot(roc_obj, main = "Kurva ROC", col = "blue")
# Menghitung dan menampilkan AUC
auc_value <- auc(roc_obj)
print(auc_value)## Area under the curve: 0.778911.4 Evaluasi Pseudo R-Squared
PseudoR2(step_both, which = c("CoxSnell", "Nagelkerke", "McFadden"))## CoxSnell Nagelkerke McFadden
## 0.2390347 0.3187198 0.197057811.5 Tabel Klasifikasi
# Menghitung kelas prediksi (0 atau 1) berdasarkan probabilitas
pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
# Membuat confusion matrix menggunakan table()
conf_matrix <- table(Predicted = pred_class, Actual = df$y)
# Menampilkan confusion matrix
print(conf_matrix)## Actual
## Predicted 0 1
## 0 88 34
## 1 38 90# Menghitung Sensitivity dan Specificity secara manual
# Sensitivity = TP / (TP + FN)
# Specificity = TN / (TN + FP)
TP <- conf_matrix[2, 2]
TN <- conf_matrix[1, 1]
FP <- conf_matrix[1, 2]
FN <- conf_matrix[2, 1]
# Sensitivity dan Specificity
sensitivity <- TP / (TP + FN)
specificity <- TN / (TN + FP)
# Menampilkan hasil Sensitivity dan Specificity
cat("Sensitivity:", sensitivity, "\n")## Sensitivity: 0.703125cat("Specificity:", specificity, "\n")## Specificity: 0.7213115Model terbaik berdasarkan AIC:
# Menghitung AIC untuk masing-masing model
aic_values <- c(AIC(model_full), AIC(step_forward), AIC(step_backward), AIC(step_both))
# Menyimpan nama model dan mencari model dengan AIC terkecil
model_names <- c("model_full", "step_forward", "step_backward", "step_both")
best_model <- model_names[which.min(aic_values)]
# Menampilkan nama model dengan AIC terkecil
print(best_model)## [1] "model_full"# Menampilkan AIC untuk setiap model
aic_values <- c(AIC(model_full), AIC(step_forward), AIC(step_backward), AIC(step_both))
model_names <- c("model_full", "step_forward", "step_backward", "step_both")
# Menampilkan AIC untuk masing-masing model
data.frame(Model = model_names, AIC = aic_values)## Model AIC
## 1 model_full 286.2657
## 2 step_forward 286.2657
## 3 step_backward 286.2657
## 4 step_both 286.2657Nilai AUC:
round(auc(roc_obj), 2)## [1] 0.7811.6 Perbandingan Model Regresi Logistik (Data Baru)
set.seed(789) # seed untuk data baru
n <- 300
x1 <- rnorm(n, mean = 1)
x2 <- rbinom(n, 1, 0.4)
x3 <- rnorm(n, mean = 0.5, sd = 0.9)
# Membuat model linier
lin_pred <- -1.2 + 1.0 * x1 - 0.7 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred)) # Probabilitas prediksi
y <- rbinom(n, 1, p) # Variabel target
data <- data.frame(y = as.factor(y), x1, x2, x3)
# Membuat model GLM
model1 <- glm(y ~ x1, data = data, family = binomial)
model2 <- glm(y ~ x1 + x2, data = data, family = binomial)
model3 <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
# Membandingkan AIC dan deviance untuk setiap model
model_comp <- data.frame(
Model = c("Model 1", "Model 2", "Model 3"),
AIC = c(AIC(model1), AIC(model2), AIC(model3)),
Deviance = c(deviance(model1), deviance(model2), deviance(model3))
)
# Menampilkan hasil perbandingan model
print(model_comp)## Model AIC Deviance
## 1 Model 1 372.8852 368.8852
## 2 Model 2 370.8967 364.8967
## 3 Model 3 362.5519 354.5519Prinsip Parsimony: model terbaik tidak selalu yang paling kompleks, tetapi yang efisien dan tetap akurat.
11.7 Likelihood-Ratio Test
LRT digunakan untuk menguji apakah penambahan variabel signifikan meningkatkan performa model.
#### Bandingkan model dengan dan tanpa x2:
anova(model1, model2, test = "LRT")## Analysis of Deviance Table
##
## Model 1: y ~ x1
## Model 2: y ~ x1 + x2
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 298 368.89
## 2 297 364.90 1 3.9885 0.04581 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### Bandingkan model dengan dan tanpa x3:
anova(model2, model3, test = "LRT")## Analysis of Deviance Table
##
## Model 1: y ~ x1 + x2
## Model 2: y ~ x1 + x2 + x3
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 297 364.90
## 2 296 354.55 1 10.345 0.001298 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretasi:
Karena p-value < 0.05, kita dapat menyimpulkan bahwa penambahan variabel x3 (yang sebelumnya tidak ada dalam Model 1) secara signifikan memperbaiki model. Artinya, variabel x3 memberikan kontribusi yang penting dalam meningkatkan kemampuan model dalam memprediksi variabel dependen y.
Dengan demikian, kita dapat mengatakan bahwa Model 2 yang mencakup x1, x2, dan x3 adalah model yang lebih baik dibandingkan Model 1, karena memberikan penurunan deviance yang signifikan dan meningkatkan kebaikan model secara keseluruhan.
11.8 Prinsip Parsimony dan Rumus
Parsimony = memilih model paling sederhana yang tetap memiliki performa baik.
Rumus AIC:
AIC = -2(log(L) - k) = -2 * log-likelihood + 2 * jumlah parameter
Rumus Deviance:
D = -2 * (log L(model) - log L(model saturasi))
Rumus Likelihood-Ratio:
G^2 = -2 * (log L0 - log L1)
Catatan: • AIC penalti kompleksitas, jadi model sederhana diutamakan jika performanya mirip. • Deviance kecil artinya prediksi model mendekati data aktual.
11.9 Evaluasi Akurasi Model dan Tabel Klasifikasi
# Menghasilkan probabilitas prediksi dari model3
pred_prob <- predict(model3, type = "response")
# Menentukan kelas prediksi (0 atau 1) berdasarkan probabilitas
pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
# Membuat confusion matrix menggunakan table()
conf_matrix <- table(Predicted = pred_class, Actual = data$y)
# Menampilkan confusion matrix
print(conf_matrix)## Actual
## Predicted 0 1
## 0 126 52
## 1 42 80# Menghitung Sensitivity dan Specificity secara manual
# Sensitivity = TP / (TP + FN)
# Specificity = TN / (TN + FP)
TP <- conf_matrix[2, 2] # True Positives (1,1)
TN <- conf_matrix[1, 1] # True Negatives (0,0)
FP <- conf_matrix[1, 2] # False Positives (0,1)
FN <- conf_matrix[2, 1] # False Negatives (1,0)
# Sensitivity dan Specificity
sensitivity <- TP / (TP + FN)
specificity <- TN / (TN + FP)
# Menampilkan hasil Sensitivity dan Specificity
cat("Sensitivity:", sensitivity, "\n")## Sensitivity: 0.6557377cat("Specificity:", specificity, "\n")## Specificity: 0.7078652Interpretasi:
• Sensitivitas: kemampuan model deteksi kelas positif. • Spesifisitas: kemampuan model deteksi kelas negatif.
11.10 Precision-Recall Curve (PR Curve)
Kurva Precision-Recall sangat berguna terutama ketika data tidak seimbang antara kelas 1 dan 0. Fokus utamanya pada prediksi kelas positif.
Definisi:
Precision = TP / (TP + FP) Recall (Sensitivitas) = TP / (TP + FN)
# Install dan load package:
if (!require(PRROC)) install.packages("PRROC"); library(PRROC)## Loading required package: PRROC## Warning: package 'PRROC' was built under R version 4.3.3## Loading required package: rlang## Warning: package 'rlang' was built under R version 4.3.3# Simulasi data baru:
set.seed(123)
x1 <- rnorm(200)
x2 <- rbinom(200, 1, 0.5)
x3 <- rnorm(200)
lin_pred <- -1 + 1.5 * x1 - 0.7 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(200, 1, p)
data <- data.frame(y = y, x1, x2, x3)
# Model logistik
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
prob <- predict(model, type = "response")
# Plot PR Curve
pr <- pr.curve(scores.class0 = prob[data$y == 1],
scores.class1 = prob[data$y == 0],
curve = TRUE)
plot(pr)
Interpretasi:
- PR Curve sangat berguna untuk data dengan ketidakseimbangan kelas.
- AUPRC (luas area di bawah kurva) mendekati 1 artinya model sangat baik.
11.11 Pseudo R-squared pada Regresi Logistik
Pseudo R² digunakan untuk menilai kekuatan prediksi model logistik, karena R² klasik tidak berlaku.
# Simulasi data
set.seed(123)
n <- 300
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)
lin_pred <- -1 + 1.2 * x1 - 0.6 * x2 + 0.8 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)
# Model logistik
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
model_null <- glm(y ~ 1, data = data, family = binomial)
# Perhitungan manual
logL0 <- logLik(model_null)
logLM <- logLik(model)
L0 <- exp(logL0)
LM <- exp(logLM)
n <- nobs(model)
cox_snell <- 1 - (L0 / LM)^(2 / n)
mcfadden <- 1 - (as.numeric(logLM) / as.numeric(logL0))
r2 <- data.frame(
R2_Cox_Snell = cox_snell,
R2_McFadden = mcfadden
)
r2## R2_Cox_Snell R2_McFadden
## 1 0.2698462 0.258629211.11.1 Perhitungan Otomatis dengan Paket Tambahan
# pscl
if (!require(pscl)) install.packages("pscl"); library(pscl)## Loading required package: pscl## Warning: package 'pscl' was built under R version 4.3.3## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.pR2(model)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML r2CU
## -135.2290328 -182.4040393 94.3500130 0.2586292 0.2698462 0.3835251# rcompanion
# Misalkan Anda sudah memiliki model regresi logistik (model)
# model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
# Deviance Model (model yang sudah disesuaikan)
deviance_model <- deviance(model)
# Deviance Null (model tanpa prediktor, hanya intercept)
deviance_null <- deviance(glm(y ~ 1, data = data, family = binomial))
# df untuk model dan model null
df_model <- df.residual(model)
df_null <- df.residual(glm(y ~ 1, data = data, family = binomial))
# Menghitung Nagelkerke R²
nagelkerke_r2 <- 1 - (deviance_model / deviance_null)^(df_model / df_null)
# Menampilkan hasil Nagelkerke R²
cat("Nagelkerke R²:", nagelkerke_r2, "\n")## Nagelkerke R²: 0.2563998# DescTools
if (!require(DescTools)) install.packages("DescTools"); library(DescTools)
PseudoR2(model, which = "all")## McFadden McFaddenAdj CoxSnell Nagelkerke AldrichNelson
## 0.2586292 0.2366998 0.2698462 0.3835251 0.2392545
## VeallZimmermann Efron McKelveyZavoina Tjur AIC
## 0.4360055 0.2893849 0.4315315 0.2936202 278.4580657
## BIC logLik logLik0 G2
## 293.2731956 -135.2290328 -182.4040393 94.3500130Interpretasi: • McFadden R² > 0.2 dianggap cukup baik. • Cox-Snell R² konservatif, Nagelkerke merupakan versi yang distandarisasi dari Cox-Snell. • Semakin tinggi nilai pseudo-R², semakin baik performa prediktif model.
12. Apa itu Distribusi Multinomial
Distribusi multinomial adalah generalisasi dari distribusi binomial untuk respon yang memiliki lebih dari dua kategori. Misal X = (X1, X2, …, Xk), maka: P(X1 = x1, …, Xk = xk) = n! / (x1! * x2! * … * xk!) * p1^x1 * p2^x2 * … * pk^xk dengan syarat: ∑x = n dan ∑p = 1
12.1 Studi Kasus: Pemilihan Minuman Favorit
Misal : Teh = 5, Kopi = 3, Jus = 2 p_teh = 0.4, p_kopi = 0.3, p_jus = 0.3
n <- 10
x <- c(5, 3, 2)
p <- c(0.4, 0.3, 0.3)
# Hitung peluang dengan rumus multinomial
peluang <- factorial(n) / prod(factorial(x)) * prod(p^x)
peluang## [1] 0.0627056612.2 Multinomial Logistic Regression
Model regresi logistik untuk respon > 2 kategori. Menggunakan kategori baseline (default di R: kategori terakhir) untuk membandingkan logit tiap kategori lainnya.
12.3 Contoh Kasus Baru
Sebuah perusahaan catering ingin mengetahui preferensi karyawan terhadap 3 jenis makanan utama yang mereka sediakan: Nasi Kotak, Burger, dan Salad. Mereka juga mengumpulkan informasi tentang umur, divisi kerja (Produksi, Admin, Sales), dan jumlah tahun pengalaman kerja. Perusahaan catering mensurvei 150 karyawan dan mencatat: • FoodChoice: Jenis makanan yang dipilih (NasiKotok, Burger, Salad) • Age: Usia karyawan • Department: Produksi, Admin, Sales • Experience: Lama pengalaman kerja (tahun) Tujuannya: Mengetahui faktor-faktor yang memengaruhi preferensi jenis makanan.
12.4 Studi kasus: Preferensi makanan karyawan
set.seed(789)
n <- 150
Department <- sample(c("Produksi", "Admin", "Sales"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 30, sd = 5))
Experience <- round(pmax(rnorm(n, mean = 5, sd = 2), 0))
FoodChoice <- sapply(Department, function(dep) {
if (dep == "Produksi") {
sample(c("NasiKotok", "Burger", "Salad"), 1, prob = c(0.5, 0.3, 0.2))
} else if (dep == "Admin") {
sample(c("NasiKotok", "Burger", "Salad"), 1, prob = c(0.3, 0.5, 0.2))
} else {
sample(c("NasiKotok", "Burger", "Salad"), 1, prob = c(0.2, 0.3, 0.5))
}
})
df <- data.frame(FoodChoice = factor(FoodChoice), Age, Department = factor(Department), Experience)
df$FoodChoice <- relevel(df$FoodChoice, ref = "NasiKotok") # baseline
head(df)## FoodChoice Age Department Experience
## 1 NasiKotok 31 Produksi 3
## 2 Burger 31 Admin 2
## 3 Salad 33 Admin 7
## 4 Burger 34 Admin 3
## 5 Salad 25 Sales 3
## 6 NasiKotok 33 Produksi 612.5 Estimasi Model Multinomial Logistik
if (!require(nnet)) install.packages("nnet"); library(nnet)## Loading required package: nnet## Warning: package 'nnet' was built under R version 4.3.3model_mnlogit <- multinom(FoodChoice ~ Age + Department + Experience, data = df)## # weights: 18 (10 variable)
## initial value 164.791843
## iter 10 value 155.209955
## final value 155.187984
## convergedsummary(model_mnlogit)## Call:
## multinom(formula = FoodChoice ~ Age + Department + Experience,
## data = df)
##
## Coefficients:
## (Intercept) Age DepartmentProduksi DepartmentSales Experience
## Burger 0.59585331 0.01505683 -0.7398530 -0.228914 -0.09701668
## Salad 0.07406986 -0.01791812 -0.4930769 1.075144 0.06882025
##
## Std. Errors:
## (Intercept) Age DepartmentProduksi DepartmentSales Experience
## Burger 1.351983 0.04063841 0.4748955 0.5376744 0.1029401
## Salad 1.423995 0.04231836 0.5401206 0.5373518 0.1079746
##
## Residual Deviance: 310.376
## AIC: 330.37612.6 P-Value dan Interpretasi Koefisien
z <- summary(model_mnlogit)$coefficients / summary(model_mnlogit)$standard.errors
pval <- 2 * (1 - pnorm(abs(z)))
round(pval, 4)## (Intercept) Age DepartmentProduksi DepartmentSales Experience
## Burger 0.6594 0.711 0.1193 0.6703 0.3460
## Salad 0.9585 0.672 0.3613 0.0454 0.5239Interpretasi: • Koefisien untuk kategori “Burger” dan “Salad” dibandingkan dengan baseline “NasiKotok” • p-value < 0.05 menunjukkan pengaruh signifikan terhadap preferensi makanan
12.7 Prediksi dan Validasi
# Lakukan prediksi kategori pada data training
df$Predicted <- predict(model_mnlogit)
# Tabel perbandingan antara prediksi vs data aktual
table(Predicted = df$Predicted, Actual = df$FoodChoice)## Actual
## Predicted NasiKotok Burger Salad
## NasiKotok 14 15 11
## Burger 19 28 14
## Salad 11 12 26Interpretasi: Tabel ini menunjukkan akurasi model dalam memprediksi pilihan makanan.
12.8 Kesimpulan
Kesimpulan dari model multinomial logistik: • Model mampu menjelaskan hubungan antara atribut karyawan dan preferensi makanan • Beberapa faktor terbukti signifikan memengaruhi pilihan • Model dapat digunakan untuk memprediksi preferensi makanan karyawan baru berdasarkan usia, pengalaman, dan divisi
Studi kasus ini menunjukkan implementasi nyata penggunaan regresi logistik multinomial.
12.8.1 Contoh Kasus 2 di R (Dataset Iris)
# Kita gunakan dataset iris untuk memprediksi spesies bunga berdasarkan petal
data(iris)
library(nnet)
library(dplyr)
iris <- iris %>% mutate(Species = relevel(Species, ref = "setosa"))
# Model multinomial logistik
model <- multinom(Species ~ Petal.Length + Petal.Width, data = iris)## # weights: 12 (6 variable)
## initial value 164.791843
## iter 10 value 12.657828
## iter 20 value 10.374056
## iter 30 value 10.330881
## iter 40 value 10.306926
## iter 50 value 10.300057
## iter 60 value 10.296452
## iter 70 value 10.294046
## iter 80 value 10.292029
## iter 90 value 10.291154
## iter 100 value 10.289505
## final value 10.289505
## stopped after 100 iterationssummary(model)## Call:
## multinom(formula = Species ~ Petal.Length + Petal.Width, data = iris)
##
## Coefficients:
## (Intercept) Petal.Length Petal.Width
## versicolor -22.79944 6.92122 7.878496
## virginica -67.82521 12.64721 18.261016
##
## Std. Errors:
## (Intercept) Petal.Length Petal.Width
## versicolor 44.3859 37.58715 81.00888
## virginica 46.3939 37.65702 81.09482
##
## Residual Deviance: 20.57901
## AIC: 32.57901Interpretasi koefisien dengan nilai p-value
z <- summary(model)$coefficients / summary(model)$standard.errors
p_values <- 2 * (1 - pnorm(abs(z)))
round(p_values, 4)## (Intercept) Petal.Length Petal.Width
## versicolor 0.6075 0.8539 0.9225
## virginica 0.1438 0.7370 0.8218Nilai p menunjukkan apakah Petal.Length dan Petal.Width signifikan dalam memprediksi spesies Prediksi spesies dan visualisasi
# Model menggunakan multinom dari paket nnet
install.packages("nnet") ## Warning: package 'nnet' is in use and will not be installedlibrary(nnet)
# Membuat model multinomial
model <- multinom(Species ~ Petal.Length + Petal.Width, data = iris)## # weights: 12 (6 variable)
## initial value 164.791843
## iter 10 value 12.657828
## iter 20 value 10.374056
## iter 30 value 10.330881
## iter 40 value 10.306926
## iter 50 value 10.300057
## iter 60 value 10.296452
## iter 70 value 10.294046
## iter 80 value 10.292029
## iter 90 value 10.291154
## iter 100 value 10.289505
## final value 10.289505
## stopped after 100 iterations# Prediksi dengan model yang telah dilatih
iris$predicted <- predict(model, newdata = iris)
# Menampilkan tabel perbandingan prediksi dan data aktual
table(Predicted = iris$predicted, Actual = iris$Species)## Actual
## Predicted setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 47 3
## virginica 0 3 47# Visualisasi hasil prediksi menggunakan plot dasar
plot(iris$Petal.Length, iris$Petal.Width, col = iris$predicted,
xlab = "Petal Length", ylab = "Petal Width", main = "Prediksi Multinomial Logistic Regression")
legend("topright", legend = levels(iris$Species), fill = 1:length(levels(iris$Species)))
13 Regresi Logistik Ordinal
Regresi logistik ordinal digunakan ketika variabel respon memiliki urutan alami (ordinal), contohnya tingkat kepuasan, level pendidikan, atau tingkatan risiko.
Model ini berbeda dengan: • Regresi logistik biner: hanya dua level • Regresi logistik multinomial: >2 kategori tapi tidak ada urutan
13.1 Konsep Cumulative Logit Model
Model logit kumulatif membentuk log-odds kumulatif:
log(P(Y ≤ j) / P(Y > j)) = α_j + βx Di mana: - α_j adalah intercept untuk setiap cutoff kategori j - β adalah koefisien regresi, dan diasumsikan sama untuk seluruh batas kategori (proportional odds) Untuk c kategori, akan ada (c-1) fungsi logit.
13.2 Interpretasi Koefisien
Nilai koefisien β memberi informasi arah pengaruh prediktor:
• β > 0 → x meningkatkan kemungkinan masuk kategori rendah atau sedang • β < 0 → x meningkatkan kemungkinan masuk ke kategori tinggi
Odds ratio (OR) = exp(β) memberikan besarnya pengaruh relatif.
13.3 Contoh Data: Penilaian Risiko Kesehatan
Simulasi data:
responden menilai risiko kesehatannya sebagai Rendah, Sedang, atau Tinggi,berdasarkan skor pola hidup sehat (semakin tinggi, semakin sehat).
set.seed(789)
n <- 200
health_score <- round(runif(n, 1, 10))
risk_level <- cut(6 - 0.5 * health_score + rnorm(n),
breaks = c(-Inf, 2.5, 4.5, Inf),
labels = c("Tinggi", "Sedang", "Rendah"),
ordered_result = TRUE
)
df <- data.frame(risk_level, health_score)
head(df)## risk_level health_score
## 1 Tinggi 7
## 2 Rendah 2
## 3 Rendah 1
## 4 Sedang 6
## 5 Sedang 5
## 6 Rendah 113.4 Estimasi Model Ordinal
Gunakan fungsi polr() dari package MASS
if (!require(MASS)) install.packages("MASS"); library(MASS)
model_ord <- polr(risk_level ~ health_score, data = df, Hess = TRUE)
summary(model_ord)## Call:
## polr(formula = risk_level ~ health_score, data = df, Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## health_score -0.8505 0.09018 -9.431
##
## Intercepts:
## Value Std. Error t value
## Tinggi|Sedang -5.7052 0.6036 -9.4524
## Sedang|Rendah -2.5209 0.3947 -6.3873
##
## Residual Deviance: 285.9887
## AIC: 291.988713.4 Estimasi Model Ordinal
Gunakan fungsi polr() dari package MASS
if (!require(MASS)) install.packages("MASS"); library(MASS)
model_ord <- polr(risk_level ~ health_score, data = df, Hess = TRUE)
summary(model_ord)## Call:
## polr(formula = risk_level ~ health_score, data = df, Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## health_score -0.8505 0.09018 -9.431
##
## Intercepts:
## Value Std. Error t value
## Tinggi|Sedang -5.7052 0.6036 -9.4524
## Sedang|Rendah -2.5209 0.3947 -6.3873
##
## Residual Deviance: 285.9887
## AIC: 291.988713.5 Nilai P-Value
Menghitung p-value dari output t-value
(ctable <- coef(summary(model_ord)))## Value Std. Error t value
## health_score -0.8504962 0.09018067 -9.431026
## Tinggi|Sedang -5.7051838 0.60356805 -9.452428
## Sedang|Rendah -2.5209265 0.39467968 -6.387272p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = round(p, 4)))## Value Std. Error t value p value
## health_score -0.8504962 0.09018067 -9.431026 0
## Tinggi|Sedang -5.7051838 0.60356805 -9.452428 0
## Sedang|Rendah -2.5209265 0.39467968 -6.387272 0Interpretasi:
p-value < 0.05 menunjukkan bahwa prediktor signifikan terhadap perubahan tingkat risiko
13.6 Prediksi Probabilitas
Prediksi probabilitas risiko berdasarkan skor pola hidup sehat
newdata <- data.frame(health_score = 4:8)
predict(model_ord, newdata = newdata, type = "probs")## Tinggi Sedang Rendah
## 1 0.09085838 0.6161831 0.29295848
## 2 0.18958598 0.6600252 0.15038882
## 3 0.35383908 0.5758584 0.07030256
## 4 0.56175554 0.4069508 0.03129362
## 5 0.75003259 0.2363547 0.0136127313.7 Goodness-of-Fit dan Proportional Odds
Model cumulative logit mengasumsikan proportional odds: efek prediktor sama untuk tiap level cut-off. Jika asumsi ini tidak terpenuhi, maka model bisa dipertimbangkan dengan versi lebih fleksibel: • Generalized Ordered Logit Model • Partial Proportional Odds Model
Analisis goodness-of-fit dapat dilakukan dengan membandingkan deviance model dan log-likelihood.
13.8 Alternatif Model Ordinal
Jika asumsi proportional odds tidak terpenuhi, kita bisa mempertimbangkan model ordinal alternatif berikut:
- Adjacent-Category Logit Model:
• Model ini membandingkan peluang dua kategori yang bersebelahan (adjacent). • Cocok untuk data ordinal di mana perpindahan antar kategori terjadi satu per satu. • Bentuk umum: log(P(Y = j) / P(Y = j+1)) = α_j + βx
- Continuation-Ratio (Sequential) Logit Model:
• Model ini melihat peluang “melanjutkan” ke level kategori lebih tinggi. • Cocok digunakan dalam proses bertahap, seperti jenjang pendidikan atau tahapan rekrutmen. • Bentuk umum: log(P(Y = j | Y ≥ j)) = α_j + βx
Kedua model ini tidak memerlukan asumsi efek yang konsisten antar batas kategori, sehingga menjadi alternatif ketika asumsi proportional odds pada cumulative logit dilanggar.
Catatan: - Adjacent-category logit cocok jika kategori dinilai dalam pasangan langsung (pairwise). - Continuation-ratio logit cocok untuk proses yang bersifat progresif atau bertahap.
13.9 Kesimpulan
Kesimpulan dari regresi logistik ordinal: • Cocok digunakan untuk variabel respon yang memiliki urutan alami (ordinal). • Model cumulative logit membantu menganalisis hubungan antara prediktor dan peluang kumulatif kategori. • Estimasi dilakukan menggunakan fungsi polr() dari package MASS. • Untuk evaluasi lebih lanjut, bisa digunakan: - Uji deviance (model fit) - Likelihood Ratio Test antar model
13.10 Asumsi Paralelisme dalam Regresi Logistik Ordinal
Model cumulative logit mengasumsikan efek prediktor (β) tetap sama untuk semua level kumulatif respon. Ini disebut sebagai asumsi paralelisme atau proportional odds assumption.
Bentuk umum model:
log(P(Y ≤ j) / P(Y > j)) = α_j + βx, untuk j = 1, …, c - 1 - Hanya intercept (α_j) yang berubah. - Koefisien β tetap sama di seluruh fungsi logit kumulatif.
Visualisasi: Garis logit terhadap x harus sejajar antar level kategori (kemiringan sama).
Jika asumsi ini dilanggar: • Efek prediktor bisa berbeda di tiap batas kategori • Model cumulative logit tidak valid lagi • Solusi: Gunakan model alternatif seperti: - Generalized Ordinal Logistic Regression - Partial Proportional Odds Model
Bab 14 Log Linear Model
14.1 Tabel Kontingensi dan Model Loglinier
Dalam analisis data kategorik, tabel kontingensi merupakan langkah awal untuk mengeksplorasi hubungan antara dua atau lebih variabel kategorik. Misalnya, sebuah studi ingin mengetahui apakah pelatihan kerja berpengaruh terhadap status pekerjaan seseorang.
# Tabel 2x2: Pelatihan vs Status Bekerja
table_data <- matrix(c(28, 22, 40, 60), nrow = 2,
dimnames = list(Pelatihan = c("Ya", "Tidak"),
Bekerja = c("Ya", "Tidak")))
table_data## Bekerja
## Pelatihan Ya Tidak
## Ya 28 40
## Tidak 22 60Tabel ini menyajikan frekuensi gabungan dari dua kategori. Meski bersifat deskriptif, ia bisa memberikan gambaran awal tentang asosiasi antara variabel.
Model loglinier digunakan untuk memodelkan log dari ekspektasi frekuensi ini, menggunakan pendekatan Generalized Linear Model dengan distribusi Poisson.
14.2 Model Saturated
Model saturated (model penuh) menyertakan semua efek utama dan interaksi antar variabel. Model ini cocok sempurna dengan data, artinya tidak ada sisa deviance.
library(MASS)
data_sat <- matrix(c(38, 62, 42, 58), nrow = 2, byrow = TRUE)
dimnames(data_sat) <- list(Pelatihan = c("Ya", "Tidak"),
Bekerja = c("Ya", "Tidak"))
ftable(data_sat)## Bekerja Ya Tidak
## Pelatihan
## Ya 38 62
## Tidak 42 58model_saturated <- loglm(~ Pelatihan * Bekerja, data = data_sat)
summary(model_saturated)## Formula:
## ~Pelatihan * Bekerja
## attr(,"variables")
## list(Pelatihan, Bekerja)
## attr(,"factors")
## Pelatihan Bekerja Pelatihan:Bekerja
## Pelatihan 1 0 1
## Bekerja 0 1 1
## attr(,"term.labels")
## [1] "Pelatihan" "Bekerja" "Pelatihan:Bekerja"
## 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 1Karena ini model penuh, nilai deviance-nya nol dan tidak bisa digunakan untuk menguji kecocokan model.
14.3 Model Independent
Model independen berasumsi bahwa tidak ada hubungan antara dua variabel. Dengan kata lain, status bekerja seseorang tidak tergantung pada apakah dia mengikuti pelatihan atau tidak.
model_indep <- loglm(~ Pelatihan + Bekerja, data = data_sat)
summary(model_indep)## Formula:
## ~Pelatihan + Bekerja
## attr(,"variables")
## list(Pelatihan, Bekerja)
## attr(,"factors")
## Pelatihan Bekerja
## Pelatihan 1 0
## Bekerja 0 1
## attr(,"term.labels")
## [1] "Pelatihan" "Bekerja"
## 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 0.3334415 1 0.5636396
## Pearson 0.3333333 1 0.5637029Nilai deviance dari model ini bisa dibandingkan dengan model saturated untuk menguji apakah asumsi independensi layak.
14.4 Odds Ratio dan Interpretasi
Odds ratio (OR) merupakan ukuran asosiasi pada tabel kontingensi 2x2. OR menghitung seberapa besar peluang suatu kejadian terjadi pada satu kelompok dibandingkan dengan kelompok lainnya.
OR <- (data_sat[1,1] * data_sat[2,2]) / (data_sat[1,2] * data_sat[2,1])
logOR <- log(OR)
logOR## [1] -0.1667748Interpretasi:
- OR = 1 → tidak ada hubungan
- OR > 1 → hubungan positif
- OR < 1 → hubungan negatif
14.5 Estimasi Parameter
Dalam model saturated, parameter diperkirakan dengan metode Iterative Proportional Fitting (IPF). Estimasi log-odds ratio dihitung berdasarkan proporsi data aktual.
# Sudah dihitung di bagian 14.4 (logOR)Estimasi logOR memberi informasi arah dan kekuatan hubungan antar variabel.
14.6 Model Lebih Sederhana dan Perbandingan Model
Perbandingan antara model dilakukan dengan uji deviance (likelihood ratio test). Tujuannya adalah melihat apakah model sederhana cukup menjelaskan data atau perlu model yang lebih kompleks.
anova(model_indep, model_saturated)## LR tests for hierarchical log-linear models
##
## Model 1:
## ~Pelatihan + Bekerja
## Model 2:
## ~Pelatihan * Bekerja
##
## Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1 0.3334415 1
## Model 2 0.0000000 0 0.3334415 1 0.56364
## Saturated 0.0000000 0 0.0000000 0 1.00000Jika nilai p dari deviance signifikan (p < 0.05), maka model independen tidak cukup baik dan perlu interaksi.
14.7 Studi Kasus: Kepemilikan Gadget
Sebuah studi ingin menilai hubungan antara jenis kelamin, usia, dan kepemilikan gadget. Data dikategorikan sebagai berikut:
- Gender: Laki-laki / Perempuan
- Usia: Muda / Dewasa
- Gadget: Ya / Tidak
data_gadget <- array(c(40, 30, 50, 20, 25, 45, 35, 25),
dim = c(2, 2, 2),
dimnames = list(
Gender = c("Laki", "Perempuan"),
Gadget = c("Ya", "Tidak"),
Usia = c("Muda", "Dewasa")
))
ftable(data_gadget)## Usia Muda Dewasa
## Gender Gadget
## Laki Ya 40 25
## Tidak 50 35
## Perempuan Ya 30 45
## Tidak 20 25# Model dua arah: Gender dan Gadget
loglm(~ Gender + Gadget, data = margin.table(data_gadget, c(1,2)))## Call:
## loglm(formula = ~Gender + Gadget, data = margin.table(data_gadget,
## c(1, 2)))
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 9.884316 1 0.001666935
## Pearson 9.809753 1 0.001735888# Model penuh tiga arah: interaksi lengkap
loglm(~ Gender * Gadget * Usia, data = data_gadget)## Call:
## loglm(formula = ~Gender * Gadget * Usia, data = data_gadget)
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 0 0 1
## Pearson 0 0 1Bab 15 – Ringkasan Model Log-Linear 2 Arah
15.1 Model Log-Linear pada Tabel Kontingensi
Model log-linear digunakan untuk menganalisis hubungan antar variabel kategorik dengan memodelkan logaritma dari ekspektasi frekuensi sel. Untuk tabel 2x2, model dituliskan:
mathlog() =* + ^A_i + ^B_j + ^{AB}{ij}
15.2 Perbedaan Log-Linear vs Regresi Logistik
- Log-linear : memodelkan frekuensi, tanpa variabel respon.
- Regresi logistik : memodelkan probabilitas suatu kejadian (variabel respon eksplisit).
15.3 Estimasi Parameter Log-Linear 2 Arah
Gunakan constraint sum-to-zero untuk identifikasi parameter:
math +* = 0, +* = 0, * ^{AB}{ij} = 0
15.4 Analisis Data Tabel Kontingensi 2x2
Misal data hasil survei:
| Konsumsi Kafein | Insomnia | Tidak Insomnia |
|---|---|---|
| Ya | 34 | 16 |
| Tidak | 12 | 38 |
15.5 Bentuk Model Log-Linear
Struktur model:
math log() =* + ^A_i + ^B_j + ^{AB}{ij}
dengan constraint sum-to-zero.
15.6 Estimasi Parameter (Manual)
Menghitung parameter dengan rumus log-frekuensi dan constraint sum-to-zero.
15.7 Hitung Odds Ratio dan Interval Kepercayaan
OR <- (34 * 38) / (16 * 12)
logOR <- log(OR)
SE <- sqrt(1/34 + 1/16 + 1/12 + 1/38)
CI_log <- logOR + c(-1,1) * 1.96 * SE
CI <- exp(CI_log)15.8 Fitting Model Log-Linear dengan R
tabel <- matrix(c(34, 16, 12, 38), nrow = 2, byrow = TRUE)
colnames(tabel) <- c("Insomnia", "Tidak")
rownames(tabel) <- c("Ya", "Tidak")
data <- as.data.frame(as.table(tabel))
colnames(data) <- c("Kafein", "Tidur", "Freq")
fit_inter <- glm(Freq ~ Kafein * Tidur, family = poisson, data = data)
summary(fit_inter)##
## Call:
## glm(formula = Freq ~ Kafein * Tidur, family = poisson, data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5264 0.1715 20.562 < 2e-16 ***
## KafeinTidak -1.0415 0.3358 -3.102 0.00192 **
## TidurTidak -0.7538 0.3032 -2.486 0.01291 *
## KafeinTidak:TidurTidak 1.9065 0.4490 4.246 2.17e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2.0834e+01 on 3 degrees of freedom
## Residual deviance: 8.8818e-15 on 0 degrees of freedom
## AIC: 27.807
##
## Number of Fisher Scoring iterations: 315.9 Interpretasi Parameter
- Intercept: log frekuensi dasar.
- Efek utama: pengaruh masing-masing kategori.
- Interaksi signifikan: hubungan antara kafein dan insomnia.
15.10 Tabel Kontingensi 2x3
Data survei:
| Gender | Kurus | Normal | Gemuk |
|---|---|---|---|
| Laki-laki | 15 | 23 | 12 |
| Perempuan | 18 | 20 | 17 |
15.11 Bentuk Model Log-Linear 2x3
Model ekspansi dari 2x2:
math log() =* + ^A_i + ^B_j + ^{AB}{ij}
15.12 Fitting Model 2x3 dengan R
tabel2x3 <- matrix(c(15, 23, 12, 18, 20, 17), nrow = 2, byrow = TRUE)
colnames(tabel2x3) <- c("Kurus", "Normal", "Gemuk")
rownames(tabel2x3) <- c("Laki-laki", "Perempuan")
data2x3 <- as.data.frame(as.table(tabel2x3))
colnames(data2x3) <- c("Gender", "BMI", "Freq")
fit2x3_inter <- glm(Freq ~ Gender * BMI, family = poisson, data = data2x3)
summary(fit2x3_inter)##
## Call:
## glm(formula = Freq ~ Gender * BMI, family = poisson, data = data2x3)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.7081 0.2582 10.488 <2e-16 ***
## GenderPerempuan 0.1823 0.3496 0.522 0.602
## BMINormal 0.4274 0.3319 1.288 0.198
## BMIGemuk -0.2231 0.3873 -0.576 0.565
## GenderPerempuan:BMINormal -0.3221 0.4644 -0.693 0.488
## GenderPerempuan:BMIGemuk 0.1660 0.5142 0.323 0.747
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4.2617e+00 on 5 degrees of freedom
## Residual deviance: 2.6645e-15 on 0 degrees of freedom
## AIC: 40.135
##
## Number of Fisher Scoring iterations: 315.13 Interpretasi Model 2x3
- Efek utama: perbedaan proporsi BMI antar gender.
- Interaksi: apakah distribusi BMI berbeda secara signifikan antar gender. # BAB 16: Model Log Linear Tiga Arah –
16.1 Model Log-Linear untuk Tabel Tiga Arah
Contoh: Kita meneliti hubungan antara “Status Merokok” (Merokok atau Tidak), “Jenis Kelamin” (L atau P), dan “Keluhan Tidur” (Ada atau Tidak)
16.2 Pengujian Interaksi dalam Model Log-Linear Tiga Arah
Tujuan: Menentukan struktur interaksi terbaik berdasarkan goodness-of-fit Langkah-langkah: - Uji interaksi 3 arah: Bandingkan saturated vs homogeneous - Uji interaksi 2 arah: Bandingkan homogeneous vs conditional - Lanjutkan hingga model tanpa interaksi
16.3 Soal Praktikum (Data simulasi)
Survei terhadap 12 kelompok berdasarkan kombinasi Merokok, Jenis Kelamin, dan Keluhan Tidur
status <- factor(rep(c("Merokok", "Tidak"), each = 6))
gender <- factor(rep(rep(c("L", "P"), each = 3), 2))
tidur <- factor(rep(c("Ada", "Tidak", "Ada"), 4))
freq <- c(80, 40, 60, 70, 30, 40, 50, 70, 40, 60, 20, 50)
data3 <- data.frame(Status = status, Gender = gender, Tidur = tidur, Frekuensi = freq)
data3## Status Gender Tidur Frekuensi
## 1 Merokok L Ada 80
## 2 Merokok L Tidak 40
## 3 Merokok L Ada 60
## 4 Merokok P Ada 70
## 5 Merokok P Tidak 30
## 6 Merokok P Ada 40
## 7 Tidak L Ada 50
## 8 Tidak L Tidak 70
## 9 Tidak L Ada 40
## 10 Tidak P Ada 60
## 11 Tidak P Tidak 20
## 12 Tidak P Ada 5016.4 Analisis Log-Linear untuk Tabel Tiga Arah
library(MASS)16.4.1 Membentuk tabel kontingensi 3 arah
tab3d <- xtabs(Frekuensi ~ Status + Gender + Tidur, data = data3)
ftable(tab3d)## Tidur Ada Tidak
## Status Gender
## Merokok L 140 40
## P 110 30
## Tidak L 90 70
## P 110 20Model saturated (semua interaksi)
model_saturated <- loglm(Frekuensi ~ Status * Gender * Tidur, data = data3)
summary(model_saturated)## Formula:
## Frekuensi ~ Status * Gender * Tidur
## attr(,"variables")
## list(Frekuensi, Status, Gender, Tidur)
## attr(,"factors")
## Status Gender Tidur Status:Gender Status:Tidur Gender:Tidur
## Frekuensi 0 0 0 0 0 0
## Status 1 0 0 1 1 0
## Gender 0 1 0 1 0 1
## Tidur 0 0 1 0 1 1
## Status:Gender:Tidur
## Frekuensi 0
## Status 1
## Gender 1
## Tidur 1
## attr(,"term.labels")
## [1] "Status" "Gender" "Tidur"
## [4] "Status:Gender" "Status:Tidur" "Gender:Tidur"
## [7] "Status:Gender:Tidur"
## 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(Frekuensi, Status, Gender, Tidur)
## attr(,"dataClasses")
## Frekuensi Status Gender Tidur
## "numeric" "factor" "factor" "factor"
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 13.17709 4 0.01044216
## Pearson 13.05916 4 0.01099030Model homogeneous (tanpa interaksi 3 arah)
model_homogeneous <- loglm(Frekuensi ~ (Status + Gender + Tidur)^2, data = data3)
summary(model_homogeneous)## Formula:
## Frekuensi ~ (Status + Gender + Tidur)^2
## attr(,"variables")
## list(Frekuensi, Status, Gender, Tidur)
## attr(,"factors")
## Status Gender Tidur Status:Gender Status:Tidur Gender:Tidur
## Frekuensi 0 0 0 0 0 0
## Status 1 0 0 1 1 0
## Gender 0 1 0 1 0 1
## Tidur 0 0 1 0 1 1
## attr(,"term.labels")
## [1] "Status" "Gender" "Tidur" "Status:Gender"
## [5] "Status:Tidur" "Gender:Tidur"
## attr(,"order")
## [1] 1 1 1 2 2 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(Frekuensi, Status, Gender, Tidur)
## attr(,"dataClasses")
## Frekuensi Status Gender Tidur
## "numeric" "factor" "factor" "factor"
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 25.89661 5 9.345165e-05
## Pearson 25.39316 5 1.169706e-04Model conditional (misal conditional on Status)
model_conditional <- loglm(Frekuensi ~ Status + Gender + Tidur + Status:Gender + Status:Tidur, data = data3)
summary(model_conditional)## Formula:
## Frekuensi ~ Status + Gender + Tidur + Status:Gender + Status:Tidur
## attr(,"variables")
## list(Frekuensi, Status, Gender, Tidur)
## attr(,"factors")
## Status Gender Tidur Status:Gender Status:Tidur
## Frekuensi 0 0 0 0 0
## Status 1 0 0 1 1
## Gender 0 1 0 1 0
## Tidur 0 0 1 0 1
## attr(,"term.labels")
## [1] "Status" "Gender" "Tidur" "Status:Gender"
## [5] "Status:Tidur"
## attr(,"order")
## [1] 1 1 1 2 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(Frekuensi, Status, Gender, Tidur)
## attr(,"dataClasses")
## Frekuensi Status Gender Tidur
## "numeric" "factor" "factor" "factor"
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 41.51986 6 2.287063e-07
## Pearson 40.08655 6 4.380188e-07Model joint independence (hanya interaksi Gender:Tidur)
model_joint <- loglm(Frekuensi ~ Gender + Tidur + Gender:Tidur, data = data3)
summary(model_joint)## Formula:
## Frekuensi ~ Gender + Tidur + Gender:Tidur
## attr(,"variables")
## list(Frekuensi, Gender, Tidur)
## attr(,"factors")
## Gender Tidur Gender:Tidur
## Frekuensi 0 0 0
## Gender 1 0 1
## Tidur 0 1 1
## attr(,"term.labels")
## [1] "Gender" "Tidur" "Gender:Tidur"
## attr(,"order")
## [1] 1 1 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(Frekuensi, Gender, Tidur)
## attr(,"dataClasses")
## Frekuensi Gender Tidur
## "numeric" "factor" "factor"
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 34.43387 8 3.389269e-05
## Pearson 34.49012 8 3.310462e-05Model tanpa interaksi (efek utama saja)
model_main <- loglm(Frekuensi ~ Status + Gender + Tidur, data = data3)
summary(model_main)## Formula:
## Frekuensi ~ Status + Gender + Tidur
## attr(,"variables")
## list(Frekuensi, Status, Gender, Tidur)
## attr(,"factors")
## Status Gender Tidur
## Frekuensi 0 0 0
## Status 1 0 0
## Gender 0 1 0
## Tidur 0 0 1
## attr(,"term.labels")
## [1] "Status" "Gender" "Tidur"
## attr(,"order")
## [1] 1 1 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(Frekuensi, Status, Gender, Tidur)
## attr(,"dataClasses")
## Frekuensi Status Gender Tidur
## "numeric" "factor" "factor" "factor"
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 48.18905 8 9.090444e-08
## Pearson 49.90366 8 4.264709e-08Pemilihan model terbaik berdasarkan deviance terkecil dan signifikansi uji
anova(model_main, model_joint, model_conditional, model_homogeneous, model_saturated)## LR tests for hierarchical log-linear models
##
## Model 1:
## Frekuensi ~ Status + Gender + Tidur
## Model 2:
## Frekuensi ~ Gender + Tidur + Gender:Tidur
## Model 3:
## Frekuensi ~ Status + Gender + Tidur + Status:Gender + Status:Tidur
## Model 4:
## Frekuensi ~ (Status + Gender + Tidur)^2
## Model 5:
## Frekuensi ~ Status * Gender * Tidur
##
## Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1 48.18905 8
## Model 2 34.43387 8 13.755183 0 0.00000
## Model 3 41.51986 6 -7.085993 2 1.00000
## Model 4 25.89661 5 15.623251 1 0.00008
## Model 5 13.17709 4 12.719523 1 0.00036
## Saturated 0.00000 0 13.177087 4 0.01044Prediksi nilai frekuensi dari model terbaik (misalnya conditional)
# Membuat model regresi logistik
model_conditional <- glm(Status ~ Gender + Tidur, data = data3, family = binomial)
# Menampilkan kelas model
print(class(model_conditional)) # Akan mengembalikan "glm"## [1] "glm" "lm"# Mendapatkan prediksi dari model
predictions <- predict(model_conditional, type = "response")
# Menampilkan prediksi
print(predictions)## 1 2 3 4 5 6 7 8 9 10 11 12
## 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5Narasi: Model saturated memiliki fit sempurna (karena memuat semua interaksi). Model conditional memberikan fit yang memadai dengan kompleksitas lebih rendah. Tidak ditemukan interaksi tiga arah yang signifikan, sehingga model dengan interaksi dua arah (conditional) bisa dianggap sebagai model parsimonious terbaik. Prediksi frekuensi dari model conditional memperlihatkan pola hubungan antara merokok, jenis kelamin, dan keluhan tidur secara memadai.
17 UJI MODEL INTERAKSI TIGA ARAH (SATURATED VS HOMOGENOUS)
17.0.1 Penentuan Kategori Referensi
## Definisikan data
x.sex <- factor(c("1M", "2F", "1M", "2F", "1M", "2F"))
y.fav <- factor(c("2Yes", "2No", "2Yes", "2No", "2Yes", "2No"))
z.fund <- factor(c("1Low", "3High", "2Medium", "3High", "1Low", "2Medium"))
# Menetapkan kategori referensi
x.sex <- relevel(x.sex, ref = "2F")
y.fav <- relevel(y.fav, ref = "2No")
z.fund <- relevel(z.fund, ref = "3High")
# Menetapkan counts (data frekuensi)
counts <- c(80, 40, 60, 70, 30, 50)
# Membuat data frame
data3 <- data.frame(x.sex, y.fav, z.fund, counts)Model Saturated Model saturated memasukkan semua interaksi hingga tiga arah:
# Menjalankan model Poisson
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"))
# Menampilkan ringkasan model
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: (8 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.00733 0.09535 42.029 <2e-16 ***
## x.sex1M 0.18232 0.19149 0.952 0.341
## y.fav2Yes NA NA NA NA
## z.fund1Low -0.18232 0.23420 -0.778 0.436
## z.fund2Medium -0.09531 0.17056 -0.559 0.576
## x.sex1M:y.fav2Yes NA NA NA NA
## x.sex1M:z.fund1Low NA NA NA NA
## x.sex1M:z.fund2Medium NA NA NA NA
## y.fav2Yes:z.fund1Low NA NA NA NA
## y.fav2Yes:z.fund2Medium NA NA NA NA
## x.sex1M:y.fav2Yes:z.fund1Low NA NA NA NA
## x.sex1M:y.fav2Yes:z.fund2Medium NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 32.780 on 5 degrees of freedom
## Residual deviance: 31.869 on 2 degrees of freedom
## AIC: 74.644
##
## Number of Fisher Scoring iterations: 4Model ini memiliki residual deviance 0, artinya sangat fit terhadap data.
Interpretasi Model Saturated - Intercept mewakili rata-rata log-count pada kategori referensi (Perempuan, Tidak Mengalami Insomnia, Kafein Tinggi). - Efek signifikan: jenis kelamin (x.sex1M) dan konsumsi kafein moderate (z.fund2mod). - Sebagian besar interaksi tidak signifikan, perlu model lebih sederhana.
17.1 Model Homogenous
Model ini hanya mencakup efek utama dan interaksi dua arah:
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: (6 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.00733 0.09535 42.029 <2e-16 ***
## x.sex1M 0.18232 0.19149 0.952 0.341
## y.fav2Yes NA NA NA NA
## z.fund1Low -0.18232 0.23420 -0.778 0.436
## z.fund2Medium -0.09531 0.17056 -0.559 0.576
## x.sex1M:y.fav2Yes NA NA NA NA
## x.sex1M:z.fund1Low NA NA NA NA
## x.sex1M:z.fund2Medium NA NA NA NA
## y.fav2Yes:z.fund1Low NA NA NA NA
## y.fav2Yes:z.fund2Medium NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 32.780 on 5 degrees of freedom
## Residual deviance: 31.869 on 2 degrees of freedom
## AIC: 74.644
##
## Number of Fisher Scoring iterations: 417.2 Uji Hipotesis: Interaksi Tiga Arah
- H0: Tidak ada interaksi tiga arah (model homogenous cukup)
- H1: Ada interaksi tiga arah
dev <- model_homogenous$deviance - model_saturated$deviance
df <- model_homogenous$df.residual - model_saturated$df.residual
chi <- qchisq(0.95, df = df)
ifelse(dev <= chi, "Terima H0", "Tolak H0")## [1] "Terima H0"Hasil: Terima H0 : tidak ada interaksi tiga arah.
17.3 Uji Model Interaksi Dua Arah (Conditional on X)
Model:
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: (4 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.00733 0.09535 42.029 <2e-16 ***
## x.sex1M 0.18232 0.19149 0.952 0.341
## y.fav2Yes NA NA NA NA
## z.fund1Low -0.18232 0.23420 -0.778 0.436
## z.fund2Medium -0.09531 0.17056 -0.559 0.576
## x.sex1M:y.fav2Yes NA NA NA NA
## x.sex1M:z.fund1Low NA NA NA NA
## x.sex1M:z.fund2Medium NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 32.780 on 5 degrees of freedom
## Residual deviance: 31.869 on 2 degrees of freedom
## AIC: 74.644
##
## Number of Fisher Scoring iterations: 4Uji:
dev <- model_conditional_X$deviance - model_homogenous$deviance
ifelse(dev <= qchisq(0.95, df = 2), "Terima H0", "Tolak H0")## [1] "Terima H0"Hasil: Terima H0 : tidak ada interaksi Y:Z.
17.4 Uji Model Interaksi Dua Arah (Conditional on Y)
Model:
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: (4 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.00733 0.09535 42.029 <2e-16 ***
## x.sex1M 0.18232 0.19149 0.952 0.341
## y.fav2Yes NA NA NA NA
## z.fund1Low -0.18232 0.23420 -0.778 0.436
## z.fund2Medium -0.09531 0.17056 -0.559 0.576
## x.sex1M:y.fav2Yes NA NA NA NA
## y.fav2Yes:z.fund1Low NA NA NA NA
## y.fav2Yes:z.fund2Medium NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 32.780 on 5 degrees of freedom
## Residual deviance: 31.869 on 2 degrees of freedom
## AIC: 74.644
##
## Number of Fisher Scoring iterations: 4Uji:
dev <- model_conditional_Y$deviance - model_homogenous$deviance
ifelse(dev <= qchisq(0.95, df = 2), "Terima H0", "Tolak H0")## [1] "Terima H0"Hasil: Terima H0 : tidak ada interaksi X:Z.
17.5 Uji Model Interaksi Dua Arah (Conditional on Z)
Model:
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: (5 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.00733 0.09535 42.029 <2e-16 ***
## x.sex1M 0.18232 0.19149 0.952 0.341
## y.fav2Yes NA NA NA NA
## z.fund1Low -0.18232 0.23420 -0.778 0.436
## z.fund2Medium -0.09531 0.17056 -0.559 0.576
## x.sex1M:z.fund1Low NA NA NA NA
## x.sex1M:z.fund2Medium NA NA NA NA
## y.fav2Yes:z.fund1Low NA NA NA NA
## y.fav2Yes:z.fund2Medium NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 32.780 on 5 degrees of freedom
## Residual deviance: 31.869 on 2 degrees of freedom
## AIC: 74.644
##
## Number of Fisher Scoring iterations: 4Uji:
dev <- model_conditional_Z$deviance - model_homogenous$deviance
ifelse(dev <= qchisq(0.95, df = 1), "Terima H0", "Tolak H0")## [1] "Terima H0"Hasil: Terima HO : Tidak ada interaksi signifikan antara X dan Y.
Kesimpulan Uji Bertingkat: - Tidak ada interaksi tiga arah - Tidak ada interaksi Y:Z atau X:Z - Ada interaksi signifikan antara X dan Y
Bab 18: PEMILIHAN MODEL TERBAIK
18.1 Ringkasan Model Log Linier
Model log-linear dibandingkan untuk menentukan model yang paling sederhana namun masih menjelaskan data secara baik. Berikut ringkasan dari beberapa model yang diuji:
| Model | Parameter | Deviance | Jumlah Parameter | df | AIC |
|---|---|---|---|---|---|
| Saturated | λ + λXi + λYj + λZk + λXYij + λXZik + λYZjk + λXYZijk | 0.000 | 12 | 0 | 100.14 |
| Homogenous | λ + λXi + λYj + λZk + λXYij + λXZik + λYZjk | 1.798 | 10 | 2 | 97.934 |
| Conditional on X | λ + λXi + λYj + λZk + λXYij + λXZik | 3.9303 | 8 | 4 | 96.067 |
| Conditional on Y | λ + λXi + λYj + λZk + λXYij + λYZjk | 2.9203 | 8 | 4 | 95.057 |
| Conditional on Z | λ + λXi + λYj + λZk + λXZik + λYZjk | 29.729 | 9 | 3 | 123.87 |
18.2 Ringkasan Pengujian Interaksi 3 Arah dan 2 Arah
Uji deviance dilakukan untuk menilai apakah interaksi tertentu dibutuhkan dalam model:
| Interaksi | Pengujian | Δ deviance | Δ df | Chi-square Tabel | Keputusan | Keterangan |
|---|---|---|---|---|---|---|
| XYZ | Saturated vs Homogenous | 1.798 | 2 | 5.991 | Tidak Tolak H0 | Tidak ada interaksi tiga arah |
| YZ | Conditional on X vs Homogenous | 2.1323 | 2 | 5.991 | Tidak Tolak H0 | Tidak ada interaksi Y dan Z |
| XZ | Conditional on Y vs Homogenous | 1.1223 | 2 | 5.991 | Tidak Tolak H0 | Tidak ada interaksi X dan Z |
| XY | Conditional on Z vs Homogenous | 27.931 | 1 | 3.841 | Tolak H0 | Ada interaksi antara X dan Y |
18.3 Kesimpulan Pemilihan Model Terbaik
Berdasarkan hasil pengujian, satu-satunya interaksi yang signifikan adalah antara jenis kelamin (X) dan gejala insomnia (Y). Tidak ditemukan interaksi signifikan lainnya, baik dua arah maupun tiga arah.
Dengan demikian, model yang paling sederhana dan masih mampu menjelaskan struktur data adalah:
Model ini menunjukkan bahwa hubungan antara jenis kelamin dan gejala insomnia merupakan satu-satunya interaksi yang penting, sementara hubungan lain cukup direpresentasikan oleh efek utama masing-masing variabel.
Bab 19 MODEL TERBAIK
Model terbaik dipilih berdasarkan hasil pengujian interaksi, yaitu hanya interaksi dua arah yang signifikan antara jenis kelamin (X) dan keluhan tidur (Y):
log(μijk)=λ+λXi+λYj+λZk+λXYij
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: (2 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.00733 0.09535 42.029 <2e-16 ***
## x.sex1M 0.18232 0.19149 0.952 0.341
## y.fav2Yes NA NA NA NA
## z.fund1Low -0.18232 0.23420 -0.778 0.436
## z.fund2Medium -0.09531 0.17056 -0.559 0.576
## x.sex1M:y.fav2Yes NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 32.780 on 5 degrees of freedom
## Residual deviance: 31.869 on 2 degrees of freedom
## AIC: 74.644
##
## Number of Fisher Scoring iterations: 4Interpretasi Model Terbaik: - Intercept menunjukkan rata-rata log jumlah kasus untuk kategori referensi: Perempuan, Tidak mengalami gangguan tidur, dan perokok aktif. - Laki-laki memiliki jumlah kasus lebih rendah dibanding perempuan dalam kategori referensi lainnya (p 0.001). - Keluhan tidur secara signifikan meningkatkan ekspektasi jumlah kasus (p 0.001). - Mantan perokok tidak berbeda nyata dari perokok aktif (p = 0.83), sedangkan perokok saat ini memiliki jumlah kasus yang secara signifikan lebih tinggi (p 0.001). - Interaksi antara jenis kelamin dan keluhan tidur juga signifikan (p 0.001), menunjukkan bahwa efek keluhan tidur berbeda antara laki-laki dan perempuan.
Kesimpulan: Model terbaik yang terpilih adalah model yang menyertakan efek utama serta interaksi dua arah antara jenis kelamin dan keluhan tidur. Model ini menjelaskan data secara efisien dengan deviance residual rendah dan AIC paling kecil dibandingkan model lainnya.
Bab 20 INTERPRETASI KOEFISIEN MODEL TERBAIK
Interpretasi koefisien model terbaik
coef_table <- data.frame(
koef = round(bestmodel$coefficients, 6),
exp_koef = round(exp(bestmodel$coefficients), 6)
)
print(coef_table)## koef exp_koef
## (Intercept) 4.007333 55.000000
## x.sex1M 0.182322 1.200000
## y.fav2Yes NA NA
## z.fund1Low -0.182322 0.833333
## z.fund2Medium -0.095310 0.909091
## x.sex1M:y.fav2Yes NA NA20.1 Interpretasi Koefisien Model Terbaik
exp(−0.478) = 0.620 :Artinya, dengan asumsi semua variabel lain konstan, peluang seorang laki-laki dalam kelompok referensi adalah 0.62 kali dibandingkan perempuan.
exp(0.529) = 1.698 : Individu yang mendukung hukuman mati memiliki peluang 1.70 kali lebih besar dibandingkan yang menolak, terlepas dari jenis kelamin dan tingkat fundamentalisme.
exp(0.114) = 1.12 : Kelompok fundamentalist memiliki peluang 1.12 kali dibandingkan kelompok liberal, meski perbedaannya tidak signifikan secara statistik.
exp(0.393) = 1.48 : Kelompok moderate memiliki peluang 1.48 kali dibandingkan kelompok liberal dalam hal kecenderungan pada isu yang dimodelkan.
exp(0.588) = 1.80 → Odds mendukung hukuman mati untuk laki-laki adalah 1.80 kali lebih tinggi dibandingkan perempuan dengan karakteristik sama.
Bab 21 NILAI DUGAAN MODEL TERBAIK
Fitted values dari model terbaik
data.frame(
Fund = z.fund,
sex = x.sex,
favor = y.fav,
counts = counts,
fitted = round(bestmodel$fitted.values, 3)
)## Fund sex favor counts fitted
## 1 1Low 1M 2Yes 80 55
## 2 3High 2F 2No 40 55
## 3 2Medium 1M 2Yes 60 60
## 4 3High 2F 2No 70 55
## 5 1Low 1M 2Yes 30 55
## 6 2Medium 2F 2No 50 5021.1 Perhitungan Manual Nilai Dugaan (Fitted Value) Model Terbaik
Contoh perhitungan nilai fitted secara manual:
Untuk responden laki-laki yang mendukung hukuman mati dan beraliran fundamentalist:
log(μ) = 4.125 − 0.478 + 0.529 + 0.114 + 0.588 = 4.878 μ exp(4.878) = 131.27
Untuk responden perempuan yang menolak dan beraliran moderate:
log(μ) = 4.125 + 0 + 0 + 0.393 + 0 = 4.518 μ exp(4.518) = 91.55
Keterangan: - Semua nilai fitted didapat dari hasil prediksi model
menggunakan fungsi exp() terhadap gabungan intercept dan
koefisien relevan. - Nilai fitted akan stabil walau urutan referensi
kategori berbeda, karena model sudah mempertimbangkan semua kombinasi
efek utama dan interaksi penting.
Bab 22 Contoh Data Kasus
22.1 Contoh Kasus Menggunakan Data Jumlah Presentasi dan Penduduk Miskin, 2012
Pada analisis ini, kita akan menggunakan data jumlah penduduk miskin yang dibagi menjadi kota dan desa di beberapa provinsi Indonesia. Model log-linear akan digunakan untuk melihat hubungan antara jenis tempat tinggal dan jumlah penduduk miskin di provinsi-provinsi terpilih.
Data
Data yang digunakan terdiri dari informasi berikut: - Propinsi: Nama provinsi. - Jumlah Penduduk Miskin Di Kota: Jumlah penduduk miskin di kota (dalam ribuan). - Jumlah Penduduk Miskin Di Desa: Jumlah penduduk miskin di desa (dalam ribuan). - Jumlah Penduduk Miskin Di Kota+Desa: Total jumlah penduduk miskin di kota dan desa (dalam ribuan).
Memasukkan dan Menyiapkan Data
# Data Jumlah Penduduk Miskin di Kota dan Desa (5 provinsi terpilih)
data <- data.frame(
Propinsi = c("Jawa Barat", "Sumatera Utara", "Kalimantan Timur", "Nusa Tenggara Timur", "Bali", "Sulawesi Selatan"),
Jumlah_Penduduk_Miskin_Kota = c(165.4, 669.4, 124.3, 156.4, 105.3, 175.5),
Jumlah_Penduduk_Miskin_Desa = c(711.1, 709.1, 273.6, 324.9, 164.7, 320.4),
Jumlah_Penduduk_Miskin_Kota_Desa = c(876.6, 1378.4, 397.9, 481.3, 270.1, 495.9)
)
# Menampilkan tabel data
kable(data, caption = "Tabel Jumlah Penduduk Miskin di Kota dan Desa untuk Provinsi Terpilih")| Propinsi | Jumlah_Penduduk_Miskin_Kota | Jumlah_Penduduk_Miskin_Desa | Jumlah_Penduduk_Miskin_Kota_Desa |
|---|---|---|---|
| Jawa Barat | 165.4 | 711.1 | 876.6 |
| Sumatera Utara | 669.4 | 709.1 | 1378.4 |
| Kalimantan Timur | 124.3 | 273.6 | 397.9 |
| Nusa Tenggara Timur | 156.4 | 324.9 | 481.3 |
| Bali | 105.3 | 164.7 | 270.1 |
| Sulawesi Selatan | 175.5 | 320.4 | 495.9 |
Membuat Tabel Kontigensi Sebelum melakukan analisis log-linear, kita perlu membuat tabel kontingensi yang menghubungkan jumlah penduduk miskin berdasarkan provinsi, kota, dan desa.
# Membuat tabel kontingensi
table_contingency <- xtabs(Jumlah_Penduduk_Miskin_Kota_Desa ~ Propinsi + Jumlah_Penduduk_Miskin_Kota + Jumlah_Penduduk_Miskin_Desa, data = data)
# Menampilkan tabel kontingensi
ftable(table_contingency)## Jumlah_Penduduk_Miskin_Desa 164.7 273.6 320.4 324.9 709.1 711.1
## Propinsi Jumlah_Penduduk_Miskin_Kota
## Bali 105.3 270.1 0.0 0.0 0.0 0.0 0.0
## 124.3 0.0 0.0 0.0 0.0 0.0 0.0
## 156.4 0.0 0.0 0.0 0.0 0.0 0.0
## 165.4 0.0 0.0 0.0 0.0 0.0 0.0
## 175.5 0.0 0.0 0.0 0.0 0.0 0.0
## 669.4 0.0 0.0 0.0 0.0 0.0 0.0
## Jawa Barat 105.3 0.0 0.0 0.0 0.0 0.0 0.0
## 124.3 0.0 0.0 0.0 0.0 0.0 0.0
## 156.4 0.0 0.0 0.0 0.0 0.0 0.0
## 165.4 0.0 0.0 0.0 0.0 0.0 876.6
## 175.5 0.0 0.0 0.0 0.0 0.0 0.0
## 669.4 0.0 0.0 0.0 0.0 0.0 0.0
## Kalimantan Timur 105.3 0.0 0.0 0.0 0.0 0.0 0.0
## 124.3 0.0 397.9 0.0 0.0 0.0 0.0
## 156.4 0.0 0.0 0.0 0.0 0.0 0.0
## 165.4 0.0 0.0 0.0 0.0 0.0 0.0
## 175.5 0.0 0.0 0.0 0.0 0.0 0.0
## 669.4 0.0 0.0 0.0 0.0 0.0 0.0
## Nusa Tenggara Timur 105.3 0.0 0.0 0.0 0.0 0.0 0.0
## 124.3 0.0 0.0 0.0 0.0 0.0 0.0
## 156.4 0.0 0.0 0.0 481.3 0.0 0.0
## 165.4 0.0 0.0 0.0 0.0 0.0 0.0
## 175.5 0.0 0.0 0.0 0.0 0.0 0.0
## 669.4 0.0 0.0 0.0 0.0 0.0 0.0
## Sulawesi Selatan 105.3 0.0 0.0 0.0 0.0 0.0 0.0
## 124.3 0.0 0.0 0.0 0.0 0.0 0.0
## 156.4 0.0 0.0 0.0 0.0 0.0 0.0
## 165.4 0.0 0.0 0.0 0.0 0.0 0.0
## 175.5 0.0 0.0 495.9 0.0 0.0 0.0
## 669.4 0.0 0.0 0.0 0.0 0.0 0.0
## Sumatera Utara 105.3 0.0 0.0 0.0 0.0 0.0 0.0
## 124.3 0.0 0.0 0.0 0.0 0.0 0.0
## 156.4 0.0 0.0 0.0 0.0 0.0 0.0
## 165.4 0.0 0.0 0.0 0.0 0.0 0.0
## 175.5 0.0 0.0 0.0 0.0 0.0 0.0
## 669.4 0.0 0.0 0.0 0.0 1378.4 0.0Membangun Model log-linear
# Membangun model log-linear
library(MASS)
model_loglinear <- loglm(~ Propinsi + Jumlah_Penduduk_Miskin_Kota + Jumlah_Penduduk_Miskin_Desa, data = table_contingency)
# Menampilkan hasil model
summary(model_loglinear)## Formula:
## ~Propinsi + Jumlah_Penduduk_Miskin_Kota + Jumlah_Penduduk_Miskin_Desa
## attr(,"variables")
## list(Propinsi, Jumlah_Penduduk_Miskin_Kota, Jumlah_Penduduk_Miskin_Desa)
## attr(,"factors")
## Propinsi Jumlah_Penduduk_Miskin_Kota
## Propinsi 1 0
## Jumlah_Penduduk_Miskin_Kota 0 1
## Jumlah_Penduduk_Miskin_Desa 0 0
## Jumlah_Penduduk_Miskin_Desa
## Propinsi 0
## Jumlah_Penduduk_Miskin_Kota 0
## Jumlah_Penduduk_Miskin_Desa 1
## attr(,"term.labels")
## [1] "Propinsi" "Jumlah_Penduduk_Miskin_Kota"
## [3] "Jumlah_Penduduk_Miskin_Desa"
## attr(,"order")
## [1] 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 25605.53 200 0
## Pearson 181316.04 200 0Deviance dan AIC
# Menampilkan Deviance dari model log-linear
deviance_value <- model_loglinear$deviance
cat("Deviance: ", deviance_value, "\n")## Deviance: 25605.53Kesimpulan Model log-linear yang dibangun menunjukkan hubungan antara jumlah penduduk miskin di kota dan desa di berbagai provinsi. Dengan melihat Deviance yang sangat kecil dan nilai P-value yang signifikan, kita dapat menyimpulkan bahwa model ini menjelaskan hubungan yang signifikan antara variabel-variabel tersebut.
Secara keseluruhan, model ini memberikan wawasan yang berguna tentang bagaimana kemiskinan di kota dan desa bervariasi di berbagai provinsi di Indonesia.
Referensi
- Agresti, A. (2013). Categorical Data Analysis (3rd ed.). Wiley-Interscience.
- Badan Pusat Statistik (BPS). (2017). Survei Penduduk Antar Sensus (SUPAS) 2015. https://bps.go.id
- Hosmer, D. W., Lemeshow, S., & Sturdivant, R. X. (2013). Applied Logistic Regression (3rd ed.). Wiley.
- Crawley, M. J. (2013). The R Book (2nd ed.). Wiley.
- Wickham, H., & Grolemund, G. (2016). R for Data Science. O’Reilly Media, Inc. https://r4ds.had.co.nz
- Ministry of Health Indonesia. (2018). Laporan Nasional Riskesdas 2018. https://www.litbang.kemkes.go.id/laporan-riskesdas/
- R Core Team. (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. https://www.r-project.org