EBOOK ADK
Mona Yola Lumban Raja (140610230081)
2025-04-14
1 Pendahuluan
Dalam dunia statistik, data diklasifikasikan ke dalam beberapa jenis, salah satunya adalah data kategori. Data kategori (categorical data) adalah jenis data yang merepresentasikan kelompok atau kategori tanpa nilai numerik yang memiliki makna matematis atau urutan tetap. Berbeda dengan data numerik, data ini tidak digunakan untuk menghitung rata-rata atau selisih, namun lebih digunakan untuk mengelompokkan informasi. Contoh umum data kategori meliputi jenis kelamin (laki-laki/perempuan), tingkat pendidikan (SD, SMP, SMA, S1), status pekerjaan, warna favorit, dan sebagainya.
Analisis terhadap data kategori sangat penting, terutama dalam ilmu sosial, kesehatan masyarakat, ekonomi, pendidikan, dan bidang lainnya. Dengan menganalisis data jenis ini, kita dapat memahami distribusi kategori dalam populasi, mengidentifikasi hubungan antara dua variabel kategori, serta membangun model prediktif yang berguna dalam pengambilan keputusan.
1.1 Tujuan Analisis Data Kategori
Analisis data kategori memiliki beberapa tujuan utama yang saling berkaitan, yaitu:
Menemukan Pola dan Tren Data kategori dapat membantu dalam mengenali kecenderungan tertentu dalam kelompok masyarakat atau sampel yang diteliti. Misalnya, kita bisa melihat bahwa sebagian besar pengguna media sosial berasal dari kelompok usia remaja dan dewasa muda.
Menganalisis Relasi antar Variabel Dengan menggunakan teknik seperti uji chi-square, kita bisa mengevaluasi apakah terdapat hubungan yang signifikan antara dua variabel kategori. Contohnya adalah hubungan antara kebiasaan merokok (ya/tidak) dengan kejadian penyakit jantung (ada/tidak).
Mendukung Pengambilan Keputusan Informasi yang diperoleh dari data kategori dapat digunakan sebagai dasar dalam membuat kebijakan atau strategi. Misalnya, pemerintah dapat menggunakan hasil survei tentang status vaksinasi masyarakat untuk menyusun strategi kampanye kesehatan.
Membangun Model Prediksi Dalam beberapa kasus, data kategori juga dapat digunakan dalam model prediksi seperti regresi logistik, khususnya saat kita ingin memprediksi kejadian biner seperti kemungkinan seseorang membeli produk atau tidak, berdasarkan kategori tertentu seperti jenis kelamin, pekerjaan, atau pendidikan.
1.2 Ruang Lingkup
Data kategori terbagi menjadi beberapa jenis utama berdasarkan karakteristiknya:
Data Nominal Merupakan data kategori yang tidak memiliki urutan atau tingkatan. Contohnya adalah jenis kelamin (laki-laki/perempuan), agama, dan warna mata. Tidak ada kategori yang “lebih tinggi” atau “lebih rendah”.
Data Ordinal Data ini memiliki tingkatan atau urutan, meskipun jaraknya tidak dapat diukur secara pasti. Contohnya adalah tingkat kepuasan (sangat puas, puas, netral, tidak puas, sangat tidak puas), atau tingkat pendidikan (SD, SMP, SMA, S1, S2, S3).
Data Biner vs Multikategori Biner (dichotomous): Hanya memiliki dua kategori, seperti ya/tidak, laki-laki/perempuan, hadir/tidak hadir.
Multikategori: Memiliki lebih dari dua kategori, contohnya adalah jenis pekerjaan (pegawai negeri, wiraswasta, pelajar, dll) atau warna favorit.
1.3 Manfaat di Berbagai Bidang
Penggunaan data kategori sangat luas dan bermanfaat di berbagai sektor, di antaranya:
Bidang Kesehatan: Untuk mengidentifikasi kelompok populasi dengan risiko tinggi terhadap penyakit tertentu berdasarkan faktor gaya hidup, kebiasaan makan, atau tingkat aktivitas fisik.
Bidang Pendidikan: Untuk mengevaluasi efektivitas metode pembelajaran berdasarkan tingkat pendidikan, persepsi siswa, dan hasil evaluasi belajar.
Bidang Sosial: Dalam survei opini publik, data kategori dapat menunjukkan bagaimana persepsi masyarakat terhadap isu-isu sosial berdasarkan kelompok usia, pendidikan, atau wilayah tempat tinggal.
Bidang Bisnis dan Pemasaran: Berguna dalam melakukan segmentasi pasar berdasarkan preferensi pelanggan, jenis kelamin, kelompok usia, dan pola konsumsi, sehingga dapat digunakan untuk strategi pemasaran yang lebih tepat sasaran.
Bidang Kriminalitas: Untuk menganalisis pola kejahatan berdasarkan jenis kelamin pelaku, usia, lokasi kejadian, atau jenis kejahatan yang dilakukan, sehingga dapat membantu perencanaan kebijakan keamanan.
2 Metode dalam Analisis Data Kategori
Analisis data kategori menyediakan berbagai pendekatan statistik untuk memahami pola, relasi, dan prediksi dalam data yang bersifat nominal atau ordinal. Beberapa metode yang paling umum dan penting dalam analisis ini antara lain:
2.1 Tabel Kontingensi dan Uji Chi-Square
Tabel kontingensi digunakan untuk menyajikan frekuensi dua variabel kategori secara bersilang. Melalui tabel ini, kita dapat mengamati bagaimana distribusi kategori satu variabel bervariasi terhadap kategori lainnya. Untuk menguji apakah ada hubungan yang signifikan antara dua variabel kategori tersebut, digunakan uji chi-square (χ² test). Uji ini membandingkan frekuensi yang diamati dengan frekuensi yang diharapkan jika tidak ada hubungan antar variabel.
2.2 Regresi Logistik
Regresi logistik merupakan teknik yang digunakan ketika variabel dependen (Y) bersifat biner atau kategorik. Model ini memprediksi probabilitas suatu kejadian berdasarkan variabel bebas, yang dapat berupa kategori maupun numerik. Contohnya, memprediksi kemungkinan seseorang akan membeli suatu produk (ya/tidak) berdasarkan jenis kelamin, status pekerjaan, dan frekuensi belanja.
2.3 Correspondence Analysis (CA)
Correspondence Analysis (CA) adalah metode eksploratori berbasis grafik yang digunakan untuk menganalisis hubungan antara dua atau lebih variabel kategori. Teknik ini menghasilkan visualisasi dalam bentuk peta dua dimensi, sehingga memudahkan interpretasi keterkaitan antar kategori.
2.4 Decision Tree dan Random Forest
Metode pohon keputusan (decision tree) dan random forest termasuk ke dalam pendekatan machine learning berbasis kategori.
Decision Tree menyusun aturan keputusan secara hierarkis berdasarkan variabel input kategorikal.
Random Forest adalah pengembangan dari decision tree yang menggunakan banyak pohon untuk menghasilkan klasifikasi atau prediksi yang lebih akurat dan stabil.
Metode ini sangat cocok untuk tujuan klasifikasi dan prediksi dalam data yang kompleks, termasuk untuk pengambilan keputusan bisnis, diagnosis medis, dan deteksi penipuan.
2.5 Pentingnya Penguasaan Metode Analisis Kategori
Analisis data kategori merupakan salah satu fondasi utama dalam statistik terapan. Relevansinya sangat tinggi terutama ketika kita bekerja dengan data yang bersifat diskrit, seperti dalam survei, kuisioner, atau data administratif. Penerapannya meluas di berbagai bidang kehidupan nyata, termasuk:
Riset kesehatan masyarakat: menilai keterkaitan antara gaya hidup dan penyakit.
Survei sosial-politik: memahami opini masyarakat terhadap isu tertentu.
Pemasaran dan bisnis: mengelompokkan konsumen berdasarkan preferensi.
Pendidikan: mengevaluasi persepsi atau kepuasan siswa.
Dengan dukungan perangkat lunak statistik seperti R, SPSS, Stata, dan Python, analisis data kategori dapat dilakukan secara efisien dan presisi. Di era digital saat ini, metode modern seperti Generalized Linear Models (GLM), model log-linear, serta pendekatan Bayesian pun mulai banyak digunakan untuk memperluas kedalaman analisis.
Kesimpulan Pemahaman terhadap metode-metode ini memungkinkan seorang analis tidak hanya untuk menggambarkan distribusi kategori dalam data, tetapi juga menguji hipotesis, memprediksi kecenderungan, dan menyusun strategi berbasis bukti. Oleh karena itu, penguasaan analisis data kategori menjadi kompetensi inti bagi peneliti, akademisi, pengambil kebijakan, maupun profesional dari berbagai sektor.
3 Distribusi Probabilitas dalam Data Kategori
Distribusi probabilitas merupakan dasar utama dalam statistik untuk memahami dan memodelkan ketidakpastian dari suatu kejadian. Dalam konteks data kategori, distribusi probabilitas berperan penting untuk menjelaskan bagaimana kemungkinan terjadinya suatu hasil dalam sejumlah kategori yang tersedia.
Distribusi ini sangat relevan saat berhadapan dengan data diskrit atau data yang terbagi ke dalam beberapa kelompok. Baik dalam survei sosial, eksperimen klinis, maupun klasifikasi produk, distribusi probabilitas kategori membantu para analis dalam memahami pola frekuensi, melakukan simulasi data, serta membangun inferensi statistik.
Beberapa distribusi penting yang digunakan untuk memodelkan data kategori antara lain:
3.1 Distribusi Bernoulli
Distribusi Bernoulli adalah distribusi probabilitas paling dasar yang hanya memiliki dua kemungkinan hasil: sukses (1) atau gagal (0). Distribusi ini digunakan untuk percobaan tunggal dengan dua hasil kemungkinan, seperti ya/tidak, hidup/mati, atau membeli/tidak membeli.
Fungsi probabilitasnya:
\[ 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)
library(kableExtra)
library(dplyr)##
## 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 perluasan dari Bernoulli untuk sejumlah percobaan 𝑛yang identik dan independen. Distribusi ini digunakan untuk menghitung jumlah sukses dari sejumlah percobaan dengan peluang yang konstan.
\[ 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 ini sangat berguna dalam pengujian hipotesis proporsi, analisis eksperimen klinis, serta studi epidemiologi.
# 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 sebuah percobaan memiliki lebih dari dua hasil kategori, kita menggunakan distribusi multinomial. Distribusi ini menggambarkan jumlah kejadian pada masing-masing kategori dari total percobaan.
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 sering digunakan dalam analisis preferensi, hasil klasifikasi, dan penelitian pasar.
# 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 suatu interval waktu atau ruang tertentu, khususnya ketika kejadian tersebut jarang terjadi (rare events). Distribusi ini hanya memiliki satu parameter, yaitu rata-rata kejadian 𝜆.
\[ P(Y = y) = \frac{e^{-\lambda} \lambda^y}{y!} \quad \text{untuk } y = 0, 1, 2, ... \]
Distribusi ini sering diterapkan dalam epidemiologi, manajemen lalu lintas, analisis antrean, dan logistik, seperti memodelkan jumlah pelanggan yang datang ke toko dalam satu jam.
# 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 analisis statistik, kualitas dan validitas hasil sangat bergantung pada bagaimana data dikumpulkan. Oleh karena itu, desain pengambilan sampel merupakan aspek krusial, terlebih saat bekerja dengan data kategori. Pemilihan teknik sampling yang tepat dapat memengaruhi keakuratan generalisasi hasil terhadap populasi yang lebih luas.
Dalam konteks analisis data kategori, dua pendekatan utama dalam desain sampling adalah sampling prospektif dan sampling retrospektif. Perbedaan utama dari kedua pendekatan ini terletak pada arah waktu pengumpulan data dan tujuan analisis.
4.1 Prospective Sampling
Sampling prospektif dilakukan dengan memilih responden terlebih dahulu, kemudian melakukan observasi ke depan dalam kurun waktu tertentu untuk melihat apakah suatu kejadian (outcome) akan muncul. Karena pengamatan dilakukan setelah paparan ditentukan, pendekatan ini memungkinkan kontrol yang lebih baik terhadap variabel-variabel yang memengaruhi hasil.
Pendekatan ini umum digunakan dalam penelitian eksperimental dan studi kohort jangka panjang.
4.1.1 Eksperimen
Dalam studi eksperimental, partisipan secara acak dibagi ke dalam kelompok perlakuan dan kontrol. Tujuannya adalah untuk mengevaluasi efek suatu intervensi secara kausal.
Beberapa teknik sampling yang umum digunakan:
Simple Random Sampling (SRS): Setiap individu dalam populasi memiliki peluang yang sama untuk terpilih sebagai sampel.
Stratified Random Sampling: Populasi dibagi ke dalam strata berdasarkan karakteristik tertentu (misalnya jenis kelamin, usia), kemudian sampel diambil secara acak dari tiap strata.
Cluster Sampling: Populasi dibagi ke dalam kelompok-kelompok (cluster), dan beberapa cluster dipilih secara acak untuk dijadikan sampel.
4.1.2 Studi Kohort
Studi kohort adalah studi observasional di mana sekelompok individu dengan karakteristik tertentu diikuti dari waktu ke waktu untuk melihat kemunculan outcome. Studi ini berguna untuk mengevaluasi hubungan antara paparan dan kejadian yang jarang namun signifikan.
Jenis teknik sampling yang umum digunakan:
Census Sampling: Seluruh anggota populasi target dijadikan sebagai sampel studi (misalnya semua pasien di rumah sakit tertentu selama 1 tahun).
Systematic Sampling: Sampel diambil berdasarkan interval tertentu dari daftar populasi, seperti memilih setiap responden ke-10 dari daftar peserta.
Matched Sampling: Subjek dari satu kelompok dipasangkan dengan subjek dari kelompok lain yang memiliki karakteristik serupa untuk mengontrol efek perancu.
4.2 Retrospective Sampling
Berbeda dengan pendekatan prospektif, sampling retrospektif dilakukan setelah kejadian atau outcome telah terjadi. Responden dikelompokkan berdasarkan status kejadian (misalnya: sakit/tidak sakit), lalu informasi tentang faktor penyebab atau paparan masa lalu dikumpulkan.
Pendekatan ini umum dalam studi kasus-kontrol dan cocok untuk menganalisis kejadian yang jarang terjadi atau membutuhkan waktu lama untuk berkembang.
4.2.1 Studi Kasus-Kontrol
Dalam studi kasus-kontrol, kelompok “kasus” (yang mengalami outcome) dibandingkan dengan kelompok “kontrol” (yang tidak mengalaminya). Studi ini bertujuan untuk mengidentifikasi faktor-faktor risiko atau penyebab yang mungkin berkaitan dengan kejadian tersebut.
Teknik sampling yang umum digunakan:
Purposive Sampling: Subjek dipilih secara sengaja berdasarkan kriteria yang relevan dengan fokus penelitian.
Snowball Sampling: Responden awal merekomendasikan subjek lain yang memenuhi kriteria, berguna saat populasi sulit dijangkau (misalnya pengguna narkoba).
Incidence Density Sampling: Kasus dan kontrol dipilih berdasarkan waktu kemunculan kejadian, biasanya untuk menjaga proporsionalitas dalam pengamatan.
Stratified dan Cluster Sampling: Dapat diterapkan untuk menjaga keberagaman dan representativitas sampel pada populasi yang luas atau tersebar secara geografis.
5 Tabel Kontingensi 2 × 2
Tabel kontingensi 2 × 2 merupakan alat analisis penting dalam statistik kategori, digunakan untuk menyajikan distribusi frekuensi dari dua variabel kategori secara silang. Analisis ini sangat umum dalam studi epidemiologi, penelitian sosial, serta kebijakan kesehatan masyarakat.
Dalam bab ini, kita akan menganalisis hubungan antara tingkat pendidikan ibu (Tinggi vs Rendah) dan status imunisasi dasar lengkap anak usia 12–23 bulan (Lengkap vs Tidak Lengkap), berdasarkan data dari Survei Demografi dan Kesehatan Indonesia (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 1500Dari tabel tersebut, kita peroleh:
Dari ibu berpendidikan tinggi, 1200 anak mendapat imunisasi lengkap dan 800 tidak.
Dari ibu berpendidikan rendah, hanya 1000 anak mendapat imunisasi lengkap dan 1500 tidak.
Total sampel: 4.500 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 dua kejadian terjadi secara bersamaan, seperti: anak dari ibu berpendidikan tinggi yang mendapat 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 menunjukkan kemungkinan suatu kejadian terjadi jika diketahui kejadian lainnya telah 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 mengevaluasi kekuatan hubungan antara tingkat pendidikan ibu dan status imunisasi anak, digunakan tiga ukuran asosiasi statistik:
5.3.1 Risk Difference (RD)
Mengukur selisih probabilitas kejadian antara dua kelompok:
\[ 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 risiko 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)
Membandingkan odds (peluang relatif) antar 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
Hasil analisis menunjukkan: Berdasarkan hasil perhitungan distribusi peluang dan ukuran asosiasi dari tabel kontingensi 2 × 2, diperoleh temuan sebagai berikut:
- Distribusi pendidikan ibu dalam sampel:
Ibu berpendidikan tinggi: 44,4%
Ibu berpendidikan rendah: 55,6%
- Distribusi status imunisasi anak secara keseluruhan:
Anak dengan imunisasi lengkap: 48,9%
Anak dengan imunisasi tidak lengkap: 51,1%
- Distribusi bersyarat (conditional probability):
Dari ibu berpendidikan tinggi, 60% anak mendapat imunisasi lengkap.
Dari ibu berpendidikan rendah, hanya 40% anak yang mendapat imunisasi lengkap.
- Ukuran asosiasi:
Risk Difference (RD) = 0,2 ➤ Terdapat selisih 20 poin persentase dalam cakupan imunisasi lengkap antara anak dari ibu berpendidikan tinggi dan rendah.
Relative Risk (RR) = 1,5 ➤ Anak dari ibu berpendidikan tinggi memiliki risiko 1,5 kali lebih besar untuk mendapatkan imunisasi lengkap dibanding anak dari ibu berpendidikan rendah.
Odds Ratio (OR) = 2,25 ➤ Peluang anak untuk mendapatkan imunisasi lengkap lebih dari dua kali lipat jika ibunya memiliki pendidikan tinggi, dibandingkan dengan ibu berpendidikan rendah.
Kesimpulan Interpretatif Temuan ini menunjukkan bahwa pendidikan ibu berperan signifikan dalam menentukan status imunisasi anak. Proporsi anak yang menerima imunisasi lengkap secara nyata lebih tinggi di antara ibu berpendidikan tinggi. Ukuran asosiasi yang kuat (RR = 1,5 dan OR = 2,25) memperkuat bukti bahwa peningkatan pendidikan perempuan dapat menjadi salah satu strategi kunci dalam meningkatkan cakupan imunisasi dasar di masyarakat.
Implikasi praktisnya, hasil ini mendukung perlunya:
Kampanye imunisasi yang mempertimbangkan faktor pendidikan, terutama di komunitas dengan latar belakang pendidikan rendah.
Intervensi kebijakan yang bersifat edukatif dan peningkatan literasi kesehatan bagi ibu dan keluarga sebagai bagian dari upaya preventif dan promotif di sektor kesehatan masyarakat.
Dengan demikian, analisis tabel kontingensi sederhana dapat memberikan bukti kuantitatif yang kuat untuk mendukung perumusan kebijakan publik yang lebih inklusif dan berbasis data.
6 Inferensi Tabel Kontingensi Dua Arah
Bab ini membahas berbagai teknik inferensi statistik yang digunakan dalam analisis tabel kontingensi dua arah. Tujuannya adalah untuk menguji apakah terdapat hubungan yang signifikan antara dua variabel kategorik. Contoh yang digunakan berasal dari data SDKI 2017, yakni hubungan antara tingkat pendidikan ibu (tinggi vs rendah) dan status imunisasi dasar lengkap anak usia 12–23 bulan (lengkap vs tidak lengkap).
6.1 Estimasi Parameter
Estimasi parameter dilakukan untuk menyimpulkan karakteristik populasi berdasarkan sampel. Dalam konteks tabel 2×2, kita menaksir proporsi tiap kelompok dan ukuran asosiasi seperti RD (Risk Difference), RR (Relative Risk), dan OR (Odds Ratio).
6.1.1 Estimasi Titik
Hasil estimasi titik menunjukkan:
Proporsi anak dengan imunisasi lengkap pada ibu berpendidikan tinggi: 0,60
Proporsi anak dengan imunisasi lengkap pada ibu berpendidikan rendah: 0,40
Selisih proporsi ini adalah dasar untuk menghitung ukuran asosiasi antara kedua variabel.
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.2288055Interpretasi: Terdapat selisih proporsi imunisasi lengkap sebesar 20% (RD = 0,20) antara dua kelompok ibu. Dengan tingkat kepercayaan 95%, kita yakin bahwa selisih tersebut berada di antara 16,6% hingga 23,4%, yang menunjukkan efek yang signifikan secara praktis dan statistik.
6.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
Analisis tabel kontingensi tiga arah memungkinkan kita mengeksplorasi hubungan antara dua variabel kategorik dengan mempertimbangkan satu variabel ketiga sebagai kontrol atau stratifikasi. Ini sangat penting dalam penelitian sosial dan kesehatan untuk mendeteksi pengaruh perancu, interaksi, atau perbedaan efek di tiap strata.
7.1 Tabel Parsial dan Marginal
Kita dapat menyusun tabel parsial untuk melihat hubungan antara dua variabel (misalnya gizi dan diare) pada masing-masing tingkat pendidikan. Sementara itu, tabel marginal menggambarkan keseluruhan data tanpa memperhitungkan strata pendidikan.
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 dalam tabel tiga arah dapat dihitung dalam tiga cara:
Peluang Bersama: probabilitas kombinasi nilai dari ketiga variabel.
Peluang Bersama: \[ P(X = x, Y = y, Z = z) = \frac{n_{xyz}}{N} \]
Peluang Marginal: Peluang Marginal: probabilitas total untuk salah satu variabel, mengabaikan variabel lainnya.
\[ 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
Peluang Bersyarat: probabilitas suatu kejadian, dengan syarat variabel kontrol tertentu.
\[ 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
Setiap strata pendidikan menghasilkan ukuran asosiasi seperti:
Risk Difference (RD): selisih risiko diare antara anak gizi kurang dan gizi baik.
Relative Risk (RR): perbandingan risiko diare antara dua kelompok gizi.
Odds Ratio (OR): perbandingan odds diare antar dua kelompok.
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)## Loading required package: Matrix## Loading required package: metadat## 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: 0x11df36608>
## <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
Studi Kasus ini menyajikan penerapan analisis data kategorik menggunakan tabel kontingensi 2 × 2 untuk mengevaluasi hubungan antara kebiasaan merokok dan kejadian kanker. Data yang digunakan bersifat observasional dan dikumpulkan dari populasi berdasarkan hasil diagnosis dan status kebiasaan merokok.
10.1. Tabel Data
Data berikut menunjukkan jumlah penderita kanker berdasarkan status merokok:
# 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 27Tabel ini menunjukkan bahwa dari 1269 perokok, terdapat 647 yang menderita kanker dan 622 yang tidak. Sementara itu, dari 29 non-perokok, hanya 2 orang yang menderita kanker dan 27 tidak. ### 10.2 Distribusi Probabilitas
Distribusi peluang bersama, marginal, dan bersyarat dihitung untuk melihat pola hubungan antara merokok dan kanker.
# 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.06896552Interpretasi:
Peluang bersama menunjukkan probabilitas seseorang berada pada kombinasi tertentu (misalnya perokok dan menderita kanker).
Peluang bersyarat P(Kanker | Perokok) ≈ 0.51, sedangkan P(Kanker | Non-perokok) ≈ 0.07, menunjukkan adanya perbedaan besar dalam proporsi kanker antara perokok dan non-perokok.
10.3 Ukuran Asosiasi
Untuk menilai kekuatan hubungan antara merokok dan kejadian kanker, digunakan ukuran asosiasi: Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR).
# 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.04Hasil:
Risk Difference (RD): 0.4409 Artinya, perokok memiliki selisih risiko kanker sekitar 44% lebih tinggi dibandingkan non-perokok.
Relative Risk (RR): 7.39 Peluang terkena kanker pada perokok sekitar 7 kali lebih tinggi dibandingkan non-perokok.
Odds Ratio (OR): 14.04 Perokok memiliki kemungkinan terkena kanker lebih dari 14 kali dibandingkan non-perokok, berdasarkan rasio odds.
10.4 Kesimpulan
Analisis ini menunjukkan adanya hubungan kuat antara status merokok dan kejadian kanker. Ketiga ukuran asosiasi (RD, RR, dan OR) memberikan bukti bahwa risiko kanker meningkat signifikan pada kelompok perokok dibandingkan non-perokok.
11 Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio
Regresi logistik digunakan untuk memodelkan hubungan antara variabel respons biner (misalnya: sukses/gagal) dengan satu atau lebih variabel prediktor. Variabel prediktor dapat berupa:
Nominal: Tidak memiliki urutan logis. Contoh: jenis pekerjaan (kantoran, lapangan).
Ordinal: Memiliki urutan logis, tetapi jarak antar kategori belum tentu sama. Contoh: kepuasan (rendah < sedang < tinggi).
Rasio: Variabel numerik kontinu. Contoh: jumlah jam kerja per minggu.
11.1 Simulasi Data
set.seed(1234)
n <- 600
pekerjaan <- sample(c("Kantor", "Lapangan"), n, replace = TRUE)
kepuasan <- sample(c("Rendah", "Sedang", "Tinggi"), n, replace = TRUE, prob = c(0.3, 0.4, 0.3))
jam_kerja <- round(rnorm(n, mean = 45, sd = 10), 1)
logit_p <- -1 + 0.7 * (pekerjaan == "Kantor") + 0.5 * as.numeric(factor(kepuasan, ordered = TRUE)) + 0.02 * jam_kerja
p <- 1 / (1 + exp(-logit_p))
hasil <- rbinom(n, 1, p)
simulasi <- tibble(hasil, pekerjaan, kepuasan, jam_kerja)
head(simulasi)## # A tibble: 6 × 4
## hasil pekerjaan kepuasan jam_kerja
## <int> <chr> <chr> <dbl>
## 1 1 Lapangan Sedang 48.1
## 2 1 Lapangan Sedang 51.1
## 3 0 Lapangan Sedang 28.1
## 4 1 Lapangan Tinggi 52.8
## 5 1 Kantor Tinggi 45.1
## 6 0 Lapangan Sedang 43.211.2 Eksplorasi Data
simulasi %>%
group_by(hasil) %>%
summarise(Jumlah = n(), Rata2_Jam = mean(jam_kerja))## # A tibble: 2 × 3
## hasil Jumlah Rata2_Jam
## <int> <int> <dbl>
## 1 0 168 43.6
## 2 1 432 45.0Interpretasi: Kelompok dengan keberhasilan (hasil = 1) memiliki rata-rata jam kerja lebih tinggi.
11.3 Perlakuan Variabel Ordinal
11.3.1 Sebagai Nominal (Dummy)
simulasi_nominal <- simulasi %>%
mutate(kepuasan = factor(kepuasan, levels = c("Rendah", "Sedang", "Tinggi")))
model_nominal <- glm(hasil ~ pekerjaan + kepuasan + jam_kerja, data = simulasi_nominal, family = binomial)
summary(model_nominal)##
## Call:
## glm(formula = hasil ~ pekerjaan + kepuasan + jam_kerja, family = binomial,
## data = simulasi_nominal)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.042833 0.488651 0.088 0.93015
## pekerjaanLapangan -0.918020 0.193746 -4.738 2.16e-06 ***
## kepuasanSedang 0.692899 0.218670 3.169 0.00153 **
## kepuasanTinggi 1.136991 0.255425 4.451 8.53e-06 ***
## jam_kerja 0.018331 0.009987 1.835 0.06643 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 711.54 on 599 degrees of freedom
## Residual deviance: 665.21 on 595 degrees of freedom
## AIC: 675.21
##
## Number of Fisher Scoring iterations: 4Interpretasi: Kategori Sedang dan Tinggi dibandingkan terhadap baseline Rendah.
Koefisien positif pada Tinggi menunjukkan peningkatan log-odds keberhasilan yang signifikan.
jam_kerja juga memiliki pengaruh positif.
11.3.2 Sebagai Rasio (Numeric Rank)
simulasi_numeric <- simulasi %>%
mutate(kepuasan_num = case_when(
kepuasan == "Rendah" ~ 1,
kepuasan == "Sedang" ~ 2,
kepuasan == "Tinggi" ~ 3
))
model_numeric <- glm(hasil ~ pekerjaan + kepuasan_num + jam_kerja, data = simulasi_numeric, family = binomial)
summary(model_numeric)##
## Call:
## glm(formula = hasil ~ pekerjaan + kepuasan_num + jam_kerja, family = binomial,
## data = simulasi_numeric)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.485030 0.533283 -0.910 0.3631
## pekerjaanLapangan -0.920873 0.193502 -4.759 1.95e-06 ***
## kepuasan_num 0.578571 0.128265 4.511 6.46e-06 ***
## jam_kerja 0.018230 0.009987 1.825 0.0679 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 711.54 on 599 degrees of freedom
## Residual deviance: 665.63 on 596 degrees of freedom
## AIC: 673.63
##
## Number of Fisher Scoring iterations: 4Interpretasi: - Koefisien kepuasan_num menunjukkan bahwa setiap peningkatan satu tingkat kepuasan meningkatkan peluang keberhasilan. - Lebih sederhana, cocok jika jarak antar kategori dianggap seimbang.
11.4 Perbandingan Model
list(
AIC_Nominal = AIC(model_nominal),
AIC_Numeric = AIC(model_numeric)
)## $AIC_Nominal
## [1] 675.2102
##
## $AIC_Numeric
## [1] 673.6273Kesimpulan: Model dengan AIC lebih rendah dipilih. Jika perbedaan kecil, model ordinal sebagai numeric lebih mudah diinterpretasikan.
11.5 Visualisasi Prediksi
library(ggplot2)
simulasi_nominal <- simulasi_nominal %>%
mutate(prediksi = predict(model_nominal, type = "response"))
ggplot(simulasi_nominal, aes(x = jam_kerja, y = prediksi, color = kepuasan)) +
geom_point(alpha = 0.5) +
labs(title = "Prediksi Probabilitas: Kepuasan sebagai Nominal", y = "Probabilitas Keberhasilan", x = "Jam Kerja") +
theme_minimal()
simulasi_numeric <- simulasi_numeric %>%
mutate(prediksi = predict(model_numeric, type = "response"))
ggplot(simulasi_numeric, aes(x = jam_kerja, y = prediksi, color = as.factor(kepuasan_num))) +
geom_point(alpha = 0.5) +
labs(title = "Prediksi Probabilitas: Kepuasan sebagai Numeric", y = "Probabilitas Keberhasilan", x = "Jam Kerja") +
theme_minimal()
11.6 Ringkasan dan Kesimpulan
Pekerjaan tipe kantor meningkatkan peluang keberhasilan.
Kepuasan kerja berpengaruh signifikan terhadap keberhasilan.
Jam kerja lebih banyak cenderung meningkatkan peluang keberhasilan.
Perlakuan variabel ordinal bisa disesuaikan dengan konteks dan asumsi data.
12. Pemilihan Model Regresi Logistik dan Evaluasi
12.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) ## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, varlibrary(DescTools)
set.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.519360312.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.265712.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.778912.4 Evaluasi Pseudo R-Squared
PseudoR2(step_both, which = c("CoxSnell", "Nagelkerke", "McFadden"))## CoxSnell Nagelkerke McFadden
## 0.2390347 0.3187198 0.197057812.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] "step_forward"# 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.7812.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.
12.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.
12.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.
12.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.
12.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## Loading required package: rlang# 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.
12.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.258629212.11.1 Perhitungan Otomatis dengan Paket Tambahan
# pscl
if (!require(pscl)) install.packages("pscl"); library(pscl)## Loading required package: pscl## 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.
13. 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
13.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.0627056613.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.
13.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.
13.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 613.5 Estimasi Model Multinomial Logistik
if (!require(nnet)) install.packages("nnet"); library(nnet)## Loading required package: nnetmodel_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.37613.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
13.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.
13.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.
13.8.1 Contoh Kasus 2 di R (Dataset Iris)
# Load semua package yang dibutuhkan
library(nnet)
library(dplyr)
library(ggplot2)
library(MASS) # Untuk regresi logistik ordinal13.8.1 Contoh Kasus 2: Regresi Logistik Multinomial dengan Dataset Iris
# Dataset Iris: prediksi spesies berdasarkan panjang dan lebar petal
data(iris)
# Ubah level referensi menjadi 'setosa'
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 iterations# Ringkasan model
summary(model)## Call:
## multinom(formula = Species ~ Petal.Length + Petal.Width, data = iris)
##
## Coefficients:
## (Intercept) Petal.Length Petal.Width
## versicolor -22.79944 6.92122 7.878496
## virginica -67.82521 12.64721 18.261016
##
## Std. Errors:
## (Intercept) Petal.Length Petal.Width
## versicolor 44.3859 37.58715 81.00888
## virginica 46.3939 37.65702 81.09482
##
## Residual Deviance: 20.57901
## AIC: 32.57901# Hitung 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.8218# Prediksi spesies
iris$predicted <- predict(model, newdata = iris)
# Tabel perbandingan prediksi vs 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
ggplot(iris, aes(x = Petal.Length, y = Petal.Width, color = predicted)) +
geom_point(size = 2) +
labs(title = "Prediksi Spesies Iris",
x = "Petal Length",
y = "Petal Width",
color = "Prediksi") +
theme_minimal()
14. Regresi Logistik Ordinal
Regresi logistik ordinal digunakan saat variabel respon bersifat urut, seperti tingkat kepuasan atau level risiko. Ini berbeda dari: - Regresi logistik biner: 2 kategori - Regresi logistik multinomial: >2 kategori tanpa urutan
14.1 Konsep Cumulative Logit Model
Model logit kumulatif:
log(P(Y ≤ j) / P(Y > j)) = α_j + βx- αj: intercept untuk setiap batas kategori
- β: koefisien regresi (sama untuk semua batas)
Untuk c kategori, akan ada (c−1) fungsi logit.
14.2 Interpretasi Koefisien
- β > 0 → meningkatkan kemungkinan masuk kategori rendah/sedang
- β < 0 → meningkatkan kemungkinan masuk kategori tinggi
Odd ratio: exp(β)
14.3 Contoh Data: Penilaian Risiko Kesehatan
set.seed(789)
n <- 200
# Skor pola hidup sehat (semakin tinggi = semakin sehat)
health_score <- round(runif(n, 1, 10))
# Simulasi level risiko berdasarkan skor + error acak
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)
# Buat data frame
df <- data.frame(risk_level, health_score)
# Lihat data awal
head(df)## risk_level health_score
## 1 Tinggi 7
## 2 Rendah 2
## 3 Rendah 1
## 4 Sedang 6
## 5 Sedang 5
## 6 Rendah 114.4 Estimasi Model Ordinal dengan polr()
# Model logistik ordinal
model_ord <- polr(risk_level ~ health_score, data = df, Hess = TRUE)
# Ringkasan model
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.988714.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
14.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.0136127314.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.
14.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.
14.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
14.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
15 Log Linear Model
15.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.
15.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.
15.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.
15.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
15.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.
15.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.
15.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 116 Ringkasan Model Log-Linear 2 Arah
16.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}
16.2 Perbedaan Log-Linear vs Regresi Logistik
- Log-linear : memodelkan frekuensi, tanpa variabel respon.
- Regresi logistik : memodelkan probabilitas suatu kejadian (variabel respon eksplisit).
16.3 Estimasi Parameter Log-Linear 2 Arah
Gunakan constraint sum-to-zero untuk identifikasi parameter:
math +* = 0, +* = 0, * ^{AB}{ij} = 0
16.4 Analisis Data Tabel Kontingensi 2x2
Misal data hasil survei:
| Konsumsi Kafein | Insomnia | Tidak Insomnia |
|---|---|---|
| Ya | 34 | 16 |
| Tidak | 12 | 38 |
16.5 Bentuk Model Log-Linear
Struktur model:
math log() =* + ^A_i + ^B_j + ^{AB}{ij}
dengan constraint sum-to-zero.
16.6 Estimasi Parameter (Manual)
Menghitung parameter dengan rumus log-frekuensi dan constraint sum-to-zero.
16.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)16.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: 1.7764e-15 on 0 degrees of freedom
## AIC: 27.807
##
## Number of Fisher Scoring iterations: 316.9 Interpretasi Parameter
- Intercept: log frekuensi dasar.
- Efek utama: pengaruh masing-masing kategori.
- Interaksi signifikan: hubungan antara kafein dan insomnia.
16.10 Tabel Kontingensi 2x3
Data survei:
| Gender | Kurus | Normal | Gemuk |
|---|---|---|---|
| Laki-laki | 15 | 23 | 12 |
| Perempuan | 18 | 20 | 17 |
16.11 Bentuk Model Log-Linear 2x3
Model ekspansi dari 2x2:
math log() =* + ^A_i + ^B_j + ^{AB}{ij}
16.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: 3.7748e-15 on 0 degrees of freedom
## AIC: 40.135
##
## Number of Fisher Scoring iterations: 316.13 Interpretasi Model 2x3
- Efek utama: perbedaan proporsi BMI antar gender.
- Interaksi: apakah distribusi BMI berbeda secara signifikan antar gender.
17 Model Log Linear Tiga Arah
17.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)
17.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
17.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 5017.4 Analisis Log-Linear untuk Tabel Tiga Arah
library(MASS)17.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.5Interpretasi : 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.
18 Uji Model Interaksi Tiga Arah
18.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.
18.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: 418.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.
18.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.
18.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.
18.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
19 Pemilihan Model Terbaik
19.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 |
19.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 |
19.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.
20 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.
21 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 NA21.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.
22 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 5022.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.
23 Contoh Data Kasus
23.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.
- 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
- Wickham, H., & Grolemund, G. (2016). R for Data Science. O’Reilly Media. https://r4ds.had.co.nz
- Dobson, A. J., & Barnett, A. G. (2018). An Introduction to Generalized Linear Models (4th ed.). CRC Press.
- James, G., Witten, D., Hastie, T., & Tibshirani, R. (2021). An Introduction to Statistical Learning with Applications in R (2nd ed.). Springer.
- Kuhn, M., & Johnson, K. (2013). Applied Predictive Modeling. Springer.
- Friendly, M., & Meyer, D. (2016). Discrete Data Analysis with R: Visualization and Modeling Techniques for Categorical and Count Data. CRC Press.
- R Core Team. (2024). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/